話說「交集」是個常常被掛在嘴邊的語句,像我和一口巾工讀生們往往都沒有交集,他們不了解我的明白,我也不了解他們的明白啊。交集在幾何上可以用下圖來簡單表示,但是在 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