PostgreSQL+PostGIS 的使用

標籤：

一、PostGIS中的幾何類型

PostGIS支援所有OGC規範的“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中的常用函數

以下內容包括比較多的角括弧，發布到blogger的時候會顯示不正常，內容太多我也無暇一個個手動改代碼，因此如有問題就去參考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)

幾何物件建構函數 ：
參考語義：
Text：WKT
WKB：WKB
Geom:Geometry
M:Multi
Bd:BuildArea
Coll:Collection ST_GeomFromText(text,[])

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)

四、 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

下面是PostGIS進行運算後返回的結果：

在這裡，Flex應用與伺服器的互動通過BlazeDS進行,下面是本樣本在伺服器端的Java代碼：

package wuyf;import java.sql.Connection;import java.sql.DriverManager;import java.sql.ResultSet;import java.sql.Statement;import java.util.ArrayList;import java.util.HashMap;public class Wind{private Connection conn = null;public Connection getConn(){if (conn==null){try{Class.forName("org.postgresql.Driver");String url = "jdbc:postgresql://localhost:5432/sde" ;conn = DriverManager.getConnection(url, "sde" , "pwd" );conn.setAutoCommit(false);}catch(Exception e){System.err.print(e);}}return conn;}public ArrayList > getWinds(){ArrayList > result = new ArrayList >();if ( this.getConn()==null )return result;try{String sql = "select *,ST_AsGeoJson(shape) from sde.wind";Statement st = this.getConn().createStatement();st.setFetchSize(0);ResultSet rs = st.executeQuery(sql);while (rs.next()){HashMap map = new HashMap ();map.put("shape", rs.getString("ST_AsGeoJson"));map.put("velocity", rs.getString("velocity"));map.put("direction", rs.getString("direction"));result.add(map);}rs.close();st.close();}catch(Exception e){System.err.print(e);}return result;} public ArrayList > getEffectZones(){ArrayList > result = new ArrayList >(); if ( this.getConn()==null )return result;try{String sql = "select *,ST_AsGeoJson(";sql+= "ST_Buffer(";sql+= "ST_PolygonFromText(";sql+= "‘POLYGON((‘";sql+= "||ST_X(shape)||‘ ‘||ST_Y(shape)||‘,‘";sql+= "||ST_X(shape)+velocity*cos((direction+15)*PI()/180)/20||‘ ‘||ST_Y(shape)+velocity*sin((direction+15)*PI()/180)/20||‘,‘";sql+= "||ST_X(shape)+velocity*cos((direction-15)*PI()/180)/20||‘ ‘||ST_Y(shape)+velocity*sin((direction-15)*PI()/180)/20||‘,‘";sql+= "||ST_X(shape)||‘ ‘||ST_Y(shape)||‘))‘";sql+= ")";sql+= ", velocity/50";sql+= ")";sql+= ") ";sql+="from sde.wind";Statement st = this.getConn().createStatement();st.setFetchSize(0);ResultSet rs = st.executeQuery(sql);while (rs.next()){HashMap map = new HashMap ();map.put("shape", rs.getString("ST_AsGeoJson"));map.put("velocity", rs.getString("velocity"));map.put("direction", rs.getString("direction"));result.add(map);}rs.close();st.close();}catch(Exception e){System.err.print(e);}return result;}}

 

PostgreSQL+PostGIS 的使用