Demand analysis

In the project, it is necessary to judge whether geometric shapes overlap and remove common edge relations. Postgis spatial analysis functions are needed, and the following useful functions are obtained by checking official documents:

The function name Function to explain
ST_Intersects If a geometry or geography shares any part of space, they will intersect.

For geography, the tolerance is 0.00001 m (so any close points are considered to intersect).

ST_Overlaps, ST_Touches, ST_Within all mean the intersection of Spaces. If any of the above returns true, then the geometry is also spatially intersected.
ST_Overlaps Returns TRUE if geometers share space and have the same dimensions, but are not completely contained by each other.

Return TRUE if there is “spatial overlap” between geometries. We mean they intersect, but one doesn’t completely contain the other.

This function call will automatically include a bounding box comparison that will take advantage of any indexes available on the geometry. To avoid using indexes, use the function _ST_Overlaps.
ST_Touches Returns TRUE if these geometries have at least one thing in common, but their interiors do not intersect.

Returns TRUE if the only common point between G1 and G2 is at the junction of g1 and G2 boundaries.

ST_Touches applies to all area/area, line/line, line/area, point/area and point/line pairs, but not point/point pairs.
ST_Simplify Returns a “simplified” version of the given geometry using the Douglas-Peucker algorithm.

It really only works with (multi) lines and (multi) polygons, but you can safely call it on any kind of geometry.

Since the simplification is done on an object by object basis, you can also enter a collection of geometry into the function.

The “retention collapse” sign will preserve objects that are too small for tolerances.

For example, a 1-meter line is simplified with a tolerance of 10 meters.

If the hold flag is given, the line will not disappear.

This flag is useful for rendering engines to prevent a large number of very small objects from disappearing from the map, leaving surprising gaps.
ST_Snap Grab the vertices and line segments of one geometry onto the vertices of another.

The tolerance of grasping distance is used to control the position of grasping.

The resulting geometry is the result of fetching the vertices of the input geometry.

If no looting occurs, the input geometry is returned unchanged.

Nailing one geometry to another can be done by eliminating edges that almost overlap

(causes problems during coding and intersection calculations) to improve the robustness of the overlay operation.

Too many captures can cause invalid topologies to be created, so the number and locations of captured vertices are heuristically used to determine when it is safe to capture.

However, this may result in some potential preemption being omitted.

As for the above function, if “ST_Intersects” is True and “ST_Touches” is False, then the geometry does overlap but does not share edges.

Project Usage framework

Frame/library version
postgresql 11
postgis 3.0.0 r17983
python 2.7.8 win32
Flask-SQLAlchemy Against 2.4.1
GeoAlchemy2 0.6.3

example

Suppose we have two spatial figures, Geojson1 and Geojson2, and their GEOJSON is:

{" 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.21389 188,21.51100818], [111.21388765, 21.51096566], [111.2138591, 21.51077055], [111.21397331, 21.51068964], [111.21420174, 21.510575 43], [111.21462052, 21.51033272], [111.21484419, 21.51021851], [111.21508214, 21.51011381], [111.21526774, 21.50996629], [111.215 39147,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.216462 22,21.50796278], [111.21640512, 21.50782001], [111.21628139, 21.50766297], [111.21613862, 21.50763441], [111.21578646, 21.507677 25], [111.21547713, 21.50777718], [111.21472046, 21.50807224], [111.21392572, 21.50839108], [111.21361163, 21.50857192], [111.213 31182,21.50867662], [111.2129787, 21.50876704], [111.21262654, 21.50895264], [111.21241238, 21.50907637], [111.21206974, 21.5095 3323], [111.21143205, 21.51041839], [111.21130832, 21.51060398], [111.21124454, 21.5108185]]]]}Copy the code
{" 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.2171787 8,21.50960221], [111.21727396, 21.50937616], [111.2171074, 21.50749639], [111.21640512, 21.50782001]]]]}Copy the code

Query directly using SQL statements, the following results will be obtained:

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') )
);

Copy the code
st_touches st_intersects st_overlaps
t t f

As you can see, these two geometry are coeding-edge and do not overlap.

The problem

In some cases, we use the above topological analysis functions to produce unexpected results. SQL > create table geojson1, geojson2, geojson2, geojson2

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) )Copy the code
st_touches st_intersects st_overlaps
f t t
f t t

According to the above query results, it is found that the relationship between the two geometric shapes should be common edge, but the relationship becomes overlapping when the spatial database is stored and topology analysis is carried out. The results are the same regardless of whether the index query is used.

why

PostGIS and ArcEngine are biased in the processing of spatial queries, ArcEngine adds tolerance and PostGIS should only add tolerance to the judgment of intersection. In addition, the default parameter of ST_Intersects in PostGIS is 0.00001 meters and cannot be changed. Following the tips on StackFlow, if you want to set the default tolerance, you can add the following function to the database:

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

The solution

Method 1

By converting geometry into GEOJSON in SQL and then conducting topological analysis, correct results can be obtained for the above two geometric relations:

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) )Copy the code
st_touches st_intersects st_overlaps
t t f

