Category Archives: SQL Spatial

沒有長進

二十幾年前,初次學寫程式時,就嘗試畫多邊形,二十多年後,沒想到還在學用 T-SQL 畫多邊形…

–Geometry 沒有方向性
SELECT *,Geom.STArea(),Geom.ToString() FROM (
SELECT ‘內圈(逆時針方向)’  AS Label,
geometry::STPolyFromText(
‘POLYGON ( (10 10, 30 10, 30 30, 10 30, 10 10))’,4326) AS Geom
UNION ALL
SELECT REPLICATE(CHAR(13)+CHAR(10),12)+ ‘方圈(內圈順時針,外圈逆時針)’  AS Label,
geometry::STPolyFromText(
‘POLYGON ((0 0, 40 0, 40 40, 0 40, 0 0),(10 10, 10 30, 30 30, 30 10, 10 10))’,4326)
UNION ALL
SELECT REPLICATE(CHAR(13)+CHAR(10),20)+’外圈’  AS Label,
geometry::STPolyFromText(
‘POLYGON ( (-10 -10, 50 -10,50 50,-10 50, -10 -10),(0 0, 0 40, 40 40, 40 0, 0 0) )’,4326).MakeValid()  — Geometry 外圍的要在外圈
) T

image

 

Geography 則 Ring 有方向性

–Geography 有方向性
SELECT *,Geom.STArea(),Geom.ToString() FROM (
SELECT ‘內圈(逆時針方向)’  AS Label,
geography::STPolyFromText(
‘POLYGON ( (10 10, 30 10, 30 30, 10 30, 10 10))’,4326) AS Geom
UNION ALL
SELECT REPLICATE(CHAR(13)+CHAR(10),12)+ ‘方圈(內圈順時針,外圈逆時針)’  AS Label,
geography::STPolyFromText(
‘POLYGON ((10 10, 10 30, 30 30, 30 10, 10 10),(0 0, 40 0, 40 40, 0 40, 0 0))’,4326)
UNION ALL
SELECT REPLICATE(CHAR(13)+CHAR(10),20)+’外圈’  AS Label,
geography::STPolyFromText(
‘POLYGON ( (0 0, 0 40, 40 40, 40 0, 0 0),(-10 -10, 50 -10,50 50,-10 50, -10 -10) )’,4326)
) T

image

透過 .NET 程式碼存取空間資料

在伺服器端定義的資料表如下:

create database SpatialLab
go

use SpatialLab
go

create table Geo(GeoKey uniqueidentifier not null, GeoData geometry)
go

create procedure up_Geo_Insert @GeoKey uniqueidentifier, @GeoData Geometry
as
insert into Geo (GeoKey,GeoData) values (@GeoKey,@GeoData)
go

.NET 程式碼

設定參照 Microsoft.SqlServer.Types.dll ,預設位置在 c:Program Files Microsoft SQL ServerMSSQL10.MSSQLServerMSSQLBinn.

專案加入
using Microsoft.SqlServer.Types;

在程式內便可以使用 SqlGeometry 等類別

SqlConnection con = new SqlConnection("Data source=.;Trusted_Connection=True;Initial Catalog=SpatialLab");
con.Open();
using (SqlCommand cmd = new SqlCommand("truncate table SpatialLab.dbo.geo", con)){
cmd.ExecuteNonQuery();}

using (SqlCommand cmd = new SqlCommand("dbo.up_Geo_Insert", con)){

cmd.CommandType = CommandType.StoredProcedure;
cmd.Parameters.Add(new SqlParameter("@GeoKey", SqlDbType.UniqueIdentifier));

//必須要指定型態為 Udt
cmd.Parameters.Add(new SqlParameter("@GeoData", SqlDbType.Udt));
//再透過 UdtTypeName 設定自訂資料型態的名稱
cmd.Parameters[1].UdtTypeName = "Geometry";

//Geoms 宣告為 BindingList<SqlGeometry> 是存放 SqlGeometry 物件執行個體的 List

foreach (SqlGeometry Geom in Geoms){

cmd.Parameters[0].Value = Guid.NewGuid();
cmd.Parameters[1].Value = Geom;
cmd.ExecuteNonQuery();
}

}
con.Close();

 

從 DB 讀回:

