一、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 2),(2 3,3 2,5 4)) ——多线
- MULTIPOLYGON(((0 0,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 0) ——插值圆弧
COMPOUNDCURVE(CIRCULARSTRING(0 0,1 0),(1 0,0 1)) ——插值复合曲线
CURVEPOLYGON(CIRCULARSTRING(0 0,3 3,3 1,1 1)) ——曲线多边形
MULTICURVE((0 0,5 5),CIRCULARSTRING(4 0,8 4)) ——多曲线
MULTISURFACE(CURVEPOLYGON(CIRCULARSTRING(0 0,((10 10,14 12,11 10,10 10),(11 11,11.5 11,11 11.5,11 11))) ——多曲面
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],"6326"]],PRIMEM["Greenwich","8901"]],UNIT["degree",0.01745329251994328,"9122"]],"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)',sans-serif; font-size:14px; orphans:2; widows:2"> 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);
另外要注意的是,空间索引只有在进行基于边界范围的查询时才起作用,比如“&&”操作。
以下内容包括比较多的尖括号,发布到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,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,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)
下面我们通过一个简单的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 >(); ; 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 = return result; } }