話說「交集」是個常常被掛在嘴邊的語句，像我和一口巾工讀生們往往都沒有交集，他們不了解我的明白，我也不了解他們的明白啊。交集在幾何上可以用下圖來簡單表示，但是在 GIS 計算交集上，除了空間上的交集外，還有的資料表的交集(ex: join)，這個我就不多花篇幅談。
以下簡單介紹用空間資料庫的概念來實作多邊形的面積以及交集
一、準備材料：
- PostgreSQL (版本 >= 9.1)
- PostGIS 支援 PostgreSQL 的空間模組 (版本 >= 2.0)
- QGIS 展圖用
二、步驟
上圖是八色鳥在的分布範圍（圖層： range_fairy_pitta，橘色虛線淺綠色底的多邊形區域，簡稱 sp），深綠色是世界的陸域地圖（圖層：map_world，資料屬性表中有國家名稱(cntry_name)以及大陸名稱(continent)，簡稱 world），淺橘色則是 sp 和 world 交集的區域。當然我知道做這個很簡單，從 ArcMap 中的工具箱選擇交集工具就可以做出來了，或是用 QGIS 的向量外掛(fTools)的 Geoprocessing Tool/Intersect 即可做出來，但是 GUI 用起來就是有不踏實感 XD 所以我們用空間資料庫來實作！
我們需要使用到 postgis 的函式有 geometry ST_Intersection(geometry_a, geometry_b) 、 boolean ST_Intersects(geometry_a, geometry_b)，以及 float ST_Area(geometry::geography) 這三個函式。先說明 PostGIS 函式的一些規則，假設有個 return_type FUNCTION(input parameter) 函式，函式名稱就叫做 FUNCTION，而 return_type 表示執行這個函式會回傳的類型，有可能是浮點數、整數或布林值(boolean)，而 input parameter 則是這個函式運上所需要輸入的參數，例如幾何欄位(geometry)，其餘詳細說明請參閱 PostGIS 文件中的地理資料類型章節。
在進行計算交集前，請先建立空間資料表，你也可以直接用 ESRI Shapefile 匯入 PostgreSQL 資料庫中，請參閱將 ESRI Shapefile 轉成 postgis 格式匯入
接下來我們如果要取交集的話，可以先用 ST_Intersects(geometry_a, geometry_b)，SQL 可以這樣表示：
SELECT * FROM spatial_tab_a a, spatial_tab_b b WHERE ST_Intersects(a.geometry_column, b.geometry_column);
詳解：邏輯上可以這樣想，當 ST_Intersects(a, b) 成立時，將空間資料表 a (spatial_tab_a) 和 b (spatial_tab_b) 交集的部份選出來，因此 a, b 兩個空間資料表所有的欄位都會被選取並回傳。
但是我們想要把交集的部份另外存成一個空間資料表，這應該怎麼做呢？因為 ST_Intersects() 只會回傳布林值，但是我們只希望有交集的多邊形幾何值，此時我們就必須要使用 ST_Intersection(geometry_a, geometry_b) 來計算交集多邊形區域的幾何值，所以我們稍微改一下剛剛的 SQL 語法：
SELECT ST_Intersection(a.geom, b.geom) geom FROM spatial_tab_a a, spatial_tab_b b WHERE ST_Intersects(a.geometry_column, b.geometry_column);
詳解：a,b 兩個圖層計算交集，並把交集的幾何值欄位(geom)計算出來，如下圖範例
接下來我們來算交集過後的面積，使用 ST_Area(geometry::geography)，架設這個已經是投影過後的圖層，我們就可以直接算面積，或是用 ST_Transform() 來轉換投影系統，但如果是 WGS84 地理座標系統，我們可以用 ST_Area(geometry::geography) 直接算面積。
SELECT ST_Intersection(a.geom, b.geom) geom, ST_Area(ST_Intersection(a.geom, b.geom)::geography) sqmeter FROM spatial_tab_a a, spatial_tab_b b WHERE ST_Intersects(a.geometry_column, b.geometry_column);
詳解：將 ST_Intersection(a.geom, b.geom) 所計算出來的幾何值丟給 ST_Area() 計算面積（平方公尺，另外重新命名欄位名稱為 sqmeter），如下圖：
但是這樣我們除了幾何上的交集，並沒有其他的屬性資料(attribute data)，也沒有空間索引(spatial index)及主鍵(primary key)，所以我們再做更細緻的處理，加入主鍵並且是具有空間索引的主鍵(使用 GIST(Generalized Search Tree Index)，在龐大的空間資料表中可以加速搜尋的速度）：
先建立一個索引序列 :
CREATE SEQUENCE sp_serial;
之後我們會用它自動建立一個 gid 主鍵，丟給他從 1 開始的序列，使用 nextval(‘sp_serial’::regclass) 來輸入 gid 值：
SELECT nextval('sp_serial'::regclass) gid, ST_Intersection(a.geom, b.geom) geom, ST_Area(ST_Intersection(a.geom, b.geom)::geography) sqmeter FROM spatial_tab_a a, spatial_tab_b b WHERE ST_Intersects(a.geometry_column, b.geometry_column);
接下來我們想要在這個交集的空間資料表中顯示 a 表格的某個欄位(column_a)，b 表格的某個欄位(column_b)：
SELECT nextval('sp_serial'::regclass) gid, a.column_a, b.column_b, ST_Intersection(a.geom, b.geom) geom, ST_Area(ST_Intersection(a.geom, b.geom)::geography) sqmeter FROM spatial_tab_a a, spatial_tab_b b WHERE ST_Intersects(a.geometry_column, b.geometry_column);
如下圖：
另外存成一個新的空間資料表，使用 CREATE TABLE AS (SELECT * FROM )：
CREATE TABLE ab_intersection AS (SELECT nextval('sp_serial'::regclass) gid, a.column_a, b.column_b, ST_Intersection(a.geom, b.geom) geom, ST_Area(ST_Intersection(a.geom, b.geom)::geography) sqmeter FROM spatial_tab_a a, spatial_tab_b b WHERE ST_Intersects(a.geometry_column, b.geometry_column));
將 gid 設成主鍵，並且建立 GIST:
ALTER TABLE ab_intersection ADD PRIMARY KEY (gid);
CREATE INDEX ab_intersection_geom_gist ON ab_intersection USING GIST ("geom");
大致上就是這樣的處理方式，如果要進行批次匯入交集，也可以使用 scripting 語言如 python (pyscopg) 或 php 來處理
2012/12/25 14:45 更新：
PostGIS 建議用 AddGeometryColumn() 來增加 geometry 欄位，是因為如果要符合 OpenGIS metadata 規範，必須要加入 SID，如果沒有使用 AddGeometryColumn() 則會產生 -1 的值，因此可以用手動加入的方式，以下使用 WGS84 (EPSG 4326) 手動在 geometry 欄位中加入 sid
ST_Transform(geometry, 4326)::geometry(Geometry, 4326)
參考資料
http://postgis.refractions.net/documentation/manual-1.3/ch03.html#id2570755
http://postgis.refractions.net/docs/using_postgis_dbmanagement.html#Manual_Register_Spatial_Column
