fixed possible div by zero, fixed problem summing up area in region f…

…illing phase 2
WZBSocialScienceCenter · Feb 2, 2021 · 1fd852e · 1fd852e
1 parent 0ea1889
commit 1fd852e
Show file tree

Hide file tree

Showing 21 changed files with 245 additions and 213 deletions.
diff --git a/geovoronoi/_geom.py b/geovoronoi/_geom.py
@@ -75,5 +75,8 @@ def calculate_polygon_areas(poly_shapes, m2_to_km2=False):
     False) or km² (otherwise).
     """
 
+    if not isinstance(poly_shapes, dict):
+        raise ValueError('`poly_shapes` must be a dict')
+
     unit_convert = 1000000 if m2_to_km2 else 1
     return {i_poly: p.area / unit_convert for i_poly, p in poly_shapes.items()}
diff --git a/geovoronoi/_voronoi.py b/geovoronoi/_voronoi.py
@@ -13,6 +13,7 @@
 from scipy.spatial import Voronoi
 from scipy.spatial.qhull import QhullError
 from shapely.geometry import box, LineString, asPoint, MultiPoint, Polygon, MultiPolygon
+from shapely.errors import TopologicalError
 from shapely.ops import cascaded_union
 
 from ._geom import line_segment_intersection
@@ -207,7 +208,8 @@ def region_polygons_from_voronoi(vor, geom, return_point_assignments=False):
                             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
+                        nonzero_coord = np.nonzero(direction)[0][0]
+                        if (far_pt - finite_pt)[nonzero_coord] / direction[nonzero_coord] > 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
@@ -241,22 +243,25 @@ def region_polygons_from_voronoi(vor, geom, return_point_assignments=False):
     uncovered_area_portion = (geom.area - covered_area) / geom.area
     polys_iter = iter(region_polys.items())
     pass_ = 1
-    while not np.isclose(uncovered_area_portion, 0) and uncovered_area_portion > 0:
+    while not np.isclose(uncovered_area_portion, 0) and 0 < uncovered_area_portion <= 1:
         try:
             i_reg, p = next(polys_iter)
         except StopIteration:
             if pass_ == 1:   # restart w/ second pass
                 polys_iter = iter(region_polys.items())
                 i_reg, p = next(polys_iter)
-                pass_ += 2
+                pass_ += 1
             else:
                 break
         logger.debug('will fill up %f%% uncovered area' % (uncovered_area_portion * 100))
 
         union_other_regions = cascaded_union([other_poly
                                               for i_other, other_poly in region_polys.items()
                                               if i_reg != i_other])
-        diff = geom.difference(union_other_regions)
+        try:
+            diff = geom.difference(union_other_regions)
+        except TopologicalError:   # may happen in rare circumstances
+            continue
         if isinstance(diff, MultiPolygon):
             diff = diff.geoms
         else:
@@ -276,7 +281,10 @@ def region_polygons_from_voronoi(vor, geom, return_point_assignments=False):
                     center_to_center = LineString([p.centroid, diff_part.centroid])
                     if not any(center_to_center.crosses(region_polys[i_neighb]) for i_neighb in neighbor_regions):
                         add.append(diff_part)
-                        covered_area += (diff_part.area - p.area)
+                        if pass_ == 1:
+                            covered_area += (diff_part.area - p.area)
+                        else:
+                            covered_area = diff_part.area
 
         if add:
             region_polys[i_reg] = cascaded_union([p] + add)

diff --git a/tests/baseline/test_issue_7b.png b/tests/baseline/test_issue_7b.png
diff --git a/tests/baseline/test_voronoi_geopandas_with_plot.png b/tests/baseline/test_voronoi_geopandas_with_plot.png
diff --git a/tests/baseline/test_voronoi_italy_with_plot.png b/tests/baseline/test_voronoi_italy_with_plot.png
diff --git a/tests/baseline/test_voronoi_italy_with_plot_10-False.png b/tests/baseline/test_voronoi_italy_with_plot_10-False.png
diff --git a/tests/baseline/test_voronoi_italy_with_plot_10-True.png b/tests/baseline/test_voronoi_italy_with_plot_10-True.png
diff --git a/tests/baseline/test_voronoi_italy_with_plot_100-False.png b/tests/baseline/test_voronoi_italy_with_plot_100-False.png
diff --git a/tests/baseline/test_voronoi_italy_with_plot_100-True.png b/tests/baseline/test_voronoi_italy_with_plot_100-True.png
diff --git a/tests/baseline/test_voronoi_italy_with_plot_1000-False.png b/tests/baseline/test_voronoi_italy_with_plot_1000-False.png
diff --git a/tests/baseline/test_voronoi_italy_with_plot_1000-True.png b/tests/baseline/test_voronoi_italy_with_plot_1000-True.png
diff --git a/tests/baseline/test_voronoi_italy_with_plot_20-False.png b/tests/baseline/test_voronoi_italy_with_plot_20-False.png
diff --git a/tests/baseline/test_voronoi_italy_with_plot_20-True.png b/tests/baseline/test_voronoi_italy_with_plot_20-True.png
diff --git a/tests/baseline/test_voronoi_italy_with_plot_50-False.png b/tests/baseline/test_voronoi_italy_with_plot_50-False.png
diff --git a/tests/baseline/test_voronoi_italy_with_plot_50-True.png b/tests/baseline/test_voronoi_italy_with_plot_50-True.png
diff --git a/tests/baseline/test_voronoi_italy_with_plot_500-False.png b/tests/baseline/test_voronoi_italy_with_plot_500-False.png
diff --git a/tests/baseline/test_voronoi_italy_with_plot_500-True.png b/tests/baseline/test_voronoi_italy_with_plot_500-True.png
diff --git a/tests/baseline/test_voronoi_spain_area_with_plot.png b/tests/baseline/test_voronoi_spain_area_with_plot.png
diff --git a/tests/baseline/test_voronoi_sweden_duplicate_points_with_plot.png b/tests/baseline/test_voronoi_sweden_duplicate_points_with_plot.png
diff --git a/tests/test_geom.py b/tests/test_geom.py
@@ -38,21 +38,29 @@ def test_line_segment_intersection(l_off, l_dir, segm_a, segm_b, expected):
 
 
 def test_calculate_polygon_areas_empty():
-    areas = calculate_polygon_areas([])
+    areas = calculate_polygon_areas({})
     assert len(areas) == 0
 
 
+def test_calculate_polygon_areas_nondict():
+    with pytest.raises(ValueError):
+        calculate_polygon_areas([])
+
+
 def test_calculate_polygon_areas_world():
     import geopandas as gpd
 
     world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
     world = world[world.continent != 'Antarctica'].to_crs(epsg=3395)  # meters as unit!
+    geoms = {i: geom for i, geom in enumerate(world.geometry)}
 
-    areas = calculate_polygon_areas(world.geometry)
+    areas = calculate_polygon_areas(geoms)
 
+    assert isinstance(areas, dict)
     assert len(areas) == len(world)
-    assert all(0 <= a < 9e13 for a in areas)
+    assert all(0 < a < 9e13 for a in areas.values())
 
-    areas_km2 = calculate_polygon_areas(world.geometry, m2_to_km2=True)
+    areas_km2 = calculate_polygon_areas(geoms, m2_to_km2=True)
+    assert isinstance(areas_km2, dict)
     assert len(areas_km2) == len(world)
-    assert all(np.isclose(a_m, a_km * 1e6) for a_m, a_km in zip(areas, areas_km2))
+    assert all(np.isclose(a_m, a_km * 1e6) for a_m, a_km in zip(areas.values(), areas_km2.values()))