plotting works

examples/random_points_across_italy.py

@@ -13,7 +13,7 @@
 import numpy as np
 import geopandas as gpd
 
-from geovoronoi import coords_to_points, points_to_coords, voronoi_regions_from_coords
+from geovoronoi import coords_to_points, voronoi_regions_from_coords
 from geovoronoi.plotting import subplot_for_map, plot_voronoi_polys_with_points_in_area
 
 
@@ -22,6 +22,8 @@
 geovoronoi_log.setLevel(logging.INFO)
 geovoronoi_log.propagate = True
 
+#%%
+
 N_POINTS = 100
 COUNTRY = 'Italy'
 
@@ -46,30 +48,32 @@
 
 # use only the points inside the geographic area
 pts = [p for p in coords_to_points(coords) if p.within(area_shape)] # converts to shapely Point
+del coords # not used any more
 
 print('will use %d of %d randomly generated points that are inside geographic area' % (len(pts), N_POINTS))
-coords = points_to_coords(pts) # convert back to simple NumPy coordinate array
 
-del pts
+#%%
 
 #
 # calculate the Voronoi regions, cut them with the geographic area shape and assign the points to them
 #
 
-poly_shapes, pts, poly_to_pt_assignments = voronoi_regions_from_coords(coords, area_shape)
+poly_shapes, poly_to_pt_assignments = voronoi_regions_from_coords(pts, area_shape)
 
+print('Voronoi region to point assignments:')
 print(poly_to_pt_assignments)
 
+#%%
+
 #
 # plotting
 #
 
 fig, ax = subplot_for_map()
 
-plot_voronoi_polys_with_points_in_area(ax, area_shape, poly_shapes, coords, poly_to_pt_assignments,
- point_labels=list(map(str, range(len(coords)))),
- voronoi_labels=list(map(str, range(len(poly_shapes)))))
-#plot_voronoi_polys_with_points_in_area(ax, area_shape, poly_shapes, coords) # monocolor
+plot_voronoi_polys_with_points_in_area(ax, area_shape, poly_shapes, pts, poly_to_pt_assignments,
+ point_labels=list(map(str, range(len(pts)))),
+ voronoi_labels=list(map(str, poly_to_pt_assignments.keys())))
 
 ax.set_title('%d random points and their Voronoi regions in %s' % (len(pts), COUNTRY))
 

examples/vregion_assignment.py

@@ -101,7 +101,7 @@
  p_vertices.extend([finite_pt, far_point])
 
  print(f'> polygon vertices: {p_vertices}')
- p = MultiPoint(p_vertices).convex_hull
+ p = MultiPoint(p_vertices).convex_hull # Voronoi regions are convex
 
  assert p.is_valid and p.is_simple, 'generated polygon is valid and simple'
  region_polys[i_reg] = p

geovoronoi/__init__.py

