fix problem with wrong intersection points from outside finite points

examples/random_points_across_italy.py

@@ -25,7 +25,7 @@
 
 #%%
 
-N_POINTS = 100
+N_POINTS = 105
 COUNTRY = 'Italy'
 
 np.random.seed(123)
@@ -56,14 +56,24 @@
 #%%
 
 #
-# calculate the Voronoi regions, cut them with the geographic area shape and assign the points to them
+# Calculate the Voronoi regions, cut them with the geographic area shape and assign the points to them
+# Note how in Sardinia there's only one point which is not assigned to any Voronoi region.
+# Since by default all sub-geometries (i.e. islands or other isolated features) in the geographic shape
+# are treated separately (set `per_geom=False` to change this) and you need more than one point to generate
+# a Voronoi region, this point is left unassigned. By setting `return_unassigned_points=True`, we can get a
+# set of unassigned point indices:
 #
 
-poly_shapes, poly_to_pt_assignments = voronoi_regions_from_coords(pts, area_shape)
+poly_shapes, poly_to_pt_assignments, unassigned_pts = voronoi_regions_from_coords(pts, area_shape,
+                                                                                  return_unassigned_points=True)
 
 print('Voronoi region to point assignments:')
 pprint(poly_to_pt_assignments)
 
+print('Unassigned points:')
+for i_pt in unassigned_pts:
+    print('#%d: %.2f, %.2f' % (i_pt, pts[i_pt].x, pts[i_pt].y))
+
 #%%
 
 #

examples/toyexample.py

@@ -28,6 +28,7 @@
 #                    [2, 0], [2, 1], [2, 2]])
 
 points = np.array([[1, 1], [1.1, 0.9], [1.15, 0.8], [2.5, 0]])
+#points = np.array([[1, 1], [1.1, 0.91], [1.2, 0.8]])
 
 # surrounding shape
 shape = Polygon([[-1, -1], [3, -1], [3, 3], [-1, 3]])

geovoronoi/_geom.py

@@ -1,8 +1,6 @@
 """
 Geometry helper functions in cartesian 2D space.
 
-"shapely" refers to the [Shapely Python package for computational geometry](http://toblerity.org/shapely/index.html).
-
 Author: Markus Konrad <markus.konrad@wzb.eu>
 """
 

geovoronoi/_voronoi.py

@@ -32,7 +32,7 @@ def points_to_coords(pts):
     return np.array([p.coords[0] for p in pts])
 
 
-def voronoi_regions_from_coords(coords, geo_shape, per_geom=True):
+def voronoi_regions_from_coords(coords, geo_shape, per_geom=True, return_unassigned_points=False):
     """
     Calculate Voronoi regions from NumPy array of 2D coordinates `coord` that lie within a shape `geo_shape`. Setting
     `shapes_from_diff_with_min_area` fixes rare errors where the Voronoi shapes do not fully cover `geo_shape`. Set this
@@ -76,6 +76,7 @@ def voronoi_regions_from_coords(coords, geo_shape, per_geom=True):
     pts_indices = set(range(len(pts)))
     region_polys = {}
     region_pts = {}
+    unassigned_pts = set()
     for i_geom, geom in enumerate(geoms):
         pts_in_geom = [i for i in pts_indices if geom.contains(pts[i])]
         logger.info('%d of %d points in geometry #%d of %d' % (len(pts_in_geom), len(pts), i_geom + 1, len(geoms)))
@@ -91,6 +92,7 @@ def voronoi_regions_from_coords(coords, geo_shape, per_geom=True):
             if exc.args and 'QH6214' in exc.args[0]:
                 logger.error('not enough input points (%d) for Voronoi generation; original error message: %s' %
                              (len(pts_in_geom), exc.args[0]))
+                unassigned_pts.update(pts_in_geom)
                 continue
             raise exc
 
@@ -107,7 +109,10 @@ def voronoi_regions_from_coords(coords, geo_shape, per_geom=True):
                      for reg_id, old_reg_id in region_ids_mapping.items()}
         )
 
-    return region_polys, region_pts
+    if return_unassigned_points:
+        return region_polys, region_pts, unassigned_pts
+    else:
+        return region_polys, region_pts
 
 
 def region_polygons_from_voronoi(vor, geom, return_point_assignments=False):
@@ -150,38 +155,39 @@ def region_polygons_from_voronoi(vor, geom, return_point_assignments=False):
                         finite_pt = vor.vertices[i]
                         p_vert_indices.add(i)
 
-                        # only add a far point if the finite end is not outside the bounding box already
-                        if geom_bb.contains(asPoint(finite_pt)):
-                            t = vor.points[pointidx[1]] - vor.points[pointidx[0]]  # tangent
-                            t /= np.linalg.norm(t)
-                            n = np.array([-t[1], t[0]])  # normal
-
-                            midpoint = vor.points[pointidx].mean(axis=0)
-                            direction = np.sign(np.dot(midpoint - center, n)) * n
-
-                            isects = []
-                            for i_ext_coord in range(len(geom_bb.exterior.coords) - 1):
-                                isect = line_segment_intersection(midpoint, direction,
-                                                                  np.array(geom_bb.exterior.coords[i_ext_coord]),
-                                                                  np.array(geom_bb.exterior.coords[i_ext_coord+1]))
-                                if isect is not None:
-                                    isects.append(isect)
-
-                            if len(isects) == 0:
-                                raise RuntimeError('far point must intersect with surrounding geometry from `geom`')
-                            elif len(isects) == 1:
-                                p_vert_farpoints.add(tuple(isects[0]))
-                            else:
-                                closest_isect_idx = np.argmin(np.linalg.norm(midpoint - isects, axis=1))
-                                p_vert_farpoints.add(tuple(isects[closest_isect_idx]))
+                        t = vor.points[pointidx[1]] - vor.points[pointidx[0]]  # tangent
+                        t /= np.linalg.norm(t)
+                        n = np.array([-t[1], t[0]])  # normal
+
+                        midpoint = vor.points[pointidx].mean(axis=0)
+                        direction = np.sign(np.dot(midpoint - center, n)) * n
+
+                        isects = []
+                        for i_ext_coord in range(len(geom_bb.exterior.coords) - 1):
+                            isect = line_segment_intersection(midpoint, direction,
+                                                              np.array(geom_bb.exterior.coords[i_ext_coord]),
+                                                              np.array(geom_bb.exterior.coords[i_ext_coord+1]))
+                            if isect is not None:
+                                isects.append(isect)
+
+                        if len(isects) == 0:
+                            raise RuntimeError('far point must intersect with surrounding geometry from `geom`')
+                        elif len(isects) == 1:
+                            far_pt = isects[0]
+                        else:
+                            closest_isect_idx = np.argmin(np.linalg.norm(midpoint - isects, axis=1))
+                            far_pt = isects[closest_isect_idx]
+
+                        if (far_pt - finite_pt)[0] / direction[0] > 0:   # only if in same direction
+                            p_vert_farpoints.add(tuple(far_pt))
 
             # create the Voronoi region polygon as convex hull of the ridge vertices and far points (Voronoi regions
             # are convex)
 
             p_pts = vor.vertices[np.asarray(list(p_vert_indices))]
             # adding the point itself prevents generating invalid hull (otherwise sometimes the hull becomes a line
             # if the vertices are almost colinear)
-            p_pts = np.vstack((p_pts, vor.points[i_pt]))
+            p_pts = np.vstack((p_pts, vor.points[pt_indices]))
             if p_vert_farpoints:
                 p_pts = np.vstack((p_pts, np.asarray(list(p_vert_farpoints))))