1. 开始
1.1 启用PostGIS扩展
-- 为数据库启用postgis扩展
CREATE EXTENSION postgis;
1.2 创建空间数据表
CREATE TABLE roads (
id SERIAL PRIMARY KEY,
name VARCHAR(64),
geom geometry(LINESTRING,3005) -- 参数:几何类型,SRID
);
2. 常用函数
2.1 ST_Boundary
文档:4.2.1.1. Dimensionally Extended 9-Intersection Model
在九交模型的部分给出了boundary的定义:
对于POINT:空集; 对于LINESTRING:两个端点;POLYGON:内外环。
2.2 构造几何类型
WKT
文档:5.8.1. Well-Known Text (WKT)
- ST_GeomFromText()构造函数,将WKT转换为geometry类型。可以指定SRID。
Constructor
- ST_Point(float x, float y) / ST_MakePoint(float x, float y, float z, float m):
ST_MakePoint比ST_GeomFromText、ST_PointFromText更快更准确,并且更适合数值类型。(以前自己拼WKT字符串真是往事不堪回首)
这两个函数没有坐标系参数,默认坐标系为SRID=0。如果需要指定坐标系,在外层套一个函数:ST_SetSRID(),如:SELECT ST_SetSRID( ST_Point( -71.104, 42.315), 4326)
- ST_Polygon(geometry lineString, integer srid)
- ST_GeomFromEWKT()
2.3 ST_Buffer
文档:ST_Buffer
签名:geometry ST_Buffer(geometry g1, float radius_of_buffer, text buffer_style_parameters = ‘’);
最后一个参数可以控制buffer的形状(拐点、端点形状、生成缓冲区的方向等)
2.4 ST_Simplify
文档: ST_Simplify
使用Douglas-Peucker算法简化线、多边形。速度非常快,对一张有4万个线要素的表只需要不到2s。
3. 关于空间索引
3.1 如何利用空间索引
文档:4.2.2. Taking Advantage of Indexes
以下为部分重要内容的翻译:
与空间查询有关的函数或操作符,必须被包含在 * WHERE* 或 * OR * 子句中。空间操作符包括&&、<->等,函数包括ST_Intersects、ST_Dwithin等。
注意:ST_Distance没有经过空间索引的优化,效率很低,请使用ST_Dwithin。(ST_Dwithin隐式地使用了操作符&&,基于bounding-box先做了筛选,再进行距离计算。
4. 空间查询
** 针对pg的语法
** 用于举例的表:
一个多边形表
字段 | 描述 |
---|---|
id | 编号 |
geom | 多边形的几何图形 |
一个点表
字段 | 描述 |
---|---|
id | 编号 |
geom | 点的几何图形 |
poly_id | 所在的多边形的id |
4.1 空间连接(Spatial join)
我的理解:把普通的join里面on Table1.column1 = Table2.column2替换成返回值为bool类型的空间查询。
查询落在多边形里的点以及对应的多边形编号:参考了这个回答
-- 使用Inner join
SELECT
poly.id,
pt.id
FROM
[poly_table_name] poly
INNER JOIN
[pt_table_name] pt
ON
St_Intersects(poly.geom, pt.geom)
-- 或者
SELECT
poly.id,
pt.id
FROM
[poly_table_name] poly, [pt_table_name] pt
WHERE
St_Intersects(poly.geom, pt.geom)
4.2 多表联合更新(Update)
4.2.1 常规方法
根据空间关系批量更新属性。在点表中更新poly_id一列,使其的值为点所在的多边形的id。(也可以遍历多边形表,使用子查询,就是有点麻烦)
参考1 参考2
注:在from里使用join似乎不是postgres的语法
update [pt_table_name] pt
set poly_id = poly.id
from [poly_table_name] poly
where st_intersects(poly.geom, pt.geom)
4.2.2 对大量数据进行更新
对于数据量较小的情况,上述方法就足够了(eg. 9.8w行16列的表,整个表进行更新,大约花费4min)。然而数据量大的时候,此方法慢到无法接受(eg. 更新250w行的表,10h过去了都没有结束……)
可能的原因
- PG中update的机制是delete+insert,即将旧数据隐藏并插入新数据,旧数据并没有删除,导致表膨胀(在pgAdmin中查看表的statistics,可以发现有大量死元组),需要用vacuum命令进行回收。
- 更新的列上有索引或触发器。
- 更新的表上有死锁。
参考1
改进
- HOT(Heap Only Tuple) Update
HOT介绍 简而言之,在每个表块上预留一定的空间给新数据,这样可以用一个Hot link指向新数据,不用重建索引。(所以主要是针对不想重建索引的情况?)
使用HOT
fillfactor是介于10到100之间的值(默认100),它确定INSERT将填充表块的百分比。-- 在创建表时 ALTER TABLE mytable SET (fillfactor = 70) -- 已经插入了数据的情况下,执行上述语句只会影响以后插入的数据。 -- 可以使用VACUUM (FULL)或CLUSTER重写表。
- 更新结果直接插入新表
连续写比间隔写要快。最后选择了这种方法,170w左右的表大约需要30min。SELECT pt.*, poly.id AS poly_id INTO [new_pt_table] FROM [poly_table_name] poly RIGHT JOIN [pt_table_name] pt ON St_Intersects(poly.geom, pt.geom)
5. 关于SRID
SRID是坐标系的编号,在spatial_ref_sys表中:
5.1 指定SRID
指定坐标系,不改变坐标数值,空间位置发生变化。假设SRID2的原点在SRID中的(1,1)位置。
原本SRID=SRID1,某点的坐标为(1, 1);指定SRID=SRID2,点的坐标仍为(1, 1),但其在SRID1中的位置变成了(2, 2)。
相关的函数有:
- ST_SetSRID
文档:ST_SetSRID - ST_UpdateGeometrySRID
文档:UpdateGeometrySRID
修改一个列中全部geometry的SRID。SELECT UpdateGeometrySRID(table_name, geom_col, new_srid)
- ST_GeomFromEWKT(text EWKT)
文档:ST_GeomFromEWKT
在创建时直接指定SRIDSELECT ST_GeomFromEWKT('SRID=4269;LINESTRING(-71.160281 42.258729,-71.160837 42.259113,-71.161144 42.25932)'); SELECT ST_GeomFromEWKT('SRID=4269;POINT(-71.064544 42.28787)');
5.2 坐标转换
根据前后两个坐标系的数学关系进行坐标转换,参考系变化导致坐标数值变化,但空间位置不变。
原本SRID=SRID1,点的坐标为(1, 1);转换到SRID2坐标系下,新坐标为(0, 0),在原坐标系中的位置没有变化。
相关的函数有:
- ST_Transform
文档:ST_TransformSELECT ST_AsText( ST_Transform( ST_GeomFromText( 'POLYGON((743238 2967416,743238 2967450, 743265 2967450,743265.625 2967416,743238 2967416))', 2249), 4326) ) As wgs_geom; wgs_geom --------------------------- POLYGON((-71.1776848522251 42.3902896512902,-71.1776843766326 42.3903829478009, -71.1775844305465 42.3903826677917,-71.1775825927231 42.3902893647987,-71.177684 8522251 42.3902896512902));
6. 拓扑检查
6.1 重叠
postgis原生的拓扑检查功能st_overlaps(部分重叠)、st_equals(完全重叠)要求非常严格,而实际上由于浮点数舍入误差,从认知上可以判断为是重叠的部分用这些方法无法检测出来。
检测road表(linestring)是否有相交的方法(考虑容差)
select a.gid, b.gid
from road a, road b
where
a.geom && b.geom and -- 使用空间索引
a.gid != b.gid and -- 排除自身
st_length(st_intersection(a.geom, st_buffer(b.geom, 0.001)))>0.1
其中a.geom && b.geom强制使用空间索引,不加这句的执行计划:
加上之后(注意Index Cond):
实测执行时间,前者10min,后者10s。
6.2 相交(ST_Relate)
官方文档:4.2.1.3. General Spatial Relationships
官方文档:ST_Relate
关于九交模型的一篇博客
6.3 自相交(ST_IsSimple)
文档:ST_IsSimple
检测是否有自相交(self intersection)或自切(self tangency),可用于拓扑检查。(处理方法:大多数情况下,和自身求intersection,再linemerge,可以修复)
和IsValid的关系:OGC Validity
IsSimple针对0、1维图形,而IsValid针对2维图形(Polygon、MultiPolygon)
类型 | Simple | Valid |
---|---|---|
Point | 一定满足 | \ |
MultiPoint | 没有重叠的点 | \ |
LineString | 不经过同一点两次(除非是端点 | \ |
MultiLinestring | 所有构成部分都Simple,且子集的相交只在结点处 | \ |
Polygon | \ | 没有Ring交叉(可以有相切的点但不能有重合的线) |
MultiPolygon | \ | 所有构成部分都Valid,且没有子集互相相交 |
7 补充一些时间运算
7.1 时间格式
- timestamp [without time zone]
- timestamp with time zone
- date:日期
- time:time of day(时间点)
- interval:时间间隔(时间作差得到)
7.2 时间运算
- 时间加减:时间字段直接加减字符串
SELECT '2022-01-01'::timestamp + '30 min'
Result: 2022-01-01 00:30:00
- 获取timestamp、interval的sub field:date_part和extract
SELECT date_part('hour', '2022-01-01 20:30'::timestamp)
Result: 20
SELECT extract('hour' from '2022-01-01 20:30'::timestamp)
Result: 20
SELECT date_part('hour', '80 mins'::interval)
Result: 1
SELECT extract(hours from '80 minutes'::interval)
Result:1
-- 周几
SELECT EXTRACT(DOW FROM TIMESTAMP '2001-02-16 20:38:40')
Result: 5
- 获取以秒计的时间差。注意,如果这里用extract或date_part,只能提取时间差的“秒”字段(0~59)。
select extract('epoch' from '2020-01-02'::timestamp - '2020-01-01'::timestamp)
Result: 86400