@@ -8,8 +8,7 @@
 """
 
 
-from ._voronoi import (coords_to_points, points_to_coords, voronoi_regions_from_coords, polygon_lines_from_voronoi,
- polygon_shapes_from_voronoi_lines, assign_points_to_voronoi_polygons,
+from ._voronoi import (coords_to_points, points_to_coords, voronoi_regions_from_coords,
  get_points_to_poly_assignments)
 from ._geom import calculate_polygon_areas
 

geovoronoi/_geom.py

@@ -10,67 +10,67 @@
 import numpy as np
 from shapely.geometry import Polygon
 
-
-def angle_between_pts(a, b, ref_vec=(1.0, 0.0)):
- """
- Angle *theta* between two points (numpy arrays) `a` and `b` in relation to a reference vector `ref_vec`.
- By default, `ref_vec` is the x-axis (i.e. unit vector (1, 0)).
- *theta* is in [0, 2Pi] unless `a` and `b` are very close. In this case this function returns NaN.
- """
- ang = inner_angle_between_vecs(a - b, np.array(ref_vec))
- if not np.isnan(ang) and a[1] < b[1]:
- ang = 2 * np.pi - ang
-
- return ang
-
-
-def inner_angle_between_vecs(a, b):
- """
- Return the inner angle *theta* between numpy vectors `a` and `b`. *theta* is in [0, Pi] if both `a` and `b` are
- not at the origin (0, 0), otherwise this function returns NaN.
- """
- origin = np.array((0, 0))
- if np.allclose(a, origin) or np.allclose(b, origin):
- return np.nan
-
- au = a / np.linalg.norm(a)
- bu = b / np.linalg.norm(b)
- ang = np.arccos(np.clip(np.dot(au, bu), -1.0, 1.0))
-
- return ang
-
-
-def polygon_around_center(points, center=None, fix_nan_angles=True):
- """
- Order numpy array of coordinates `points` around `center` so that they form a valid polygon. Return that as
- shapely `Polygon` object. If no valid polygon can be formed, return `None`.
- If `center` is None (default), use midpoint of `points` as center.
- """
- if center is None:
- center = points.mean(axis=0)
- else:
- center = np.array(center)
-
- # angle between each point in `points` and `center`
- angles = np.apply_along_axis(angle_between_pts, 1, points, b=center)
-
- # sort by angles and generate polygon
- if fix_nan_angles:
- for repl in (0, np.pi):
- tmp_angles = angles.copy()
- tmp_angles[np.isnan(tmp_angles)] = repl
- poly = Polygon(points[np.argsort(tmp_angles)])
- if poly.is_simple and poly.is_valid:
- return poly
-
- return None
- else:
- poly = Polygon(points[np.argsort(angles)])
-
- if poly.is_simple and poly.is_valid:
- return poly
- else:
- return None
+#
+# def angle_between_pts(a, b, ref_vec=(1.0, 0.0)):
+#  """
+#  Angle *theta* between two points (numpy arrays) `a` and `b` in relation to a reference vector `ref_vec`.
+#  By default, `ref_vec` is the x-axis (i.e. unit vector (1, 0)).
+#  *theta* is in [0, 2Pi] unless `a` and `b` are very close. In this case this function returns NaN.
+#  """
+#  ang = inner_angle_between_vecs(a - b, np.array(ref_vec))
+#  if not np.isnan(ang) and a[1] < b[1]:
+#  ang = 2 * np.pi - ang
+#
+#  return ang
+#
+#
+# def inner_angle_between_vecs(a, b):
+#  """
+#  Return the inner angle *theta* between numpy vectors `a` and `b`. *theta* is in [0, Pi] if both `a` and `b` are
+#  not at the origin (0, 0), otherwise this function returns NaN.
+#  """
+#  origin = np.array((0, 0))
+#  if np.allclose(a, origin) or np.allclose(b, origin):
+#  return np.nan
+#
+#  au = a / np.linalg.norm(a)
+#  bu = b / np.linalg.norm(b)
+#  ang = np.arccos(np.clip(np.dot(au, bu), -1.0, 1.0))
+#
+#  return ang
+
+
+# def polygon_around_center(points, center=None, fix_nan_angles=True):
+#  """
+#  Order numpy array of coordinates `points` around `center` so that they form a valid polygon. Return that as
+#  shapely `Polygon` object. If no valid polygon can be formed, return `None`.
+#  If `center` is None (default), use midpoint of `points` as center.
+#  """
+#  if center is None:
+#  center = points.mean(axis=0)
+#  else:
+#  center = np.array(center)
+#
+#  # angle between each point in `points` and `center`
+#  angles = np.apply_along_axis(angle_between_pts, 1, points, b=center)
+#
+#  # sort by angles and generate polygon
+#  if fix_nan_angles:
+#  for repl in (0, np.pi):
+#  tmp_angles = angles.copy()
+#  tmp_angles[np.isnan(tmp_angles)] = repl
+#  poly = Polygon(points[np.argsort(tmp_angles)])
+#  if poly.is_simple and poly.is_valid:
+#  return poly
+#
+#  return None
+#  else:
+#  poly = Polygon(points[np.argsort(angles)])
+#
+#  if poly.is_simple and poly.is_valid:
+#  return poly
+#  else:
+#  return None
 
 
 def calculate_polygon_areas(poly_shapes, m2_to_km2=False):