quicker filling of preliminary region polygons

geovoronoi/_voronoi.py

@@ -291,7 +291,9 @@ def region_polygons_from_voronoi(vor, geom, return_point_assignments=False,
 
  # --- phase 1: preliminary polygons ---
 
- covered_area = 0 # keeps track of area so far covered by generated Voronoi polygons
+ covered_area = 0 # keeps track of area so far covered by generated Voronoi polygons
+ border_regions = set() # keeps track of regions that don't have fully finite bounds, i.e. that are at the edge
+ inner_regions = set() # keeps track of regions that have fully finite bounds, i.e. "inner" regions
  # iterate through regions; `i_reg` is the region ID, `reg_vert` contains vertex indices of this region into
  # `ridge_vert`
  for i_reg, reg_vert in enumerate(vor.regions):
@@ -304,10 +306,13 @@ def region_polygons_from_voronoi(vor, geom, return_point_assignments=False,
 
  # construct the theoretical Polygon `p` of this region
  if np.all(np.array(reg_vert) >= 0): # fully finite-bounded region
+ inner_regions.add(i_reg)
  p = Polygon(vor.vertices[reg_vert])
  else:
  # this region has a least one infinite bound aka "loose ridge"; we go through each ridge and we will need to
  # calculate a finite end ("far point") for loose ridges
+ border_regions.add(i_reg)
+
  p_vert_indices = set() # collect unique indices of vertices into `ridge_vert`
  p_vert_farpoints = set() # collect unique calculated finite far points
 
@@ -396,37 +401,50 @@ def region_polygons_from_voronoi(vor, geom, return_point_assignments=False,
  covered_area += p.area
  region_polys[i_reg] = p
 
- # --- phase 2: filling preliminary polygons ---
+ # --- phase 2: filling preliminary polygons at the edges ---
 
- logger.debug('filling preliminary polygons to fully cover the surrounding area')
+ logger.debug('filling %d preliminary polygons to fully cover the surrounding area' % len(border_regions))
 
  # plot_voronoi_polys_with_points_in_area(ax, geom, region_polys, vor.points, region_pts)
  # fig.savefig(plotfile_fmt % 'beforefill')
 
  uncovered_area_portion = (geom.area - covered_area) / geom.area # initial portion of uncovered area
- polys_iter = iter(region_polys.items())
+ polys_iter = iter(border_regions)
  # the loop has two passes: in the first pass, only areas are considered that intersect with the preliminary
  # Voronoi region polygon; in the second pass, also areas that don't intersect are considered, i.e. isolated features
  pass_ = 1
  # the loop stops when all area is covered or there are no more preliminary Voronoi region polygons left after both
  # passes
+ inner_regions_poly = None
  while not np.isclose(uncovered_area_portion, 0) and 0 < uncovered_area_portion <= 1:
  try:
- i_reg, p = next(polys_iter)
+ i_reg = next(polys_iter)
  except StopIteration:
  if pass_ == 1: # restart w/ second pass
- polys_iter = iter(region_polys.items())
- i_reg, p = next(polys_iter)
+ polys_iter = iter(border_regions)
+ i_reg = next(polys_iter)
  pass_ += 1
  else: # no more preliminary Voronoi region polygons left after both passes
  break
  logger.debug('pass %d / region %d: will fill up %f%% uncovered area'
  % (pass_, i_reg, uncovered_area_portion * 100))
 
- # generate polygon from all other regions' polygons other then the current region `i_reg`
- union_other_regions = cascaded_union([other_poly
- for i_other, other_poly in region_polys.items()
- if i_reg != i_other])
+ p = region_polys[i_reg]
+
+ # generate polygon from all inner regions *once* because this stays constant and is re-used in the loop
+ if not inner_regions_poly and inner_regions:
+ if len(inner_regions) == 1:
+ inner_regions_poly = region_polys[next(iter(inner_regions))]
+ else:
+ inner_regions_poly = cascaded_union([region_polys[i_reg] for i_reg in inner_regions])
+
+ # generate polygon from all other regions' polygons other than the current region `i_reg`
+ other_regions_polys = [region_polys[i_other] for i_other in border_regions if i_reg != i_other]
+ if inner_regions_poly:
+ other_regions_polys.append(inner_regions_poly)
+
+ union_other_regions = cascaded_union(other_regions_polys)
+
  # generate difference between geom and other regions' polygons -- what's left is the current region's area
  # plus any so far uncovered area
  try:

geovoronoi/plotting.py

@@ -14,13 +14,16 @@
 from ._voronoi import points_to_coords, points_to_region
 
 
-def subplot_for_map(show_x_axis=False, show_y_axis=False, aspect='equal', **kwargs):
+def subplot_for_map(show_x_axis=False, show_y_axis=False, show_spines=None, aspect='equal', **kwargs):
  """
  Helper function to generate a matplotlib subplot Axes object suitable for plotting geographic data, i.e. axis
  labels are not shown and aspect ratio is set to 'equal' by default.
 
  :param show_x_axis: show x axis labels
  :param show_y_axis: show y axis labels
+ :param show_spines: controls display of frame around plot; if set to None, this is "auto" mode, meaning
+ that the frame is removed when `show_x_axis` and `show_y_axis` are both set to False;
+ if set to True/False, the frame is always shown/removed
  :param aspect: aspect ratio
  :param kwargs: additional parameters passed to `plt.subplots()`
  :return: tuple with (matplotlib Figure, matplotlib Axes)
@@ -31,6 +34,12 @@ def subplot_for_map(show_x_axis=False, show_y_axis=False, aspect='equal', **kwar
  ax.get_xaxis().set_visible(show_x_axis)
  ax.get_yaxis().set_visible(show_y_axis)
 
+ if show_spines is None:
+ show_spines = show_x_axis or show_y_axis
+
+ for sp in ax.spines.values():
+ sp.set_visible(show_spines)
+
  if show_x_axis:
  fig.autofmt_xdate()
 
@@ -155,11 +164,11 @@ def plot_points(ax, points, markersize=1, marker='o', color=None, labels=None, l
  :param label_color: point labels color
  :param label_draw_duplicates: if False, suppress drawing labels on duplicate coordinates
  :param kwargs: additional parameters passed to matplotlib's `scatter` function
- :return: None
+ :return: return value from `ax.scatter()`
  """
  x, y = xy_from_points(points)
 
- ax.scatter(x, y, s=markersize, marker=marker, color=color, **kwargs)
+ scat = ax.scatter(x, y, s=markersize, marker=marker, color=color, **kwargs)
 
  if labels:
  # plot labels using matplotlib's text()
@@ -176,6 +185,8 @@ def plot_points(ax, points, markersize=1, marker='o', color=None, labels=None, l
  ax.text(x_i, y_i, labels[i], fontsize=label_fontsize, color=_color_for_labels(label_color, color, i))
  drawn_coords.add(pos)
 
+ return scat
+
 
 def plot_line(ax, points, linewidth=1, color=None, **kwargs):
  """
@@ -187,11 +198,11 @@ def plot_line(ax, points, linewidth=1, color=None, **kwargs):
  :param linewidth: line width
  :param color: line color
  :param kwargs: additional parameters passed to matplotlib's `plot` function
- :return: None
+ :return: return value from `ax.plot()`
  """
  x, y = xy_from_points(points)
 
- ax.plot(x, y, linewidth=linewidth, color=color, **kwargs)
+ return ax.plot(x, y, linewidth=linewidth, color=color, **kwargs)
 
 
 def plot_polygon(ax, polygon, facecolor=None, edgecolor=None, linewidth=1, linestyle='solid',
@@ -209,15 +220,17 @@ def plot_polygon(ax, polygon, facecolor=None, edgecolor=None, linewidth=1, lines
  :param label_fontsize: label font size
  :param label_color: label color
  :param kwargs: additonal parameters passed to matplotlib `PatchCollection` object
- :return: None
+ :return: return value from `ax.add_collection()`
  """
- ax.add_collection(PatchCollection([PolygonPatch(polygon)],
- facecolor=facecolor, edgecolor=edgecolor,
- linewidth=linewidth, linestyle=linestyle,
- **kwargs))
+ coll = ax.add_collection(PatchCollection([PolygonPatch(polygon)],
+  facecolor=facecolor, edgecolor=edgecolor,
+  linewidth=linewidth, linestyle=linestyle,
+  **kwargs))
  if label:
  ax.text(polygon.centroid.x, polygon.centroid.y, label, fontsize=label_fontsize, color=label_color)
 
+ return coll
+
 
 def plot_voronoi_polys_with_points_in_area(ax, area_shape, region_polys, points, region_pts=None,
  area_color=(1,1,1,1), area_edgecolor=(0,0,0,1),

tests/test_main.py

@@ -216,7 +216,7 @@ def test_voronoi_italy_with_plot(n_pts, per_geom):
  assert 0 < len(region_polys) <= n_pts
 
  # generate plot
- fig, ax = subplot_for_map()
+ fig, ax = subplot_for_map(show_spines=True)
  plot_voronoi_polys_with_points_in_area(ax, area_shape, region_polys, coords, region_pts,
  point_labels=list(map(str, range(len(coords)))))
 
@@ -279,7 +279,7 @@ def test_voronoi_geopandas_with_plot():
  assert len(region_polys) == len(region_pts) == len(coords)
 
  # generate plot
- fig, ax = subplot_for_map()
+ fig, ax = subplot_for_map(show_spines=True)
  plot_voronoi_polys_with_points_in_area(ax, south_am_shape, region_polys, coords, region_pts)
 
  return fig
@@ -317,7 +317,7 @@ def test_voronoi_sweden_duplicate_points_with_plot():
  for i_poly, pt_indices in region_pts.items()}
 
  # generate plot
- fig, ax = subplot_for_map()
+ fig, ax = subplot_for_map(show_spines=True)
  plot_voronoi_polys_with_points_in_area(ax, area_shape, region_polys, distinct_pt_coords,
  plot_voronoi_opts={'alpha': 0.2},
  plot_points_opts={'alpha': 0.4},
@@ -363,7 +363,7 @@ def test_issue_7b():
 
  assert all([len(pts_in_region) == 1 for pts_in_region in region_pts.values()]) # no duplicates
 
- fig, ax = subplot_for_map()
+ fig, ax = subplot_for_map(show_spines=True)
  plot_voronoi_polys_with_points_in_area(ax, polygon, region_polys, centroids, region_pts)
 
  return fig