[postgis-tickets] [PostGIS] #4711: ST_Union loses precision on complex multilinestring geometries

[PostGIS]

#4711: ST_Union loses precision on complex multilinestring geometries
-------------------------+---------------------------
  Reporter:  dannytoone  |      Owner:  pramsey
      Type:  defect      |     Status:  new
  Priority:  medium      |  Milestone:  PostGIS 3.1.0
 Component:  postgis     |    Version:  2.5.x
Resolution:              |   Keywords:
-------------------------+---------------------------

Comment (by dannytoone):

 I believe I have a reproducible example using generated data. Please let
 me know if you have any problems with running or generating the dataset.

 Unfortunately, the geometries are not simple, and simple geometries did
 not reproduce this bug, so creating a dataset that recreates the bug is
 not simple.

 For background, the data I'm working on, that I'm trying to emulate with
 this generated data set, is roughly analogous to the license geometries in
 the MDS/ITFS frequency bands (2.5-2.7ghz), as regulated by the FCC. They
 do not overlap when viewing the same frequencies, but across frequencies
 they may overlap. The following temp tables set up a data structure that
 is somewhat similar to these license boundaries.

 {{{
 DROP TABLE IF EXISTS areas;
 CREATE TEMP TABLE areas AS
 WITH centerpoints AS (
     SELECT
         row_number() over () as id,
         floor(random() * 10 + 1) as lic_group,
         ST_SetSRID(ST_MakePoint(random() * 10 + 1, random() * 10 +
 1),4326) as centerpoint
     FROM generate_series(1,10000)
 )
 SELECT
     id,
     lic_group,
     centerpoint,
     ST_Buffer(centerpoint::geography,10000)::geometry as geom
 FROM centerpoints
 ;
 CREATE INDEX idx_geom ON areas USING GIST (geom);
 CREATE INDEX idx_group ON areas (lic_group);

 DROP TABLE IF EXISTS area_exclusions;
 CREATE TEMP TABLE area_exclusions AS
 WITH split_lines AS (
     SELECT
         a.id,
         b.id as b_id,
         a.lic_group,
         a.centerpoint,
         a.geom as a_geom,
         b.geom as b_geom,
 ST_LineFromMultiPoint(ST_Intersection(ST_ExteriorRing(a.geom),ST_ExteriorRing(b.geom)))
 as split_line
     FROM areas a
     INNER JOIN areas b
         ON ST_Intersects(a.geom,b.geom)
         AND a.lic_group = b.lic_group
         AND a.id != b.id
 ), split_areas AS (
     SELECT
         id,
         b_id,
         lic_group,
         centerpoint,
         ST_Split(ST_Union(a_geom,b_geom),split_line) as split_geom
     FROM split_lines
 )
 SELECT
     id,
     b_id,
     lic_group,
     CASE ST_ContainsProperly(ST_GeometryN(split_geom,1),centerpoint)
         WHEN false THEN ST_GeometryN(split_geom,1)
         WHEN true THEN ST_GeometryN(split_geom,2)
     END as geom
 FROM split_areas
 ;
 CREATE INDEX idx_e_id ON area_exclusions (id,lic_group);

 DROP TABLE IF EXISTS processed_areas;
 CREATE TEMP TABLE processed_areas AS
 WITH exclusions AS (
     SELECT
         id,
         lic_group,
         ST_Union(geom) as geom
     FROM area_exclusions
     GROUP BY
         id,
         lic_group
 )
 SELECT
     id,
     lic_group,
     ST_Difference(i.geom,COALESCE(e.geom,ST_GeomFromEWKT('SRID=4326;
 POLYGON EMPTY'))) as geom
 FROM areas i
 LEFT JOIN exclusions e
     USING (id,lic_group)
 ;
 }}}

 If I am to query a single `lic_group` I can see that they do not have
 overlapping geometries:

 {{{
 SELECT geom
 FROM processed_areas
 WHERE lic_group = 1;
 }}}

 However, across several `lic_group`s, there is plenty of overlap:

 {{{
 SELECT geom
 FROM processed_areas;
 }}}
 .

 =======================
 Okay, that is the underlying data. Now onto the actual bug. This table is
 merely the processed geometries as an exterior ring linestring.


 {{{
 DROP TABLE IF EXISTS rings;
 CREATE TEMP TABLE rings AS
 SELECT id, lic_group, ST_ExteriorRing(geom) as geom
 FROM processed_areas
 ;
 }}}

 Viewing the rings themselves, we can see there is plenty of precision in
 the underlying points contributing to the linestrings:

 {{{
 SELECT Left(ST_AsText(geom),200)
 FROM rings;
 =====================
 LINESTRING(10.8643881174031 3.77734188849076,10.8626209909011
 3.75970533281878,10.8574629500733 3.74275385983207,10.8491123276953
 3.72713876384977,10.8378900600966 3.71345998467201,10.8242273509016 3.
 LINESTRING(7.57876514569725 2.81432363450241,7.57705816723587
 2.7966769867886,7.5719601905669 2.77970431048747,7.56366729736813
 2.76405788049345,7.55249832730074 2.75033896339531,7.53888260377069 2.73
 LINESTRING(8.47788820510008 6.81387172768232,8.46467044867379
 6.8125559523933,8.44701315022753 6.81427415953191,8.43003023487508
 6.81940188537274,8.41437401669197 6.82774214951577,8.40064591887797 6.8
 ...
 }}}

 However, once unioned, that precision is lost:


 {{{
 SELECT Left(ST_AsText(ST_Union(geom)),200)
 FROM rings;
 ==================
 MULTILINESTRING((8 3,7 3),(8 7,9 7),(2 10,2 9),(10 4,9 4),(4 8,3 8),(2 1,2
 2),(5 10,5 11),(9 11,8 11),(9 9,8 9),(3 9,2 9),(11 1,10 1),(10 1,11 2),(11
 2,11 1),(3 5,3 4),(5 2,5 1),(10 6,9 6),(9 10,8 10)...
 }}}

 Again, please let me know if this example works for you to be able to
 reproduce.

-- 
Ticket URL: <https://trac.osgeo.org/postgis/ticket/4711#comment:2>
PostGIS <http://trac.osgeo.org/postgis/>
The PostGIS Trac is used for bug, enhancement & task tracking, a user and developer wiki, and a view into the subversion code repository of PostGIS project.