private void DrawFromDB_Click(object sender, RoutedEventArgs e)
{
    SqlConnection con = new SqlConnection("Data source=.;Trusted_Connection=True;Initial Catalog=SpatialLab");
    con.Open();
    using (SqlCommand cmd = new SqlCommand("SELECT GeoData FROM Geo", con))
    {
        SqlDataReader drGeoms = cmd.ExecuteReader();
        while (drGeoms.Read())
        {
            Geoms.Add((SqlGeometry)drGeoms.GetValue(0));

 

若要將地址轉成經緯度,也可以在 Google Map API 註冊取得 Key 後,參考 blog: 利用 Google Maps 查詢地址經緯度 – Geocoding via HTTP 簡易範例,即可利用 .NET 程式將 DB 內的地址加上經緯度微笑

我所蒐集到存取空間資料的範例程式碼

空間索引(Spatial Index)

空間資料型態是用 .NET 所撰寫,提供豐富功能的類別,則空間索引則是該資料型態進資料庫的重要考量。以下節錄 SQL Server 2008 對空間索引的說明:ms-help://MS.SQLCC.v10/MS.SQLSVR.v10.en/s10de_1devconc/html/16a32f47-961e-4cba-b60c-11303acee1ab.htm (強迫自己好好讀一遍微笑 )

格階層(Grid Hierarchy)

在 SQL Server 2008 空間索引依然採用 B-Tree 結構,這代表 2 維的的空間資料要以線性的 B-Tree 來表現。因此 SQL Server 將空間資料以階層(level)的方式解構。例如先以 Grid 結構均分 16 格,每格再切分 16 格

image

如圖 Grid 階層經四層切割後,是 16 的四次方,總 65536 格。

Grid 階層中,所有的格子(cell)透過型變(variation of)過的 Hibert space-filling curve 來線性編號。但為了簡化說明,此處改以列(row-wise)為主的順序編號,從左上方開始編。在圖形中,以多邊形代表建築物,而線條代表道路,此為 4*4 階層 1 的圖示:

image

格密度(Grid Density)

軸上所劃分的格子數定義了格密度,數量越大,格密度越大,而 T-SQL 的 CREATE SPATIAL INDEX 語法的 GRIDS 子句可以讓你一各層定義不同的格密度。其關鍵字意義如下表:

關鍵字

格設定

格子(cell)的數量

Low

4*4

16

Medium

8*8

64

High

16*16

256

若未透過 GRIDS 子句指定格密度,預設為 Medium

語法範例如下:

CREATE SPATIAL INDEX SIndx_SpatialTable_geography_col3
   ON SpatialTable2(object)
   WITH ( GRIDS = ( LEVEL_3 = HIGH, LEVEL_2 = HIGH ) );

棋盤(tessellation)

當切割要建索引的空間成格階層後,就執行棋盤程序(tessellation process),將空間資料套入格階層,從第一層開始執行 breadth first,找空間物件有碰到哪些格子(touched cells),逐層往下做,一次一層。棋盤程序將物件所碰觸的格子記錄進空間索引,藉由這些記錄下來的格子,透過空間索引可以比較兩個空間物件的關係。

棋盤規則(Tessellation Rules)

為了局限紀錄空間物件的被碰觸格子(touched cell)之數量,棋盤程序使用了如下的棋盤規則,以決定棋盤程序的深度(depth),和索引要記錄哪些被碰觸格子:

  • 覆蓋(covering)規則:如果空間物件完整覆蓋某個格子,則該格子會被記錄下來,且不再進行進一步的棋盤程序。如圖中第一層的第 15 格,在第二層的完整覆蓋第 11 格。覆蓋規則將簡化棋盤程序,並減少空間索引紀錄的資料量。

image

  • Cell-per-object 規則:強制每個物件可使用的格子最大數量。除了階層 1 外,Cell-per-object 規則將會控制每個物件可以記錄的資訊量。在建立索引時,可以透過 CELLS_PER_OBJECT = n 子句來定義,預設是 16,以取得容量和精確度的平衡,其值域為 1 – 8192。T-SQL 語法範例如下:

    CREATE SPATIAL INDEX SIndx_SpatialTable_geography_col2
       ON SpatialTable2(object)
       USING GEOGRAPHY_GRID
       WITH (
        GRIDS = (MEDIUM, LOW, MEDIUM, HIGH ),
        CELLS_PER_OBJECT = 64,
        PAD_INDEX  = ON );

    由於第 1 層可以不接受限制,所以當第 1 層的數量達到或超過限制時,棋盤程序就會停止,而不再做下一階層的掃描。

    由於格子的數量要小於 CELLS_PER_OBJECT 的設定,所以當執行某個格子的棋盤程序後,其記錄的子格子數量將會超過限制,則該程序就會放棄,僅記錄父格子。以上圖為例,如 CELLS_PER_OBJECT 設定為 9,則於第一層紀錄 15 後,針對格子 15 進行棋盤程序,當要進入再下一層時,需要 9 個格子來記錄,超過總數量,因此放棄下一層,而在索引中僅記錄格子 15。

     

  • deepest-cell 規則:每一個低階的格子隸屬於高階的格子,例如 4.4.10.13 隸屬 4.4.10,4.4.10 隸屬 4.4,4.4 隸屬 4。由於 Query Processor 可以理解格子的階層關係,因此索引中只需要記載最下層的格子,以減少須存放的資訊。

image

以上圖為例,索引的 CELLS_PER_OBJECT 選項採用了預設的 16,而記錄物件所花的格子並未超過這個數量,因此可進行棋盤程序到最下一層,套用 deepest-cell 規則後,僅需要記載 12 個階層四的格子:4.4.10.13-15 、 4.4.14.1-3、 4.4.14.5-7、 4.4.14.9-11

Bounding Box

幾何資料(geometric)可以定義無限大空間,但 SQL Server 2008 的空間索引需要有限空間(finite space),為了解構(decomposition)而定義一個有限空間,geometry grid tessellation scheme 需要一個長方形的 bounding box;以座標的右上和左下點 (x-min,y-min) 和 (x-max,y-max) 來定義。

 image

語法範例如下:
CREATE SPATIAL INDEX SIndx_SpatialTable_geometry_col1
   ON SpatialTable(geometry_col)
   WITH ( BOUNDING_BOX = ( 0, 0, 500, 200 ) );

建立了空間索引後,以下的 geometry 方法可以較有效率地處理集合導向資料:

  • geometry1.STContains( geometry2 ) = 1
  • geometry1.STDistance( geometry2 ) < number
  • geometry1.STDistance( geometry2 ) <= number
  • geometry1.STEquals( geometry2 ) = 1
  • geometry1.STIntersects( geometry2 ) = 1
  • geometry1.STOverlaps(geometry2) = 1
  • geometry1.STTouches( geometry2 ) = 1
  • geometry1.STWithin( geometry2 ) = 1

使用方式必須是在 WHERE 子句後,以下述格式設定條件:
geometry1.方法( geometry2 )  比較運算子  數值

geography 則是 STIntersects(), STEquals(), and STDistance() 等方法可以使用空間索引

簡單的使用範例

搭配以下網頁 http://mikeo.co.uk/demo/sqlspatial/default.aspx 取得臺灣各地理座標(取多邊形的點時,要逆時針取)

image

在此建立存放地理資訊的資料表,並建立空間索引,然後輸入台北市各地標的經緯度資料,再比對各物件是否有交錯(STIntersects)

USE tempdb
go
CREATE TABLE tbGeography(
PK INT IDENTITY(1,1) PRIMARY KEY,
Descriptions NVARCHAR(100),
Geographys GEOGRAPHY);
go
CREATE SPATIAL INDEX SIndx ON tbGeography(Geographys);
–DROP INDEX SIndx ON tbGeography;

SELECT * FROM sys.spatial_reference_systems
WHERE well_known_text LIKE ‘%Taiwan%’

INSERT tbGeography(Descriptions, Geographys)VALUES (N’台北 101 大樓’,
geography::STGeomFromText(‘POINT(25.034127684732052 121.56410694122314)’, 4326))

INSERT into tbGeography (Descriptions, Geographys) VALUES (N’基隆路’,
geography::STGeomFromText(‘LINESTRING(25.042681800625725 121.56569480895997,
25.042759562579235 121.56586647033695, 25.028217219856525 121.55702590942383,
25.009628508097456 121.53462409973144)’, 4326))

INSERT into tbGeography (Descriptions, Geographys) VALUES (N’大安森林公園’,
geography::STGeomFromText(‘POLYGON((25.033369450224278 121.53292894363406,
25.026214596461695 121.53513908386233, 25.025884063244827 121.53743505477908,
25.033233356355037 121.53756380081182, 25.033369450224278 121.53292894363406))’, 4326))

DECLARE @g GEOGRAPHY = geography::STGeomFromText –建國南路
(‘LINESTRING(25.037335548101733 121.5373492240906, 25.036693982138722 121.5373492240906,
25.026972875183723 121.537070274353)’, 4326)

SELECT * FROM tbGeography WHERE Geographys.STIntersects(@g)= 1

結果如下:

image

TRUNCATE TABLE tbGeography

SQL Server 2008 的空間資料

以 .NET 實做的空間(Spatial)資料格式:提供兩種資料型態。可在資料庫內結合地理資訊,並搭配空間索引(Spatial Index)的特殊階層式索引結構,以有效存取資料。

  • Geometry 用在 2 維描述平面地球,符合如 Open Geospatial Consortium(OGC)對空間資料的定義標準,存放多邊形(polygon)、線(line)和點(point)的資訊。
  • Geography 則以 3 維立體的描述圓地球,存放緯度(latitude)和經度(longitude)的座標資料,符合產業標準 WGS84,提供 GPS 對地球表面定位的資訊。

存放的幾何形狀有:POINT, MULTIPOINT, LINESTRING, MULTILINESTRING, POLYGON, MULTIPOLYGON, GEOMETRYCOLLECTION

image

雖結構有 11 ,但可成為執行個體的只有 7 種,雖然圖中的根為 Geometry,但 Geography 亦是如此

而表現資料的格式有三種:

  • Well Know Text(WKT):易讀,但不容易用程式產生,因為要注意格式。
    • 以文字描述幾何圖形,不分大小寫
    • 點是以空白相隔的兩個數字,如:(X Y)
      • 若是 Geography,則 X 代表經度,Y 代表緯度
    • 可以在 X Y 之後添加 Z 和 M 屬性,Z 和 M 可以不填,但 X Y 不行
    • 多個點間以小括號括起,並以逗號相隔
    • 若是 POLYGON, MULTIPOLYGON, GEOMETRYCOLLECTION,再以小括號括起多個點所成集合
    • 一組資料可以有多個多邊形,但邊界不能重疊
    • 範例
      • declare @point geometry = ‘POINT(100 10)’
      • declare @line geometry = ‘LINESTRING(0 0,10 10)’
      • declare @square geometry = ‘POLYGON((0 0,10 0,10 10,0 10,0 0))’
      • declare @donught geometry = ‘POLYGON((0 0,10 0,10 10,0 10,0 0),(4 4,4 6,6 6,6 4,4 4))’
      • declare @vs geometry = ‘MULTILINESTRING((0 10,5 0,10 10),(20 10,25 0,30 10))’
      • declare @rockets geometry = ‘MULTIPOLYGON(((0 0,5 10,5 30,15 40,25 30,25 10,30 0,20 0,15 5, 5 0,0 0)),((10 0,15 10,15 30,25 40,35 30,35 10,40 0,40 0,25 5, 15 0,10 0)))’
      • declare @glasses geometry = ‘GEOMETRYCOLLECTION(LINESTRING(0 20,5 10),POLYGON((5 5,5 15,10 20,20 20,25 15,25 5,20 0,10 0,5 5)),POLYGON((35 5,35 15,40 20,50 20,55 15,55 5,50 0,40 0,35 5)),LINESTRING(55 10,60 20),LINESTRING(25 15,35 15))’
  • Well Known Binary(WKB):以 2 進位的位元組表現,透過編碼解釋結構,加上數值資料,各種幾何結構定義如下

    Point = [Byte order][Type][X value][Y value]

    MultiPoint = [Byte order][Type][PointCount]<Point1><Point2>

    Line = [Byte order][Type][PointCount][X Value][YValue][X value][Y value]…..

    MultiLine = [Byte order][Type][LineCount]<Line1><Line2> …..

    Polygon = [Byte order][Type][Polygon Count][PointCount][X Value][YValue][X value][Y value]…..

    MultiPolygon = [Byte order][Type][Polygon Count]<Polygon1><Polygon2>….

    Geometry Collection = [Byte order][Type][Shape Count]<Shape1><Shape2>…..

    • 而 Type 的定義如下:

      • Unknown = 0
      • Point = 1
      • LineString = 2
      • Polygon = 3
      • MultiPoint = 4
      • MultiLineString = 5
      • MultiPolygon = 6
      • GeometryCollection = 7

      範例如下:

      declare @p geometry = geometry::STGeomFromWKB

            AABBBBBBBBCCCCCCCCDDDDDDDD
      (0x00000000030000000200000005

      X1X1X1X1X1X1X1X1Y1Y1Y1Y1Y1Y1Y1Y1X2X2X2X2X2X2X2X2Y2Y2Y2Y2Y2Y2Y2Y2
      0000000000000000000000000000000000000000000000004010000000000000

      X3X3X3X3X3X3X3X3Y3Y3Y3Y3Y3Y3Y3Y3X4X4X4X4X4X4X4X4Y4Y4Y4Y4Y4Y4Y4Y4
      4010000000000000401000000000000040100000000000000000000000000000

      X5X5X5X5X5X5X5X5Y5Y5Y5Y5Y5Y5Y5Y5
      00000000000000000000000000000000

      DDDDDDDD
      00000005

      X1….

      3FF00000000000003FF00000000000003FF00000000000004008000000000000

      4008000000000000400800000000000040080000000000003FF0000000000000

      3FF00000000000003FF0000000000000,0)

      print @p.ToString()

      上述資料經 SQL Server 2008 @p.ToString() 印出的 WKT 表示法為:
      POLYGON ((0 0, 0 4, 4 4, 4 0, 0 0), (1 1, 1 3, 3 3, 3 1, 1 1))

      範例中上方的 ABC 各代表意義如下:

      • A:位元組順序(byte order)
      • B:幾何形狀的格式(type),此處 3 代表多邊形(polygon)
      • C:子形狀(sub shape)的數量,此處 2 代表有兩個多邊形
      • D:子形狀的點數數量,此處 5 代表由五個點組成多邊形
      • X1 Y1…Xn Yn 分別代表 X Y 所標示的點
      • 5 個點過後,再度標示下一個多邊形的點數,然後就是各點的資料

      declare @p geometry = geometry::STGeomFromWKB
      (0x000000000300000002000000050000000000000000000000000000000000000000000000004010000000000000401000000000000040100000000000004010000000000000000000000000000000000000000000000000000000000000000000053FF00000000000003FF00000000000003FF000000000000040080000000000004008000000000000400800000000000040080000000000003FF00000000000003FF00000000000003FF0000000000000,0)

      Geography Markup Language(GML) 以 XML 格式描述:.

      declare @p geometry;

      set @p = geometry::GeomFromGml(‘<Point xmlns="http://www.opengis.net/gml"><pos>10 10</pos></Point>’,0);

       

      透過 T-SQL 組 WKB 的內容

      declare @i int = cast(cast(1 as binary(8)) as int)

      declare @OuterX1 float = 0, @OuterY1 float = 0
      declare @OuterX2 float = 400, @OuterY2 float = 0
      declare @OuterX3 float = 400, @OuterY3 float = 400
      declare @OuterX4 float = 0, @OuterY4 float = 400
      declare @OuterX5 float = 0, @OuterY5 float = 0

      declare @InnerX1 float = 100, @InnerY1 float = 100
      declare @InnerX2 float = 300, @InnerY2 float = 100
      declare @InnerX3 float = 300, @InnerY3 float = 300
      declare @InnerX4 float = 100, @InnerY4 float = 300
      declare @InnerX5 float = 100, @InnerY5 float = 100

      declare @p geometry


      print 0x000000000300000002
        + 0x00000005
        + cast(@OuterX1 as binary(8)) + cast(@OuterY1 as binary(8))
        + cast(@OuterX2 as binary(8)) + cast(@OuterY2 as binary(8))
        + cast(@OuterX3 as binary(8)) + cast(@OuterY3 as binary(8))
        + cast(@OuterX4 as binary(8)) + cast(@OuterY4 as binary(8))
        + cast(@OuterX5 as binary(8)) + cast(@OuterY5 as binary(8))
        + 0x00000005
        + cast(@InnerX1 as binary(8)) + cast(@InnerY1 as binary(8))
        + cast(@InnerX2 as binary(8)) + cast(@InnerY2 as binary(8))
        + cast(@InnerX3 as binary(8)) + cast(@InnerY3 as binary(8))
        + cast(@InnerX4 as binary(8)) + cast(@InnerY4 as binary(8))
        + cast(@InnerX5 as binary(8)) + cast(@InnerY5 as binary(8))
      /*
      0x00000000030000000200000005000000000000000000000000000000004079000000000000000000000000000040790000000000004079000000000000000000000000000040790000000000000000000000000000000000000000000000000005405900000000000040590000000000004072C0000000000040590000000000004072C000000000004072C0000000000040590000000000004072C0000000000040590000000000004059000000000000
      */
      set @p = geometry::STGeomFromWKB( 
          0x000000000300000002
        + 0x00000005
        + cast(@OuterX1 as binary(8)) + cast(@OuterY1 as binary(8))
        + cast(@OuterX2 as binary(8)) + cast(@OuterY2 as binary(8))
        + cast(@OuterX3 as binary(8)) + cast(@OuterY3 as binary(8))
        + cast(@OuterX4 as binary(8)) + cast(@OuterY4 as binary(8))
        + cast(@OuterX5 as binary(8)) + cast(@OuterY5 as binary(8))
        + 0x00000005
        + cast(@InnerX1 as binary(8)) + cast(@InnerY1 as binary(8))
        + cast(@InnerX2 as binary(8)) + cast(@InnerY2 as binary(8))
        + cast(@InnerX3 as binary(8)) + cast(@InnerY3 as binary(8))
        + cast(@InnerX4 as binary(8)) + cast(@InnerY4 as binary(8))
        + cast(@InnerX5 as binary(8)) + cast(@InnerY5 as binary(8)),0)

      print @p.ToString()
      /*
      POLYGON ((0 0, 400 0, 400 400, 0 400, 0 0), (100 100, 300 100, 300 300, 100 300, 100 100))
      */

      透過 Simon Sabin 提供的 spatialViewer 檢視該多邊形如下

      image

      Geography 資料型態是以經緯度來表示,但相關函數需要換算長、寬距離…等。例如從 http://mikeo.co.uk/demo/sqlspatial/default.aspx 取得台北與高雄兩點的緯度與經度:

      image

      直接透過 Geography 資料型態的函數換算距離:

      DECLARE @Taipei GEOGRAPHY=’POINT(25.015928763367846 121.58569335937504)’  — 宣告台北的緯度與經度
      DECLARE @Kaohsiung GEOGRAPHY=’POINT(22.614010874370265 120.33325195312501)’
      SELECT @Kaohsiung.STDistance(@Taipei) — 取回在地球上兩點間的距離
      回傳值:295047.011615369(單位:公尺)

      而非簡單的兩個座標點:

      DECLARE @x FLOAT=(25.015928763367846-22.614010874370265)
      DECLARE @y FLOAT=(121.58569335937504-120.33325195312501)
      SELECT sqrt(square(@x)+square(@y))
      回傳值:2.70884090001169

      若僅是簡單的兩個座標點,也就是一般的平面幾何

      DECLARE @Taipei GEOMETRY=’POINT(25.015928763367846 121.58569335937504)’
      DECLARE @Kaohsiung GEOMETRY=’POINT(22.614010874370265 120.33325195312501)’
      SELECT @Kaohsiung.STDistance(@Taipei)
      其回傳值就是上述的 2.70884090001169

      另外,在標示 Geometry 的座標點時,雖然可以標示 Z 座標和度量單位(M measure),但實際在運算時,並沒有作用,僅讓使用者可以透過屬性參考到

      DECLARE @T GEOMETRY=’POINT(0 0 0 10)’
      SELECT @T.STX
      SELECT @T.STY
      SELECT @T.Z  — 取出 Z 的資料
      SELECT @T.M  –取出 m 的資料

      DECLARE @K GEOMETRY=’POINT(10 10 10 10)’
      SELECT @K.STDistance(@T)
      結果依然是 200 開根號:14.142135623731