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_Simplify	使用Douglas-Peucker算法返回给定几何体的 "简化 "版本。 实际上只对（多）线和（多）多边形进行处理，但你可以安全地对任何种类的几何体调用它。 由于简化是在逐个对象的基础上进行的，你也可以给这个函数输入一个几何体集合。 "保留塌陷 "标志将保留那些在公差条件下会太小的对象。 例如，一条1米长的线以10米的公差简化。 如果给了保留标志，这条线就不会消失。 这个标志对渲染引擎很有用，可以避免大量非常小的物体从地图上消失，留下令人惊讶的空隙。
ST_Snap	将一个几何体的顶点和线段抓取到另一个几何体的顶点上。 抓取距离的公差用于控制抓取的位置。 结果几何体是输入几何体的顶点被抓取后的结果。 如果没有抢夺发生，那么输入的几何体将被返回，没有任何变化。 将一个几何体钉在另一个几何体上，可以通过消除几乎重合的边缘 （在编码和交集计算过程中会导致问题）来提高叠加操作的稳健性。 过多的捕捉会导致无效的拓扑结构被创建，所以捕捉顶点的数量和位置是通过启发式方法来决定何时捕捉是安全的。 然而，这可能会导致一些潜在的抢占被省略掉。

对于以上函数，我们可以得到，如果两个几何之间”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和ArcEngine在空间查询的处理上是有偏差的，ArcEngine加上了容差而PostGIS应该只是在相交的判断上加了容差。并且PostGIS中的 ST_Intersects默认参数是0.00001米，不能更改。 根据stackflow上的提示，如果要设置默认容差，可以在数据库添加如下函数：

CREATE OR REPLACE FUNCTION ST_Intersects(geography, geography,float8)
    RETURNS boolean
    AS 'SELECT $1 && $2 AND _ST_Distance($1, $2, 0.0, false) < $3'
    LANGUAGE 'sql' IMMUTABLE ;

解决方法

方法1

在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)
)

st_touches	st_intersects	st_overlaps
t	t	f

方法2

通过方法1可以结局部分几何之间的拓扑关系错误问题，但是对于某些几何图形如下，则不能得到正确的关系：

{"type":"MultiPolygon","coordinates":[[[[111.0723406,21.5049103],[111.07233446,21.50490513],[111.07222508,21.50475759],[111.07203,21.50430971],[111.07201278,21.5042166],[111.07199821,21.50413781],[111.07196294,21.50380627],[111.07198383,21.50358818],[111.07199469,21.5034747],[111.07216043,21.50298178],[111.07222792,21.5028999],[111.07227687,21.5028405],[111.0721827,21.50283679],[111.072163,21.50283601],[111.07200164,21.50282964],[111.07196158,21.50282119],[111.07161041,21.50272766],[111.07126649,21.50252217],[111.07081693,21.50209146],[111.07071911,21.50193103],[111.07018377,21.50105303],[111.06959261,21.4999927],[111.06939792,21.4996435],[111.069187,21.49927205],[111.0689942,21.49882584],[111.0689278,21.49867218],[111.06891811,21.49864401],[111.06889065,21.4985061],[111.06891686,21.49822628],[111.06892447,21.49810298],[111.06887527,21.4977503],[111.06877668,21.49728526],[111.06867775,21.49710206],[111.06853058,21.49699869],[111.06849433,21.49697322],[111.06842451,21.49700493],[111.06826067,21.49718995],[111.06800183,21.4973277],[111.06785324,21.49736503],[111.0676817,21.49747264],[111.06738154,21.49771649],[111.0669772,21.49794749],[111.06638289,21.49820724],[111.06596621,21.49827766],[111.06557092,21.49820899],[111.0651526,21.49804141],[111.06509439,21.49803426],[111.06499081,21.49807019],[111.06495233,21.49809976],[111.06487276,21.49820586],[111.06482882,21.49862675],[111.06442551,21.499164],[111.06480719,21.49996883],[111.06497851,21.50019726],[111.06533543,21.50048279],[111.0660564,21.501061],[111.06639771,21.50145521],[111.06674533,21.50185672],[111.06682716,21.50211476],[111.06696776,21.50245422],[111.06700431,21.50254245],[111.06744565,21.50334276],[111.06856891,21.50366411],[111.0693734,21.50399354],[111.06960711,21.50412606],[111.06997871,21.50426928],[111.07178306,21.50493163],[111.07199759,21.50501039],[111.07221207,21.50498162],[111.0723406,21.5049103]]]]}