Method 2

The problem of incorrect topological relations between partial geometries can be solved by method 1. However, for some geometries as follows, correct relations cannot be obtained:

{" 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.502 72766], [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.498672 18], [111.06891811, 21.49864401], [111.06889065, 21.4985061], [111.06891686, 21.49822628], [111.06892447, 21.49810298], [111.0688 7527,21.4977503], [111.06877668, 21.49728526], [111.06867775, 21.49710206], [111.06853058, 21.49699869], [111.06849433, 21.49697 322], [111.06842451, 21.49700493], [111.06826067, 21.49718995], [111.06800183, 21.4973277], [111.06785324, 21.49736503], [111.067 6817,21.49747264], [111.06738154, 21.49771649], [111.0669772, 21.49794749], [111.06638289, 21.49820724], [111.06596621, 21.49827 766], [111.06557092, 21.49820899], [111.0651526, 21.49804141], [111.06509439, 21.49803426], [111.06499081, 21.49807019], [111.064 95233,21.49809976], [111.06487276, 21.49820586], [111.06482882, 21.49862675], [111.06442551, 21.499164], [111.06480719, 21.49996 883], [111.06497851, 21.50019726], [111.06533543, 21.50048279], [111.0660564, 21.501061], [111.06639771, 21.50145521], [111.06674 533,21.50185672], [111.06682716, 21.50211476], [111.06696776, 21.50245422], [111.06700431, 21.50254245], [111.06744565, 21.50334 276], [111.06856891, 21.50366411], [111.0693734, 21.50399354], [111.06960711, 21.50412606], [111.06997871, 21.50426928], [111.071 78306,21.50493163], [111.07199759, 21.50501039], [111.07221207, 21.50498162], [111.0723406, 21.5049103]]]]}Copy the code
{" 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.066483 09,21.50518754], [111.06638725, 21.50555749], [111.06637814, 21.50573432], [111.06641093, 21.50597327], [111.0664285, 21.5064186 8], [111.06647659, 21.50661938], [111.06664221, 21.50693856], [111.06688954, 21.50727205], [111.06690897, 21.50731601], [111.0669 3152,21.50749761], [111.06697444, 21.5075387], [111.06708649, 21.50763182], [111.06713997, 21.50767627], [111.06721965, 21.50770 3], [111.06729902, 21.5076846], [111.06731463, 21.50768098], [111.06743209, 21.50767137], [111.06747602, 21.50768952], [111.06751 282,21.50774207], [111.06770396, 21.50804213], [111.06773769, 21.50809086], [111.06785444, 21.50825057], [111.06814327, 21.50885 501], [111.06850504, 21.50952327], [111.06872551, 21.50993052], [111.06878477, 21.51003999], [111.06894204, 21.51013331], [111.06 913484,21.50993983], [111.06921857, 21.50989487], [111.06932225, 21.50984064], [111.06933172, 21.50983372], [111.0693588, 21.509 82153], [111.06989211, 21.50958139], [111.06996352, 21.5095603], [111.0700009, 21.50954927], [111.07024062, 21.50947849], [111.07 072184,21.50922021], [111.07125417, 21.50903838], [111.07160774, 21.50884706], [111.0718273, 21.50875019], [111.07226134, 21.508 60901], [111.0727081, 21.50838181], [111.07312861, 21.50820371], [111.07321359, 21.50816772], [111.07326952, 21.50814916], [111.0 7330527,21.50814914], [111.07334008, 21.50815529], [111.07331602, 21.50813069], [111.07327309, 21.50806475], [111.07324632, 21.5 0785974], [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.506 1078], [111.07290856, 21.50538889], [111.0723406, 21.5049103]]]]}Copy the code
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') )
);

Copy the code
st_touches st_intersects st_overlaps
f t t

There is no “overlap” between the input geometries. Use st_OVERLaps and still get True. Therefore, ST_Intersects was first used to judge the spatial relationship, and then Shapely was used to judge whether there was a “common edge” relationship between geometries.

Method 3 Final solution

Arc length formula: L = n PI r / 180 ° or L = | alpha | r radius of the earth is about 6400 kilometers Latitude 0.000001, for example, arc length = (0.000001/180) * 3.14 * 6400 = 0.000111644444 km Approximately equal to 0.1 m

ST_Simplify was used to simplify the map spots under a certain accuracy, and ST_Snap was used to make the adjacent coordinate points between two contrasting geometers capture the same position under a certain accuracy, and then topology checking function was used to check

Tolerance = 0.1 * 180 / (6356752.31414 * 3.14)Copy the code
st_overlaps (
	( SELECT ST_Simplify ( multipolygon, tolerance  ) ),
	(
	SELECT
		st_snap (
			ST_Simplify ( multipolygon1, tolerance ),
			( SELECT ST_Simplify ( multipolygon2,tolerance  ) ),
			tolerance 
		) 
	) 
)
Copy the code