用 SQL 分析 GIS
MS SQL Geoprocess
Geography & Geometry
這兩個都是地理資訊上常見名詞，也是常常搞混的兩個單詞，簡單來說當我們要計算距離、周長及面積時，空間不會是一個平面，因此隨著點跟點之間亦或者是物件之間的距離越來越遠，所產生的誤差就越來越大，因此我們聰明的地理學家麥卡托（Gerardus Mercator），就透過他天才的理論，將空間壓成了平面，透過一連串的演算來簡化誤差的問題。
透過下面兩張圖得比較，對於兩者的差異已經能有初步的認知，那麼在資料上呢?廣泛用於資料處理的經緯度座標，如:121.3345678,23.478587 這種資料是屬於 Geometry 還是 Geography 呢? 那二度分帶座標，如:248172,2652130 又是屬於哪種呢? 答: 經緯度座標屬於 Geography; 二度分帶座標屬於 Geometry
＊ 基於上面這個概念，因此當我們在使用 MS SQL Spatial 時，要特別注意將取得的經緯度座標轉換為二度分帶座標，再來計算距離、周長及面積。
一、資料格式 (Data Format)
使用 Geomertry 資料格式，支援以下兩種格式
- WKB(well-known binary)
- WKT(well-known text)
上述兩種格式，皆為開放式地理空間協會 (OGC) 已知表示字串形式，主要是用來表示平面幾何執行個體。因此，此兩種資料格式將不會包含此執行個體所夾帶的任何 Z (高度) 或 M (測量) 值。
WKB 格式範例
| 礦場名稱 | GeomertryBinaryResult | |
|---|---|---|
| 1 | 白沙屯 8 號探井 | 0x000000000105080000003AD6AE50B2315E4094FD868C889938406C67AD8298315E409A75405DB39838403B623CA07C2F5E400E3568284A953840C596F5C52E2F5E404DD71F3F0B973840B3253B40E3305E40A1AC03E0E899384078AF64C457315E4009BB56BC5299384050C7519A86315E40C97F9461109A38403AD6AE50B2315E4094FD868C8899384000000000B8804D3F00000000587D4D3F00000000746F4D3F0000000094764D3F000000003E824D3F00000000DC7F4D3F00000000DE824D3F00000000B8804D3F01000000020000000001000000FFFFFFFF0000000003 |
WKT 格式範例
依照不同的資料樣態可分為 POINT,LINESTRING,POLYGON
| 礦場名稱 | GeomertryWKTResult | |
|---|---|---|
| 1 | 白沙屯 8 號探井 | POLYGON ((120.77650849412831 24.599739821392703, 120.77493349966034 24.596486881482271, 120.74198156258437 24.583162808839056, 120.73722981437011 24.590015359185873, 120.763870294343 24.601209641358427, 120.77098188240586 24.598918696569886, 120.7738405035609 24.6018124568507, 120.77650849412831 24.599739821392703)) |
二、建置工具準備
環境準備
- SQL SERVER 2008 以上版本(含)，需支援 geomertry 資料屬性。
- .NET Framework 3.5 以上版本
- IDE VS (Visual Studio) , VSC (Visual Studio Code)
輔助工具
- SHP2SQL(OpenSource)，GIS 檔案匯入資料庫。
- QGIS(OpenSource)，各式 GIS 檔案格式轉換、資料展示、空間分析及應用。
三、資料比較方法 (Support Method)
透過 MS SQL SERVER 提供的 OGC 幾何項目方法，可以進行效率及便捷的空間幾何與平面運算，由於資料演算及邏輯部分已由微軟這邊支援，其中學術演算及原理就不詳細說明，這邊僅講求如何使用微軟支援的方法說明。
此方法回傳的資料類型可分為空間幾何與位元(Boolean)。下段文字將針對常用的空間分析演算進行說明與範例實作，以下所述的方法不論是 Geography 或是 Geometry 皆有支援，兩者的差別在於對空間資料的詮釋及計算結果的差異。共通的部分為皆須要轉換為 ST Text Format。
起手式 ST Text Format
-- obj: 欲操作的物件 -- srs: 坐標系統:4326,3826,3857,0,etc..... geometry::STGeomFromText(obj,srs) geography::STGeomFromText(obj,srs)
二進位字元資料轉換
- STAsText() 主要用於將透過輔助工具將 GIS Shape File 匯入資料庫的 Binary geometry 轉換為可讀的文字模式。
| 礦場名稱 | Shape | WKTFormat |
|---|---|---|
| 新寮山石礦 | 0xF20E00000107... | POLYGON ((326164.1672 2715456.7046000008, 326163.97360000014 2715457.4999, 326010.38920000009 2716088.5636, 325616.17520000041 2716014.5116000008, 325660.25050000008 2715358.0406, 325660.2592000002 2715357.9105999991, 326164.1672 2715456.7046000008)) |
SELECT top 1 [礦場名稱],Shape, Shape.STAsText() as WKTFormat FROM [dbo].[MIN_MineArea_107]
關於單位
在空間幾何上的資料一定有其單位，舉凡距離、周長和面積因應不同座標系統都有其換算單位，要如何得知 SQL Geomertry 計算出來的單位呢? 答: 透過資料參照的座標系統與以下此段查詢結果即可得知。
SELECT spatial_reference_id , well_known_text , unit_of_measure , unit_conversion_factor FROM sys.spatial_reference_systems
一般廣泛使用的座標系統為 WGS84 其可對應的資訊為 EGPS:4326。
| spatial_reference_id | well_known_text | unit_of_measure | |
|---|---|---|---|
| 1 | 4326 | GEOGCS["WGS 84", DATUM["World Geodetic System 1984", ELLIPSOID["WGS 84", 6378137, 298.257223563]], PRIMEM["Greenwich", 0], UNIT["Degree", 0.0174532925199433]] | metre |
距離計算
- STLength() 可用於線段長度計算及周長計算，當輸入為兩點則計算線段，若為兩點以上則計算為周長，只有一個點則回傳 0。
-- 線段長度計算 Declare @l geometry set @l = geometry::STGeomFromText('LINESTRING(326164.1672 2715456.7046000008, 326163.97360000014 2715457.4999)',0) select @l.STLength(); -- 結果驗算 select sqrt(SQUARE(326163.97360000014-326164.1672)+SQUARE(2715457.4999-2715456.7046000008)) -- 周長計算 SELECT top 1 [礦場名稱],geometry::STGeomFromText(Shape.STAsText(),4326).STLength() as distance FROM [dbo].[MIN_MineArea_107]
物件之間距離計算
- STDistance()
用於計算 A 物件(點、線、面)到 B 物件(點、線、面) 之間的距離，一旦接觸即回傳。
如果 other_geometry 是空的集合，
STDistance()會傳回 Null，如果是兩條交疊的線段則回傳 0。
-- A 線段 到 B 點 的距離 Declare @l2 geometry Declare @p geometry set @l2 = geometry::STGeomFromText('LINESTRING(326164.1672 2715456.7046000008, 326163.97360000014 2715457.4999)',0) set @p = geometry::STGeomFromText('POINT(328564.1672 2716456.7046000008)',0) SELECT @l2.STDistance(@p) as distance;
-- B 點 到 A 區塊的 的距離 Declare @p2 geometry set @p2 = geometry::STGeomFromText('POINT(328564.1672 2716456.7046000008)',0) SELECT top 1 geometry::STGeomFromText(Shape.STAsText(),0).STDistance(@p2) as distance FROM [dbo].[MIN_MineArea_107]
-- 路線偏移計算運用 Declare @l3 geometry Declare @l4 geometry set @l3 = geometry::STGeomFromText('LINESTRING(326164.1672 2715456.7046000008, 326163.97360000014 2715457.4999)',0) set @l4 = geometry::STGeomFromText('LINESTRING(326164.1672 2715456.7046000008, 326163.97360000014 2715457.4999)',0) SELECT @l3.STDistance(@l4) as distance;
面積計算
- STArea()
如果執行個體不一個封閉多變形，或它是空的，
STArea()會傳回 0。 如果 geometry 執行個體尚未初始化，STArea()會傳回 Null。
-- Geography SELECT top 1 geom, geom.STAsText() as WKTFormat,[AREA_LAW],geometry::STGeomFromText(geom.STAsText(),4326).STArea()*10000,geography::STGeomFromText(geom.STAsText(),4326).STArea()/1000000 FROM [dbo].[WSEPA_WQ]
-- Geometry SELECT top 1 [礦場名稱],Shape, Shape.STAsText() as WKTFormat,geometry::STGeomFromText(Shape.STAsText(),0).STArea()/10000,[面積] FROM [dbo].[MIN_MineArea_107]
交集
- STIntersection() 回傳兩個物件交集的部分
- STIntersects() 用於邏輯運算，當兩個物件有交集回傳 1 ，沒有交集回傳 0。
-- 面 DECLARE @a geometry;
DECLARE @a2 geometry;
SET @a = geometry::STGeomFromText('POLYGON((326164.1672 2715456.7046000008, 326163.97360000014 2715457.4999, 327133.97360000014 2756531.4999, 326164.1672 2715456.7046000008))', 0);
SET @a2 = geometry::STGeomFromText('POLYGON ((326164 2715456.7046000008, 329163.97360000014 2715457.4999, 326010.38920000009 2716088.5636, 325616.17520000041 2716014.5116000008, 325660.25050000008 2715358.0406, 325660.2592000002 2715357.9105999991, 326164 2715456.7046000008))', 0);
-- 面交集的範圍 SELECT @a.STIntersection(@a2) as IntersectionResult;
-- 面交集判斷 --SET @a4 = geometry::STGeomFromText('POLYGON ((1 1, 2.0000000000000036 1.0000000000000071, 2 2, 1.0000000000000071 2.0000000000000036, 1 1))', 0);
SELECT @a.STIntersects(@a2) as IntersectionResult;
-- 線 DECLARE @l5 geometry;
DECLARE @l6 geometry;
SET @l5 = geometry::STGeomFromText('LINESTRING(0 0, 5 5)', 0);
SET @l6 = geometry::STGeomFromText('LINESTRING(0 5, 5 0)', 0);
-- 線交集的範圍 SELECT @l5.STIntersection(@l6).STAsText() as IntersectionResult;
-- 線交集判斷 --SET @l6 = geometry::STGeomFromText('POLYGON ((1 1, 2.0000000000000036 1.0000000000000071, 2 2, 1.0000000000000071 2.0000000000000036, 1 1))', 0);
SELECT @l5.STIntersects(@l6) as IntersectionResult;
重疊
- STWithin( ) 與 Intersect 的意義相近，但是定義上更為嚴謹，兩個比對的個體須完全重疊，當符合條件時，回傳 1，反之 0。
DECLARE @g geometry;
DECLARE @h geometry;
SET @g = geometry::STGeomFromText('POLYGON((0 0, 2 0, 2 2, 0 2, 0 0))', 0);
SET @h = geometry::STGeomFromText('POLYGON((0 0, 2 0, 2 1.9, 0 2, 0 0))', 0);
SELECT @g.STWithin(@h);
環域分析
- STContains( ) 用於判斷資料是否位於任意多邊形範圍內，專案上的應用包含自行繪製任意多邊形、矩形及圓形等等。與直接使用數學運算的方式相比，透過此方法計算出來的成果更加精準。
DECLARE @g geometry;
DECLARE @h geometry;
SET @g = geometry::STGeomFromText('POLYGON((0 0, 2 0, 2 2, 0 2, 0 0))', 0);
SET @h = geometry::STGeomFromText('LINESTRING(1 1,0 0)', 0);
SELECT @g.STContains(@h);