{"type":"MultiPolygon","coordinates":[[[[111.0723406,21.5049103],[111.07221207,21.50498162],[111.07199759,21.50501039],[111.07178306,21.50493163],[111.06997871,21.50426928],[111.06960711,21.50412606],[111.0693734,21.50399354],[111.06856891,21.50366411],[111.06744565,21.50334276],[111.06736479,21.50348293],[111.06724553,21.50364625],[111.06698112,21.50376352],[111.06681264,21.50388306],[111.06676165,21.50405513],[111.06686264,21.50429273],[111.06668303,21.50461874],[111.06648309,21.50518754],[111.06638725,21.50555749],[111.06637814,21.50573432],[111.06641093,21.50597327],[111.0664285,21.50641868],[111.06647659,21.50661938],[111.06664221,21.50693856],[111.06688954,21.50727205],[111.06690897,21.50731601],[111.06693152,21.50749761],[111.06697444,21.5075387],[111.06708649,21.50763182],[111.06713997,21.50767627],[111.06721965,21.507703],[111.06729902,21.5076846],[111.06731463,21.50768098],[111.06743209,21.50767137],[111.06747602,21.50768952],[111.06751282,21.50774207],[111.06770396,21.50804213],[111.06773769,21.50809086],[111.06785444,21.50825057],[111.06814327,21.50885501],[111.06850504,21.50952327],[111.06872551,21.50993052],[111.06878477,21.51003999],[111.06894204,21.51013331],[111.06913484,21.50993983],[111.06921857,21.50989487],[111.06932225,21.50984064],[111.06933172,21.50983372],[111.0693588,21.50982153],[111.06989211,21.50958139],[111.06996352,21.5095603],[111.0700009,21.50954927],[111.07024062,21.50947849],[111.07072184,21.50922021],[111.07125417,21.50903838],[111.07160774,21.50884706],[111.0718273,21.50875019],[111.07226134,21.50860901],[111.0727081,21.50838181],[111.07312861,21.50820371],[111.07321359,21.50816772],[111.07326952,21.50814916],[111.07330527,21.50814914],[111.07334008,21.50815529],[111.07331602,21.50813069],[111.07327309,21.50806475],[111.07324632,21.50785974],[111.07324632,21.50783884],[111.07324632,21.5077349],[111.07332396,21.50744451],[111.07334137,21.5073794],[111.0735869,21.50688346],[111.07411805,21.5062581],[111.07405529,21.50615429],[111.07398575,21.50611467],[111.0739737,21.5061078],[111.07290856,21.50538889],[111.0723406,21.5049103]]]]}

ST_Touches(
(SELECT st_geomfromgeojson('geojson3')  ),
(SELECT st_geomfromgeojson('geojson4') )
),
st_intersects(
(SELECT st_geomfromgeojson('geojson3')  ),
(SELECT st_geomfromgeojson('geojson4') )
),
st_overlaps(
(SELECT st_geomfromgeojson('geojson3')  ),
(SELECT st_geomfromgeojson('geojson4') )
);

st_touches	st_intersects	st_overlaps
f	t	t

输入几何图形之间并无“重叠”关系，使用st_overlaps仍然得到True。 所以先使用st_intersects进行空间关系的判断，再使用shapely判断几何之间是否存在”共边“关系。

方法3 最终解决方法

弧长公式：L=nπr/180°或l=|α|r 地球半径大致是6400千米 以纬度0.000001为例： 弧长＝（0.000001/180）×3.14×6400 = 0.000111644444千米 约等于0.1米

使用ST_Simplify在一定精度下对图斑进行简化，再在在一定精度下使用ST_Snap使两个相互对比的几何之间相邻的坐标点抓取到同一个位置，之后使用拓扑检查函数进行检查

# 0.1米精度
tolerance =  0.1 * 180 / (6356752.31414 * 3.14)

st_overlaps (
	( SELECT ST_Simplify ( multipolygon, tolerance  ) ),
	(
	SELECT
		st_snap (
			ST_Simplify ( multipolygon1, tolerance ),
			( SELECT ST_Simplify ( multipolygon2,tolerance  ) ),
			tolerance 
		) 
	) 
)