PostGIS空间查询学习笔记

1. 开始

1.1 启用PostGIS扩展

官方文档：Getting Started

-- 为数据库启用postgis扩展
CREATE EXTENSION postgis;

1.2 创建空间数据表

官方文档：Spatial Tables

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

文档：5.3. Geometry Constructors

• 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

  -- 在创建表时
  ALTER TABLE mytable SET (fillfactor = 70)
  -- 已经插入了数据的情况下，执行上述语句只会影响以后插入的数据。
  -- 可以使用VACUUM (FULL)或CLUSTER重写表。
  

  fillfactor是介于10到100之间的值（默认100），它确定INSERT将填充表块的百分比。
• 更新结果直接插入新表
  连续写比间隔写要快。最后选择了这种方法，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
  在创建时直接指定SRID

  SELECT 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_Transform

  SELECT 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