利用空間資料庫計算多邊形交集以及計算面積

話說「交集」是個常常被掛在嘴邊的語句,像我和一口巾工讀生們往往都沒有交集,他們不了解我的明白,我也不了解他們的明白啊。交集在幾何上可以用下圖來簡單表示,但是在 GIS 計算交集上,除了空間上的交集外,還有的資料表的交集(ex: join),這個我就不多花篇幅談。Image

以下簡單介紹用空間資料庫的概念來實作多邊形的面積以及交集

一、準備材料:

二、步驟

qgis_layout

上圖是八色鳥在的分布範圍(圖層: 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)計算出來,如下圖範例

inter_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),如下圖:

inter_geom2

但是這樣我們除了幾何上的交集,並沒有其他的屬性資料(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);

如下圖:

inter_geom3

另外存成一個新的空間資料表,使用 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

Comments

comments powered by Disqus