I have a workflow where I am creating Voronoi polygons from a series of points. However since the polygons on the exterior extend really far out I have clipped them and then joined them back to points to get the IDs back. However when I join them back based on a ST_Intersects about 1-2% of them get lost. Here is the original Voronoi polygons and then after the intersection.

Voronoi before cliping

Voronoi after clipping and ST_intersects

Here is the code that I am using

WITH voronoi_temp2 AS (
    SELECT (ST_dump(ST_VoronoiPolygons(ST_Collect(the_geom)))).geom as the_geom
    FROM ways_vertices_pgr)
SELECT id, voronoi_temp2.the_geom
FROM voronoi_temp2  
JOIN ways_vertices_pgr
ON ST_Intersects(ways_vertices_pgr.the_geom, voronoi_temp2.the_geom);

Upon further inspection of the geometries those that drop out have duplicated points. So something like [(0, 0), (1, 0.5), (1, 0.5), (0, 1), (0, 0)]. So they are self-intersecting and invalid. I have tried MakeValid and RemoveRepeatedPoints but both end up dropping the invalid polygons.

This only occurs when the data is in a geographic coordinate system and not when in a projected coordinate system.

  • Do you get the missing triangles if you use ST_Disjoint instead of ST_Intersects? – user30184 Dec 11 '17 at 13:28
  • I do get the triangle with disjoint, but I am not sure why. If I create voronoi polygons from points they 100% should intersect with the points. – James Dec 11 '17 at 14:02
  • I agree that they should. Disjoint is opposite to Intersects so at least PostGIS seems to behave in a logical way. I would try to select the geometries of the disjoint triangles as WKT and check if geometries are closed. – user30184 Dec 11 '17 at 15:15
  • 1
    Upon further inspection I found that they are self-intersecting. I have tried MakeValid and buffering them by 0 to fix them but both of those drop some of the rows. – James Dec 11 '17 at 16:02
  • 2
    If you have a repeatable test case that generates invalid geometry, you should open a ticket (or at least, add that info to this page). I cannot repeat your result from synthetic data. with polys as (select (st_dump(ST_VoronoiPolygons(st_collect(st_makepoint(s/0.2, t/0.2))))).geom as geom from generate_series(1,50) s, generate_series(1,50) t) select st_astext(geom), st_isvalid(geom) from polys; – Paul Ramsey Dec 12 '17 at 20:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.