postgis拓扑分析避坑-一一网

需求分析

项目中需要判断几何图形之间是否重叠，并剔除共边关系，需要用到postgis的空间分析函数，查看官方文档，得到以下有用的函数：

函数名	函数解释
ST_Intersects	如果一个几何体或地理体共享空间的任何部分，那么它们就会相交。对于地理学来说，公差是0.00001米（所以任何接近的点都被认为是相交的）。 ST_Overlaps, ST_Touches, ST_Within都意味着空间上的相交。如果上述任何一个返回真，那么几何体也是空间相交的。
ST_Overlaps	如果几何体共享空间，具有相同的维度，但不完全被对方包含，则返回TRUE。如果几何之间存在”空间重叠”，返回TRUE。我们的意思是它们相交，但一个并不完全包含另一个。这个函数调用将自动包括一个包围盒比较，它将利用几何体上的任何可用索引。要避免使用索引，请使用函数_ST_Overlaps。
ST_Touches	如果这些几何体至少有一个共同点，但它们的内部不相交，则返回TRUE。如果g1和g2之间唯一的共同点位于g1和g2的边界的结合处，返回TRUE。 ST_Touches关系适用于所有区域/区域、线/线、线/区域、点/区域和点/线对的关系，但不适用于点/点对。

对于以上函数，我们可以得到，如果两个几何之间”ST_Intersects“为True，并且“ST_Touches ”为False，则证明几何之间存在”重叠“关系，不存在”共边“关系。

项目使用框架

框架/库	版本
postgresql	11
postgis	3.0.0 r17983
python	2.7.8 win32
Flask-SQLAlchemy	2.4.1
GeoAlchemy2	0.6.3

例子

假设有两个空间图形geojson1和geojson2，他们的geojson分别是：

{"type":"MultiPolygon","coordinates":[[[[111.21124454,21.5108185],[111.21184285,21.51176974],[111.21185173,21.51178386],[111.21192795,21.51176727],[111.21195947,21.51176041],[111.21221703,21.51168379],[111.21229219,21.5116613],[111.21237885,21.51164951],[111.21241984,21.51163725],[111.21253281,21.51162187],[111.21267234,21.51160124],[111.21296984,21.51151721],[111.21311642,21.51148691],[111.21322617,21.51144904],[111.21333346,21.51137632],[111.21376959,21.51108607],[111.21389188,21.51100818],[111.21388765,21.51096566],[111.2138591,21.51077055],[111.21397331,21.51068964],[111.21420174,21.51057543],[111.21462052,21.51033272],[111.21484419,21.51021851],[111.21508214,21.51011381],[111.21526774,21.50996629],[111.21539147,21.5098378],[111.2155961,21.50979972],[111.21591971,21.509814],[111.21612434,21.50980448],[111.21621476,21.5097902],[111.21636705,21.5096284],[111.21644795,21.50915251],[111.21637656,21.50876228],[111.21644639,21.50811051],[111.21646222,21.50796278],[111.21640512,21.50782001],[111.21628139,21.50766297],[111.21613862,21.50763441],[111.21578646,21.50767725],[111.21547713,21.50777718],[111.21472046,21.50807224],[111.21392572,21.50839108],[111.21361163,21.50857192],[111.21331182,21.50867662],[111.2129787,21.50876704],[111.21262654,21.50895264],[111.21241238,21.50907637],[111.21206974,21.50953323],[111.21143205,21.51041839],[111.21130832,21.51060398],[111.21124454,21.5108185]]]]}
复制代码

{"type":"MultiPolygon","coordinates":[[[[111.21640512,21.50782001],[111.21646222,21.50796278],[111.21644639,21.50811051],[111.21637656,21.50876228],[111.21644795,21.50915251],[111.21636705,21.5096284],[111.21648874,21.50984016],[111.21717878,21.50960221],[111.21727396,21.50937616],[111.2171074,21.50749639],[111.21640512,21.50782001]]]]}
复制代码

直接使用sql语句进行查询，将会得到以下结果：

ST_Touches(
(SELECT st_geomfromgeojson('geojson1')  ),
(SELECT st_geomfromgeojson('geojson2') )
),
st_intersects(
(SELECT st_geomfromgeojson('geojson1')  ),
(SELECT st_geomfromgeojson('geojson2') )
),
st_overlaps(
(SELECT st_geomfromgeojson('geojson1')  ),
(SELECT st_geomfromgeojson('geojson2') )
);

复制代码

st_touches	st_intersects	st_overlaps
t	t	f

可以看到，这两个几何图形是共边关系，并不存在重叠关系。

问题

对于某种情况下，我们使用以上的拓扑分析函数会出现意料之外的结果。
当使用geoalchemy-2和sqlalchemy创建multipolygon类型字段之后，建立GIST索引，分别将geojson1和geojson2的存入表中，直接使用以下的sql语句会发生问题：

SELECT ST_Touches(
(SELECT feature.multipolygon from feature where objectid = 1),
(SELECT feature.multipolygon from feature where objectid = 2)
),st_intersects(
(SELECT feature.multipolygon from feature where objectid = 1),
(SELECT feature.multipolygon from feature where objectid = 2)
),
st_overlaps(
(SELECT feature.multipolygon from feature where objectid = 1),
(SELECT feature.multipolygon from feature where objectid = 2)
)
UNION ALL
SELECT _ST_Touches(
(SELECT feature.multipolygon from feature where objectid = 1),
(SELECT feature.multipolygon from feature where objectid = 2)
),_st_intersects(
(SELECT feature.multipolygon from feature where objectid = 1),
(SELECT feature.multipolygon from feature where objectid = 2)
),
_st_overlaps(
(SELECT feature.multipolygon from feature where objectid = 1),
(SELECT feature.multipolygon from feature where objectid = 2)
)
复制代码

st_touches	st_intersects	st_overlaps
f	t	t
f	t	t

根据以上查询结果发现，本来两个几何图形之间应该是共边关系，但是存入空间数据库再进行拓扑分析变成了重叠关系，无论是否使用索引查询都是同样的结果。

解决方法

由于问题来源不明，可能是postgis或GeoAlchemy2 版本导致，因此暂定在sql中将几何转为geojson，再进行拓扑分析，可以得到正确的结果：

SELECT ST_Touches(
(SELECT st_geomfromgeojson(st_asgeojson(feature.multipolygon)) from feature where objectid = 1),
(SELECT st_geomfromgeojson(st_asgeojson(feature.multipolygon)) from feature where objectid = 2)
),st_intersects(
(SELECT st_geomfromgeojson(st_asgeojson(feature.multipolygon)) from feature where objectid = 1),
(SELECT st_geomfromgeojson(st_asgeojson(feature.multipolygon)) from feature where objectid = 2)
),
st_overlaps(
(SELECT st_geomfromgeojson(st_asgeojson(feature.multipolygon)) from feature where objectid = 1),
(SELECT st_geomfromgeojson(st_asgeojson(feature.multipolygon)) from feature where objectid = 2)
)
复制代码