Opened 9 years ago
Last modified 2 years ago
#3447 new defect
PG::InternalError: ERROR: Error computing Voronoi diagram.
Reported by: | loarie | Owned by: | dbaston |
---|---|---|---|
Priority: | medium | Milestone: | PostGIS GEOS |
Component: | postgis | Version: | master |
Keywords: | ST_Voronoi | Cc: |
Description
I'm running ST_Voronoi on lots of different multipoints and its very odd that sometimes very complex multipoints run successfully with very small tolerance settings while other very simple multipoints fail unless unusably high tolerances are set.
Seems like there must be a bug here. Happy to help diagnose/provide failing/passing multipoint examples if someone will help fix it
Attachments (1)
Change History (13)
comment:2 by , 9 years ago
Owner: | changed from | to
---|
by , 9 years ago
Attachment: | voronoi_example added |
---|
the voronoi_test_pt table used in the failed ST_Voronoi example in my comment
comment:3 by , 9 years ago
Great thanks dbaston - sorry for the delay, but here's a really simple table where ST_voronoi doesn't work
bounds2=# SELECT postgis_full_version();;
postgis_full_version
POSTGIS="2.3.0dev r14623" GEOS="3.5.0-CAPI-1.9.0 r4084" PROJ="Rel. 4.9.2, 08 September 2015" GDAL="GDAL 1.11.3, released 2015/09/16" LIBXML="2.9.0" LIBJSON="0.12" TOPOLOGY RASTER
(1 row)
bounds2=# SELECT name, ST_astext(ST_Voronoi(geom, null,0.002)) as geom FROM voronoi_test_pt WHERE name = 'all'; ERROR: Error computing Voronoi diagram.
COPY voronoi_test_pt TO '/voronoi_example'
comment:4 by , 9 years ago
Thanks for the example loarie. Here are my observations so far:
- The problem occurs within GEOS.
- JTS has no problem producing a Voronoi diagram for this input.
- ST_DelaunayTriangles has no problem with this input; I'd have expected both or none to work
I'll keep looking around to see if I can turn up anything obvious.
GEOS error is IllegalArgumentException: Points of LinearRing do not form a closed linestring
, from
#0 0x00007f93f2c498b0 in __cxa_throw () from /usr/lib/x86_64-linux-gnu/libstdc++.so.6 #1 0x00007f93f2fae617 in geos::geom::LinearRing::validateConstruction (this=this@entry=0x7f94068cbc00) at LinearRing.cpp:70 #2 0x00007f93f2fae71b in geos::geom::LinearRing::LinearRing (this=0x7f94068cbc00, newCoords=0x7f94068cbbe0, newFactory=<optimized out>, __in_chrg=<optimized out>, __vtt_parm=<optimized out>) at LinearRing.cpp:46 #3 0x00007f93f2fac115 in geos::geom::GeometryFactory::createLinearRing (this=this@entry=0x7f94067fa4a0, newCoords=newCoords@entry=0x7f94068cbbe0) at GeometryFactory.cpp:469 #4 0x00007f93f303531b in geos::triangulate::quadedge::QuadEdgeSubdivision::getVoronoiCellPolygon (this=this@entry=0x7f940689efa0, qe=<optimized out>, geomFact=...) at QuadEdgeSubdivision.cpp:582 #5 0x00007f93f3036c0e in geos::triangulate::quadedge::QuadEdgeSubdivision::getVoronoiCellPolygons (this=0x7f940689efa0, geomFact=...) at QuadEdgeSubdivision.cpp:528 #6 0x00007f93f3036cc9 in geos::triangulate::quadedge::QuadEdgeSubdivision::getVoronoiDiagram (this=<optimized out>, geomFact=...) at QuadEdgeSubdivision.cpp:505 #7 0x00007f93f303265b in geos::triangulate::VoronoiDiagramBuilder::getDiagram (this=this@entry=0x7ffd18db06d0, geomFact=...) at VoronoiDiagramBuilder.cpp:106 #8 0x00007f93f414042e in GEOSVoronoiDiagram_r (extHandle=0x7f94067fa050, g1=0x7f94068bb7e0, env=0x0, tolerance=0, onlyEdges=0) at geos_ts_c.cpp:6519 #9 0x00007f93f43c7fab in lwgeom_voronoi_diagram (g=g@entry=0x7f940683a670, env=env@entry=0x0, tolerance=tolerance@entry=0, output_edges=0) at lwgeom_geos.c:2013 #10 0x00007f93f4374977 in ST_Voronoi (fcinfo=0x7f94068159e0) at lwgeom_geos.c:3611 #11 0x00007f94042882f3 in ?? () #12 0x00007f940428debc in ExecProject () #13 0x00007f9404293f90 in ExecAgg () #14 0x00007f9404287068 in ExecProcNode () #15 0x00007f94042845e2 in standard_ExecutorRun () #16 0x00007f9404367edf in ?? () #17 0x00007f94043694aa in PortalRun () #18 0x00007f9404366e65 in PostgresMain () #19 0x00007f940415a280 in ?? () #20 0x00007f9404321ab1 in PostmasterMain () #21 0x00007f940415af03 in main ()
comment:5 by , 9 years ago
I take that back - JTS produced a result for this case only because of the rounding that occurred when I transferred the geometry using WKT. In fact, the problem can be produced in both PostGIS and JTS with a smaller example of only 7 points:
SELECT ST_Voronoi('01040000000700000001010000000f8b33e3d97742c038c453588d0423c001010000001171d6d1b45d42c06adc1693e78c22c001010000001c8b33e3d97742c062c453588d0423c00101000000afa5c71fda7742c04b93c61d8e0423c00101000000b0cddcb4b57942c026476887d7b122c00101000000e0678421dc7642c0f7736021e1fb22c00101000000e32fd565018d42c0c7ea1222167c22c0', tolerance := 0.02);
This should be probably reported upstream to JTS; I don't think there's much that can be done at the PostGIS level.
For your particular case, it seems like the problem can be solved by running ST_SnapToGrid on your inputs before sending them to ST_Voronoi. I was able to get a diagram on your inputs with a tolerance of 0.0 if I ran ST_SnapToGrid(geom, 1e-18) beforehand:
SELECT name, ST_AsText(ST_Voronoi(ST_SnapToGrid(geom, 1e-18))) as geom FROM voronoi_test_pt WHERE name = 'all';
comment:7 by , 9 years ago
Milestone: | PostGIS 2.3.0 → PostGIS GEOS |
---|
dbaston so sounds like this is a bug in GEOS (to be fixed after fixed in GEOS). Switching to a GEOS milestone.
comment:9 by , 9 years ago
Also ticket on JTS as strk noted in IRC - https://github.com/dr-jts/jts/ since it seems to stemp from JTS logic.
comment:10 by , 9 years ago
OK, I created tickets with GEOS (https://trac.osgeo.org/geos/ticket/769) and JTS (https://github.com/dr-jts/jts/issues/20).
robe2, should this ticket be closed? There is no action required in PostGIS after the issue is fixed in GEOS.
comment:11 by , 9 years ago
No. We don't close the ticket until the issue is fixed upstream, then we close out with a note about the versions that have the fix. We just move the milestone to PostGIS GEOS or PostGIS GDAL etc, to so it is not confused with actual fixes we do in our own code.
Yes, please post any examples that don't behave as you'd expect, as well as details of your platform using SELECT postgis_full_version()