[postgis-tickets] [PostGIS] #5392: Aggregate ST_Union fails to join some "circles" when radius is too large - lwgeom_unaryunion_prec: GEOS Error: TopologyException

[PostGIS]

#5392: Aggregate ST_Union fails to join some "circles" when radius is too large -
lwgeom_unaryunion_prec: GEOS Error: TopologyException
----------------------+---------------------------
  Reporter:  rotten   |      Owner:  pramsey
      Type:  defect   |     Status:  closed
  Priority:  medium   |  Milestone:  PostGIS 3.3.4
 Component:  postgis  |    Version:  3.3.x
Resolution:  invalid  |   Keywords:
----------------------+---------------------------
Changes (by pramsey):

 * status:  new => closed
 * resolution:   => invalid

Comment:

 This is less interesting than it appears to be on first inspection (but it
 is still pretty interesting). Because the SQL starts by taking buffers of
 points, which intuitively should be easy-to-process circles, the result is
 surprising to our eyes.

 However, taking the SQL apart into something a little more step wise, the
 real problem is just the normal boring one of feeding invalid inputs to
 the overlay functions.

 {{{
 create table polys (name text, geom geometry);

 insert into polys values ('kpc', ST_Buffer(
             ST_Point(-166.9324011::float8, 65.154889::float8,
 4326)::geography,
             290086::integer
         )::geometry);

 insert into polys values ('tonga', ST_Buffer(
             ST_Point(-177.2519031::float8, -18.5973127::float8,
 4326)::geography,
             290086::integer
         )::geometry);

 select name, st_isvalid(geom) from polys;
 }}}

 The main thing to note is that taking large buffers of objects in
 geography may yield invalid outputs, since the buffer function in
 geography is a nasty hack. In geometry of course a buffer of a point is a
 circle and never invalid, but that's not what is happening here.

 There is still an "interesting" bit left, and that is is that the pairwise
 and the aggregate versions of union are returning different results.

 {{{
 select st_union(geom) from polys;

 select st_union(a.geom, b.geom)
 from polys a
 cross join polys b
 where a.name = 'kpc' and b.name = 'tonga';
 }}}

 If you backtrace from inside the overlay code where that TopologyException
 is thrown, you'll find that the process is actually going through
 different code lines. The aggregate goes through the cascaded union, which
 is very complex, while the pairwise just goes straight into the overlay
 code.

 If you make valid inputs, you can get a result from both the aggregate and
 pairwise union code lines.

 {{{
 insert into polys
 select name || ' valid', st_makevalid(geom) as geom
 from polys
 where name in ('kpc', 'tonga');

 select st_area(st_union(geom)) from polys where name ~ 'valid';

 select st_area(st_union(a.geom, b.geom))
 from polys a
 cross join polys b
 where a.name = 'kpc valid'
 and b.name = 'tonga valid'
 }}}

 Note however, that if you compare the result of the pairwise union of the
 invalid inputs to the result for the valid inputs, the answer is wrong.

 {{{
 select st_area(st_union(a.geom, b.geom))
 from polys a
 cross join polys b
 where a.name = 'kpc'
 and b.name = 'tonga'
 }}}

 So the moral of the story is not "pairwise union is more resistant to
 invalid inputs than aggregate union", it is "invalid inputs can fail in
 different ways through different code lines, so you really want to avoid
 invalid inputs".

 I sometimes wonder if we should do automagical makevalid as part of the
 try/catch block on the overlay code, so the case of the TopologyException
 would end up actually working, automagically (catch exception, make valid,
 retry). But seeing the silent failure of the pairwise code makes me
 question that. Automagical makevalid would make those cases of silent
 failure even more likely to just slide by people's notice.
-- 
Ticket URL: <https://trac.osgeo.org/postgis/ticket/5392#comment:5>
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.