一、PostGIS介紹
PostGIS是對象關系型數據庫系統PostgreSQL的一個擴展,PostGIS提供如下空間信息服務功能:空間對象、空間索引、空間操作函數和空間操作符。
同時,PostGIS遵循OpenGIS的規范。
二、 PostGIS中的幾何類型
PostGIS支持所有OGC(Open Geospatial Consortium) 規范的“Simple Features”類型,同時在此基礎上擴展了對3DZ、3DM、4D坐標的支持。
1. OGC的WKB和WKT格式
OGC定義了兩種描述幾何對象的格式,分別是WKB(Well-Known Binary)和WKT(Well-Known Text)。
在SQL語句中,用以下的方式可以使用WKT格式定義幾何對象:
POINT(0 0) //點 LINESTRING(0 0,1 1,1 2) //線 POLYGON((0 0,4 0,4 4,0 4,0 0),(1 1, 2 1, 2 2, 1 2,1 1)) //面 MULTIPOINT(0 0,1 2) //多點 MULTILINESTRING((0 0,1 1,1 2),(2 3,3 2,5 4)) //多線 MULTIPOLYGON(((0 0,4 0,4 4,0 4,0 0),(1 1,2 1,2 2,1 2,1 1)), ((-1 -1,-1 -2,-2 -2,-2 -1,-1 -1)))//多面 GEOMETRYCOLLECTION(POINT(2 3),LINESTRING((2 3,3 4))) //幾何集合
以下語句可以使用WKT格式插入一個點要素到一個表中,其中用到的GeomFromText等函數在后面會有詳細介紹:
INSERT INTO table ( SHAPE, NAME ) VALUES ( GeomFromText('POINT(116.39 39.9)', 4326), '北京');
2. EWKT、EWKB和Canonical格式
EWKT和EWKB相比OGC,WKT和WKB格式主要的擴展有3DZ、3DM、4D坐標和內嵌空間參考支持。
以下以EWKT語句定義了一些幾何對象:
POINT(0 0 0) //3D點 SRID=32632;POINT(0 0) //內嵌空間參考的點 POINTM(0 0 0) //帶M值的點 POINT(0 0 0 0) //帶M值的3D點 SRID=4326;MULTIPOINTM(0 0 0,1 2 1) //內嵌空間參考的帶M值的多點
以下語句可以使用EWKT格式插入一個點要素到一個表中:
INSERT INTO table ( SHAPE, NAME ) VALUES ( GeomFromEWKT('SRID=4326;POINTM(116.39 39.9 10)'), '北京' )
Canonical格式是16進制編碼的幾何對象,直接用SQL語句查詢出來的就是這種格式。
3. SQL-MM格式
SQL-MM格式定義了一些插值曲線,這些插值曲線和EWKT有點類似,也支持3DZ、3DM、4D坐標,但是不支持嵌入空間參考。
以下以SQL-MM語句定義了一些插值幾何對象:
CIRCULARSTRING(0 0, 1 1, 1 0) //插值圓弧 COMPOUNDCURVE(CIRCULARSTRING(0 0, 1 1, 1 0),(1 0, 0 1)) //插值復合曲線 CURVEPOLYGON(CIRCULARSTRING(0 0, 4 0, 4 4, 0 4, 0 0),(1 1, 3 3, 3 1, 1 1)) //曲線多邊形 MULTICURVE((0 0, 5 5),CIRCULARSTRING(4 0, 4 4, 8 4)) //多曲線 MULTISURFACE(CURVEPOLYGON(CIRCULARSTRING(0 0, 4 0, 4 4, 0 4, 0 0),(1 1, 3 3, 3 1, 1 1)),((10 10, 14 12, 11 10, 10 10),(11 11, 11.5 11, 11 11.5, 11 11))) //多曲面
三、 PostGIS中空間信息處理的實現
1. spatial_ref_sys表
在基於PostGIS模板創建的數據庫的public模式下,有一個spatial_ref_sys表,它存放的是OGC規范的空間參考。我們取我們最熟悉的4326參考看一下:
它的srid存放的就是空間參考的Well-Known ID,對這個空間參考的定義主要包括兩個字段,srtext存放的是以字符串描述的空間參考,proj4text存放的則是以字符串描述的PROJ.4 投影定義(PostGIS使用PROJ.4實現投影)。
4326空間參考的srtext內容:
GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]
4326空間參考的proj4text內容:
+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
2. geometry_columns表
geometry_columns表存放了當前數據庫中所有幾何字段的信息,比如我當前的庫里面有兩個空間表,在geometry_columns表中就可以找到這兩個空間表中幾何字段的定義:
其中f_table_schema字段表示的是空間表所在的模式,f_table_name字段表示的是空間表的表名,f_geometry_column字段表示的是該空間表中幾何字段的名稱,srid字段表示的是該空間表的空間參考。
3. 在PostGIS中創建一個空間表
在PostGIS中創建一個包含幾何字段的空間表分為2步:第一步創建一個一般表,第二步給這個表添加幾何字段。
以下先在test模式下創建一個名為cities的一般表:
create table test.cities (id int4, name varchar(20))
再給cities添加一個名為shape的幾何字段(二維點):
select AddGeometryColumn('test', 'cities', 'shape', 4326, 'POINT', 2)
4. PostGIS對幾何信息的檢查
PostGIS可以檢查幾何信息的正確性,這主要是通過IsValid函數實現的。
以下語句分辨檢查了2個幾何對象的正確性,顯然,(0, 0)點和(1,1)點可以構成一條線,但是(0, 0)點和(0, 0)點則不能構成,這個語句執行以后的得出的結果是TRUE,FALSE。
select IsValid('LINESTRING(0 0, 1 1)'), IsValid('LINESTRING(0 0,0 0)')
默認PostGIS並不會使用IsValid函數檢查用戶插入的新數據,因為這會消耗較多的CPU資源(特別是復雜的幾何對象)。當你需要使用這個功能的時候,你可以使用以下語句為表新建一個約束:
ALTER TABLE cities ADD CONSTRAINT geometry_valid CHECK (IsValid(shape))
這時當我們往這個表試圖插入一個錯誤的空間對象的時候,會得到一個錯誤:
INSERT INTO test.cities ( shape, name ) VALUES ( GeomFromText('LINESTRING(0 0,0 0)', 4326), '北京'); ERROR: new row for relation "cities" violates check constraint "geometry_valid" SQL 狀態: 23514
5. PostGIS中的空間索引
數據庫對多維數據的存取有兩種索引方案,R-Tree和GiST(Generalized Search Tree),在PostgreSQL中的GiST比R-Tree的健壯性更好,因此PostGIS對空間數據的索引一般采用GiST實現。
以下的語句給sde模式中的cities表添加了一個空間索引shape_index_cities,在pgAdmin中也可以通過圖形界面完成相同的功能。
CREATE INDEX shape_index_cities ON sde.cities USING gist (shape);
另外要注意的是,空間索引只有在進行基於邊界范圍的查詢時才起作用,比如“&&”操作。
四、 PostGIS中的常用函數
首先需要說明一下,這里許多函數是以ST_[X]yyy形式命名的,事實上很多函數也可以通過xyyy的形式訪問,在PostGIS的函數庫中我們可以看到這兩種函數定義完全一樣。
1. OGC標准函數
管理函數:
添加幾何字段 AddGeometryColumn(, , , , , ) 刪除幾何字段 DropGeometryColumn(, , ) 檢查數據庫幾何字段並在geometry_columns中歸檔 Probe_Geometry_Columns() 給幾何對象設置空間參考(在通過一個范圍做空間查詢時常用) ST_SetSRID(geometry, integer)
幾何對象關系函數 :
獲取兩個幾何對象間的距離 ST_Distance(geometry, geometry) 如果兩個幾何對象間距離在給定值范圍內,則返回TRUE ST_DWithin(geometry, geometry, float) 判斷兩個幾何對象是否相等(比如LINESTRING(0 0, 2 2)和LINESTRING(0 0, 1 1, 2 2)是相同的幾何對象) ST_Equals(geometry, geometry) 判斷兩個幾何對象是否分離 ST_Disjoint(geometry, geometry) 判斷兩個幾何對象是否相交 ST_Intersects(geometry, geometry) 判斷兩個幾何對象的邊緣是否接觸 ST_Touches(geometry, geometry) 判斷兩個幾何對象是否互相穿過 ST_Crosses(geometry, geometry) 判斷A是否被B包含 ST_Within(geometry A, geometry B) 判斷兩個幾何對象是否是重疊 ST_Overlaps(geometry, geometry) 判斷A是否包含B ST_Contains(geometry A, geometry B) 判斷A是否覆蓋 B ST_Covers(geometry A, geometry B) 判斷A是否被B所覆蓋 ST_CoveredBy(geometry A, geometry B) 通過DE-9IM 矩陣判斷兩個幾何對象的關系是否成立 ST_Relate(geometry, geometry, intersectionPatternMatrix) 獲得兩個幾何對象的關系(DE-9IM矩陣) ST_Relate(geometry, geometry)
幾何對象處理函數:
獲取幾何對象的中心 ST_Centroid(geometry) 面積量測 ST_Area(geometry) 長度量測 ST_Length(geometry) 返回曲面上的一個點 ST_PointOnSurface(geometry) 獲取邊界 ST_Boundary(geometry) 獲取緩沖后的幾何對象 ST_Buffer(geometry, double, [integer]) 獲取多幾何對象的外接對象 ST_ConvexHull(geometry) 獲取兩個幾何對象相交的部分 ST_Intersection(geometry, geometry) 將經度小於0的值加360使所有經度值在0-360間 ST_Shift_Longitude(geometry) 獲取兩個幾何對象不相交的部分(A、B可互換) ST_SymDifference(geometry A, geometry B) 從A去除和B相交的部分后返回 ST_Difference(geometry A, geometry B) 返回兩個幾何對象的合並結果 ST_Union(geometry, geometry) 返回一系列幾何對象的合並結果 ST_Union(geometry set) 用較少的內存和較長的時間完成合並操作,結果和ST_Union相同 ST_MemUnion(geometry set)
幾何對象存取函數:
獲取幾何對象的WKT描述 ST_AsText(geometry) 獲取幾何對象的WKB描述 ST_AsBinary(geometry) 獲取幾何對象的空間參考ID ST_SRID(geometry) 獲取幾何對象的維數 ST_Dimension(geometry) 獲取幾何對象的邊界范圍 ST_Envelope(geometry) 判斷幾何對象是否為空 ST_IsEmpty(geometry) 判斷幾何對象是否不包含特殊點(比如自相交) ST_IsSimple(geometry) 判斷幾何對象是否閉合 ST_IsClosed(geometry) 判斷曲線是否閉合並且不包含特殊點 ST_IsRing(geometry) 獲取多幾何對象中的對象個數 ST_NumGeometries(geometry) 獲取多幾何對象中第N個對象 ST_GeometryN(geometry,int) 獲取幾何對象中的點個數 ST_NumPoints(geometry) 獲取幾何對象的第N個點 ST_PointN(geometry,integer) 獲取多邊形的外邊緣 ST_ExteriorRing(geometry) 獲取多邊形內邊界個數 ST_NumInteriorRings(geometry) 同上 ST_NumInteriorRing(geometry) 獲取多邊形的第N個內邊界 ST_InteriorRingN(geometry,integer) 獲取線的終點 ST_EndPoint(geometry) 獲取線的起始點 ST_StartPoint(geometry) 獲取幾何對象的類型 GeometryType(geometry) 類似上,但是不檢查M值,即POINTM對象會被判斷為point ST_GeometryType(geometry) 獲取點的X坐標 ST_X(geometry) 獲取點的Y坐標 ST_Y(geometry) 獲取點的Z坐標 ST_Z(geometry) 獲取點的M值 ST_M(geometry)
幾何對象構造函數 :
ST_PointFromText(text,[]) ST_LineFromText(text,[]) ST_LinestringFromText(text,[]) ST_PolyFromText(text,[]) ST_PolygonFromText(text,[]) ST_MPointFromText(text,[]) ST_MLineFromText(text,[]) ST_MPolyFromText(text,[]) ST_GeomCollFromText(text,[]) ST_GeomFromWKB(bytea,[]) ST_GeometryFromWKB(bytea,[]) ST_PointFromWKB(bytea,[]) ST_LineFromWKB(bytea,[]) ST_LinestringFromWKB(bytea,[]) ST_PolyFromWKB(bytea,[]) ST_PolygonFromWKB(bytea,[]) ST_MPointFromWKB(bytea,[]) ST_MLineFromWKB(bytea,[]) ST_MPolyFromWKB(bytea,[]) ST_GeomCollFromWKB(bytea,[]) ST_BdPolyFromText(text WKT, integer SRID) ST_BdMPolyFromText(text WKT, integer SRID)
參考語義:
Text:WKT WKB:WKB Geom:Geometry M:Multi Bd:BuildArea Coll:Collection ST_GeomFromText(text,[])
2. PostGIS擴展函數
管理函數:
刪除一個空間表(包括geometry_columns中的記錄) DropGeometryTable([], ) 更新空間表的空間參考 UpdateGeometrySRID([], , , ) 更新空間表的統計信息 update_geometry_stats([, ])
postgis_version()
postgis_lib_version()
postgis_lib_build_date()
postgis_script_build_date()
postgis_scripts_installed()
postgis_scripts_released()
postgis_geos_version()
postgis_jts_version()
postgis_proj_version()
postgis_uses_stats()
postgis_full_version()
參考語義:
Geos:GEOS庫
Jts:JTS庫
Proj:PROJ4庫
幾何操作符:
A范圍=B范圍 A = B A范圍覆蓋B范圍或A范圍在B范圍左側 A &<> B A范圍在B范圍左側 A <<>> B A范圍覆蓋B范圍或A范圍在B范圍下方 A &<| B A范圍覆蓋B范圍或A范圍在B范圍上方 A |&> B A范圍在B范圍下方 A <<| B A范圍在B范圍上方 A |>> B A=B A ~= B A范圍被B范圍包含 A @ B A范圍包含B范圍 A ~ B A范圍覆蓋B范圍 A && B
幾何量測函數:
量測面積 ST_Area(geometry)
根據經緯度點計算在地球曲面上的距離,單位米,地球半徑取值6370986米 ST_distance_sphere(point, point)
類似上,使用指定的地球橢球參數 ST_distance_spheroid(point, point, spheroid)
量測2D對象長度 ST_length2d(geometry)
量測3D對象長度 ST_length3d(geometry)
根據經緯度對象計算在地球曲面上的長度 ST_length_spheroid(geometry,spheroid)
ST_length3d_spheroid(geometry,spheroid)
量測兩個對象間距離 ST_distance(geometry, geometry)
量測兩條線之間的最大距離 ST_max_distance(linestring,linestring)
量測2D對象的周長 ST_perimeter(geometry)
ST_perimeter2d(geometry)
量測3D對象的周長 ST_perimeter3d(geometry)
量測兩點構成的方位角,單位弧度 ST_azimuth(geometry, geometry)
幾何對象輸出:
ST_AsBinary(geometry,{'NDR'|'XDR'})
ST_AsEWKT(geometry)
ST_AsEWKB(geometry, {'NDR'|'XDR'})
ST_AsHEXEWKB(geometry, {'NDR'|'XDR'})
ST_AsSVG(geometry, [rel], [precision])
ST_AsGML([version], geometry, [precision])
ST_AsKML([version], geometry, [precision])
ST_AsGeoJson([version], geometry, [precision], [options])
參考語義:
NDR:Little Endian XDR:big-endian HEXEWKB:Canonical SVG:SVG 格式 GML:GML 格式 KML:KML 格式 GeoJson:GeoJson 格式
幾何對象創建:
ST_GeomFromEWKB(bytea) ST_MakePoint(, , [], []) ST_MakePointM(, , ) ST_MakeBox2D(, ) ST_MakeBox3D(, ) ST_MakeLine(geometry set) ST_MakeLine(geometry, geometry) ST_LineFromMultiPoint(multipoint) ST_MakePolygon(linestring, [linestring[]]) ST_BuildArea(geometry) ST_Polygonize(geometry set) ST_Collect(geometry set) ST_Collect(geometry, geometry) ST_Dump(geometry) ST_DumpRings(geometry)
參考語義:
Dump:轉儲 ST_GeomFromEWKT(text)
幾何對象編輯:
給幾何對象添加一個邊界,會使查詢速度加快 ST_AddBBOX(geometry) 刪除幾何對象的邊界 ST_DropBBOX(geometry) 添加、刪除、設置點 ST_AddPoint(linestring, point, []) ST_RemovePoint(linestring, offset) ST_SetPoint(linestring, N, point) 幾何對象類型轉換 ST_Force_collection(geometry) ST_Force_2d(geometry) ST_Force_3dz(geometry), ST_Force_3d(geometry), ST_Force_3dm(geometry) ST_Force_4d(geometry) ST_Multi(geometry) 將幾何對象轉化到指定空間參考 ST_Transform(geometry,integer) 對3D幾何對象作仿射變化 ST_Affine(geometry, float8, float8, float8, float8, float8, float8, float8, float8, float8, float8, float8, float8) 對2D幾何對象作仿射變化 ST_Affine(geometry, float8, float8, float8, float8, float8, float8) 對幾何對象作偏移 ST_Translate(geometry, float8, float8, float8) 對幾何對象作縮放 ST_Scale(geometry, float8, float8, float8) 對3D幾何對象作旋轉 ST_RotateZ(geometry, float8) ST_RotateX(geometry, float8) ST_RotateY(geometry, float8) 對2D對象作偏移和縮放 ST_TransScale(geometry, float8, float8, float8, float8) 反轉 ST_Reverse(geometry) 轉化到右手定則 ST_ForceRHR(geometry) 參考IsSimple函數 使用Douglas-Peuker算法 ST_Simplify(geometry, tolerance) ST_SimplifyPreserveTopology(geometry, tolerance) 講幾何對象頂點捕捉到網格 ST_SnapToGrid(geometry, originX, originY, sizeX, sizeY) ST_SnapToGrid(geometry, sizeX, sizeY), ST_SnapToGrid(geometry, size) 第二個參數為點,指定原點坐標 ST_SnapToGrid(geometry, geometry, sizeX, sizeY, sizeZ, sizeM) 分段 ST_Segmentize(geometry, maxlength) 合並為線 ST_LineMerge(geometry)
線性參考:
根據location(0-1)獲得該位置的點 ST_line_interpolate_point(linestring, location) 獲取一段線 ST_line_substring(linestring, start, end) 根據點獲取location(0-1) ST_line_locate_point(LineString, Point) 根據量測值獲得幾何對象 ST_locate_along_measure(geometry, float8) 根據量測值區間獲得幾何對象集合 ST_locate_between_measures(geometry, float8, float8)
雜項功能函數:
幾何對象的摘要 ST_Summary(geometry) 幾何對象的邊界 ST_box2d(geometry) ST_box3d(geometry) 多個幾何對象的邊界 ST_extent(geometry set) 0=2d, 1=3dm, 2=3dz, 3=4d ST_zmflag(geometry) 是否包含Bounding Box ST_HasBBOX(geometry) 幾何對象的維數:2、3、4 ST_ndims(geometry) 子對象的個數 ST_nrings(geometry) ST_npoints(geometry) 對象是否驗證成功 ST_isvalid(geometry) 擴大幾何對象 ST_expand(geometry, float) 計算一個空間表的邊界范圍 ST_estimated_extent([schema], table, geocolumn) 獲得空間參考 ST_find_srid(, , ) 幾何對象使用的內存大小,單位byte ST_mem_size(geometry) 點是否在圓上 ST_point_inside_circle(,,,) 獲取邊界的X、Y、Z ST_XMin(box3d) ST_YMin(box3d) ST_ZMin(box3d) ST_XMax(box3d) ST_YMax(box3d) ST_ZMax(box3d) 構造一個幾何對象的數組 ST_Accum(geometry set)
長事務支持:
啟用/關閉長事務支持,重復調用無副作用 EnableLongTransactions() DisableLongTransactions() 檢查對行的update和delete操作是否已授權 CheckAuth([], , ) 鎖定行 LockRow([], , , , []) 解鎖行 UnlockRows() 在當前事務中添加授權ID AddAuth()
其它還有SQL-MM和ArcSDE樣式的函數支持,可以參考http://postgis.refractions.net/documentation/manual-1.3/ch06.html#id2750611。
五、 PostGIS示例
下面我們通過一個簡單的Flex應用示例來看一下PostGIS的用法:
假想現在發生了恐怖襲擊,導致在一些城市有污染物出現,現在我們要根據污染物和當地風力、風向情況,計算污染擴散范圍,針對這些區域及時進行警報和疏散。
首先我們希望獲得所有發生污染的城市的當前風速、風向等信息,在我們的PostGIS數據庫中有一個空間表保存着這些信息,我們構造這樣的SQL語句進行查詢:
select *,ST_AsGeoJson(shape) from sde.wind
這里會獲取所有風相關的信息,並且附加了以JSON格式返回的幾何信息,這有助於我們在Flex中進行解析。如下圖是關於風的查詢結果:
下面我們希望PostGIS幫助我們實現一些空間分析。我們以污染發生的城市為起點,當地風向為主方向,構造一個30度開角的范圍;這個范圍將是污染擴散的主要方向,擴散的范圍主要和風的強度有關;在構造這個區域以后,為了保險起見,我們在對其進行一定范圍的緩沖,最后得到每個污染源可能擴散的范圍。我們構造的SQL語句如下:
select *,ST_AsGeoJson( ST_Buffer( ST_PolygonFromText( 'POLYGON((' ||ST_X(shape)||' '||ST_Y(shape)||',' ||ST_X(shape)+velocity*cos((direction+15)*PI()/180)/20||' '||ST_Y(shape)+velocity*sin((direction+15)*PI()/180)/20||',' ||ST_X(shape)+velocity*cos((direction-15)*PI()/180)/20||' '||ST_Y(shape)+velocity*sin((direction-15)*PI()/180)/20||',' ||ST_X(shape)||' '||ST_Y(shape)||'))' ) , velocity/50 ) ) from sde.wind
參考:
