PostGIS: Correct way of using ST_Union/ST_Difference with or without gridsize

Question

I'm designing a system where there is a set of polygons representing different postal codes (not overlapping with each other). An user can either draw his own polygon OR fully select/deselect particular postal code. That means that each postal code can be included in the final area fully or partially.

I will use example (brown - user drawn area, light blue - postal code)

I found that when user draws some area, selects a postal code and then deselects it, the output geometry isn't the same as if he would deselect the postal code right away (g1 and g2 columns in below example). Mathematically it doesn't make sense. The g1 figure have weird spikes on the edges of the deselected postal code. How to fix it?

I found out that to fix it I can pass grid size (eg. 1e-9) to the st_union / st_difference, but then st_touches doesn't work correctly, because I can't pass grid size to this function. I can pass it to st_intersection but it is less efficient.

You can see it there:

WITH source AS (
    SELECT st_geomfromgeojson('{"type":"MultiPolygon","coordinates":[[[[-85.92333,40.63875],[-85.92333,40.63392],[-85.91539,40.63492],[-85.91467,40.63875],[-85.92333,40.63875]]]]}') area,
           st_geomfromgeojson('{"type":"MultiPolygon","coordinates":[[[[-85.82047, 40.65114],[-85.82046, 40.6297],[-85.92046, 40.6297],[-85.92068, 40.63807],[-85.92047, 40.65114],[-85.82047, 40.65114]]]]}') postalcode
),
difference AS (
SELECT st_difference(st_union(postalcode, area), postalcode)             g1,
       st_difference(area, postalcode)                                   g2,
       st_difference(st_union(postalcode, area, 1e-9), postalcode, 1e-9) g3,
       area,
       postalcode
FROM source
)
SELECT *, st_touches(g1, postalcode), st_touches(g2, postalcode), st_touches(g3, postalcode)
FROM difference

I understand that this can be because of how float is represented in binary and when the PostGIS calculates points of intersection it can be some irrational number which can't be saved in binary form exactly as it should. Also I found that both of those points from g1 and g2 are the same (POINT (-85.92066907421577 40.63875)).

What should I do?

EDIT: Also forgot to mention, would switching to different srid (some without floating operations) help me?

Answer 1

Using a grid size argument in the overlay functions is a good way to avoid the spikes that you are seeing in the full-precision results.

ST_Touches is a very specific spatial situation. It would be nice if it supported a distance tolerance - perhaps that will be available in some future PostGIS version.

If you want to find postal code polygons which might touch or intersect a polygon, you can use ST_DWithin with a small distance tolerance.

EDIT: To emulate a "touches with tolerance" relationship, there are a few options:

• test whether ST_Intersection(area, postalcode, tol) is linear (as the OP suggested)
• test whether ST_Intersection(area, postalcode) has a small area (based on a fraction of the areas of the input geometries)
• the exact equivalent of "touches with tolerance" is to test whether the boundary of area which lies within postalcode lies within the tolerance distance of the boundary of postalcode. This can be computed as:

ST_Covers( ST_Buffer(ST_Boundary(postalcode), tol),
           ST_Intersection(ST_Boundary(area), postalcode) )

(Note: it would be computationally more efficient to express this as:

ST_OrientedHausdorffDistance( 
     ST_Intersection(ST_Boundary(area), postalcode), 
     ST_Boundary(postalcode) ) 
  <= tol

but unfortunately PostGIS does not (yet) have a function to compute the oriented Hausdorff distance.)

Answer 2

I found some algorithms to remove the spikes that can be seen in the question, and for provided example it worked, but for more complex examples still st_touches is not working. I will leave the links, maybe someone could leverage on that.

https://gasparesganga.com/labs/postgis-normalize-geometry/

https://trac.osgeo.org/postgis/wiki/UsersWikiExamplesSpikeRemover