diff --git a/examples/duplicate_points.py b/examples/duplicate_points.py index b7337f9..682398f 100644 --- a/examples/duplicate_points.py +++ b/examples/duplicate_points.py @@ -4,7 +4,7 @@ Duplicate points, i.e. points with exactly the same coordinates will belong to the same Voronoi region. Author: Markus Konrad -January 2020 +January 2021 """ diff --git a/examples/random_points_across_italy.py b/examples/random_points_across_italy.py index 7a32297..78e7884 100644 --- a/examples/random_points_across_italy.py +++ b/examples/random_points_across_italy.py @@ -3,7 +3,7 @@ and their points will be plotted using the `plotting` sub-module of `geovoronoi`. Author: Markus Konrad -March 2018 +January 2021 """ diff --git a/examples/random_points_and_area.png b/examples/random_points_and_area.png index f1ad380..b9d0997 100644 Binary files a/examples/random_points_and_area.png and b/examples/random_points_and_area.png differ diff --git a/examples/random_points_and_area.py b/examples/random_points_and_area.py index c71584d..022a585 100644 --- a/examples/random_points_and_area.py +++ b/examples/random_points_and_area.py @@ -6,10 +6,11 @@ Note that it is important to use an *equal area* projection before calculating the areas of the Voronoi regions! Author: Markus Konrad -March 2018 +January 2021 """ import logging +from pprint import pprint import matplotlib.pyplot as plt import numpy as np @@ -24,6 +25,8 @@ geovoronoi_log.setLevel(logging.INFO) geovoronoi_log.propagate = True +#%% + N_POINTS = 20 COUNTRY = 'Spain' @@ -48,26 +51,31 @@ # 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 +n_pts = len(pts) -print('will use %d of %d randomly generated points that are inside geographic area' % (len(pts), N_POINTS)) +print('will use %d of %d randomly generated points that are inside geographic area' % (n_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(coords, area_shape) # calculate area in km², too poly_areas = calculate_polygon_areas(poly_shapes, m2_to_km2=True) # converts m² to km² print('areas in km²:') -print(poly_areas) +pprint(poly_areas) print('sum:') -print(sum(poly_areas)) +print(sum(poly_areas.values())) + +#%% # # plotting @@ -75,12 +83,12 @@ fig, ax = subplot_for_map(show_x_axis=True, show_y_axis=True) -voronoi_labels = ['%d km²' % round(a) for a in poly_areas] +voronoi_labels = {poly_i: '%d km²' % round(a) for poly_i, a in poly_areas.items()} plot_voronoi_polys_with_points_in_area(ax, area_shape, poly_shapes, coords, poly_to_pt_assignments, voronoi_labels=voronoi_labels, voronoi_label_fontsize=7, voronoi_label_color='gray') -ax.set_title('%d random points and their Voronoi regions in %s' % (len(pts), COUNTRY)) +ax.set_title('%d random points and their Voronoi regions in %s' % (n_pts, COUNTRY)) plt.tight_layout() plt.savefig('random_points_and_area.png') diff --git a/examples/toyexample.py b/examples/toyexample.py index 35e9f85..38ab652 100644 --- a/examples/toyexample.py +++ b/examples/toyexample.py @@ -1,3 +1,11 @@ +""" +Artificial mini-example useful for debugging. + +Author: Markus Konrad +January 2021 +""" + + import logging import matplotlib.pyplot as plt @@ -21,16 +29,17 @@ points = np.array([[1, 1], [1.1, 0.9], [1.15, 0.8], [2.5, 0]]) +# surrounding shape shape = Polygon([[-1, -1], [3, -1], [3, 3], [-1, 3]]) -#%% +#%% Voronoi region generation region_polys, region_pts = voronoi_regions_from_coords(points, shape) print('Voronoi region to point assignments:') print(region_pts) -#%% +#%% plotting fig, ax = subplot_for_map(show_x_axis=True, show_y_axis=True) diff --git a/examples/using_geopandas.png b/examples/using_geopandas.png index 1272452..7f6dd38 100644 Binary files a/examples/using_geopandas.png and b/examples/using_geopandas.png differ diff --git a/examples/using_geopandas.py b/examples/using_geopandas.py index ac70c55..2fea3fa 100644 --- a/examples/using_geopandas.py +++ b/examples/using_geopandas.py @@ -2,7 +2,7 @@ Example script to show how to incorporate GeoPandas with geovoronoi. Author: Markus Konrad -March 2018 +January 2021 """ import logging @@ -20,6 +20,8 @@ geovoronoi_log.setLevel(logging.INFO) geovoronoi_log.propagate = True +#%% + # # load geo data # @@ -35,6 +37,7 @@ south_am_shape = cascaded_union(south_am.geometry) south_am_cities = cities[cities.geometry.within(south_am_shape)] # reduce to cities in South America +#%% # # calculate the Voronoi regions, cut them with the geographic area shape and assign the points to them @@ -44,8 +47,10 @@ coords = points_to_coords(south_am_cities.geometry) # calculate the regions -poly_shapes, pts, poly_to_pt_assignments = voronoi_regions_from_coords(coords, south_am_shape) +poly_shapes, poly_to_pt_assignments = voronoi_regions_from_coords(coords, south_am_shape) + +#%% # # Plotting @@ -53,7 +58,8 @@ fig, ax = subplot_for_map() -plot_voronoi_polys_with_points_in_area(ax, south_am_shape, poly_shapes, pts, poly_to_pt_assignments) +plot_voronoi_polys_with_points_in_area(ax, south_am_shape, poly_shapes, coords, poly_to_pt_assignments, + point_labels=south_am_cities.name.tolist()) ax.set_title('Cities data for South America from GeoPandas\nand Voronoi regions around them') diff --git a/examples/vregion_assignment.py b/examples/vregion_assignment.py deleted file mode 100644 index 9ccb45e..0000000 --- a/examples/vregion_assignment.py +++ /dev/null @@ -1,212 +0,0 @@ -from collections import defaultdict -from matplotlib import pyplot as plt - -from scipy.spatial import Voronoi, voronoi_plot_2d -from shapely.geometry import LineString, LinearRing, asPoint, MultiPoint, Polygon -from shapely.ops import polygonize - -import numpy as np - -from geovoronoi._geom import polygon_around_center -from geovoronoi.plotting import subplot_for_map, plot_points, plot_lines, plot_polygon - -#%% - -geo_shape = Polygon([[-1, -1], [3, -1], [3, 3], [-1, 3], [-1, -1]]) - -points = np.array([[0, 0], [1, 0], [2, 0], - [0, 1], [1, 1], [2, 1], - [0, 2], [1, 2], [2, 2]]) - -fig, ax = subplot_for_map(show_x_axis=True, show_y_axis=True) - -plot_polygon(ax, geo_shape, facecolor=(0, 0, 0, 0.0), edgecolor='r') -plot_points(ax, points, color='r', markersize=1, labels=list(map(str, range(len(points)))), label_color='r') - -fig.show() - - -#%% - -vor = Voronoi(points) - -vorfig = voronoi_plot_2d(vor) -vorfig.show() -del vorfig - -vor.vertices -vor.regions -vor.point_region -vor.ridge_vertices -vor.ridge_points -vor.ridge_dict - -ridge_vert = np.array(vor.ridge_vertices) - -#%% - -farpoints_max_extend_factor = 1 -max_dim_extend = vor.points.ptp(axis=0).max() * farpoints_max_extend_factor -center = np.array(MultiPoint(vor.points).convex_hull.centroid) - -region_pts = defaultdict(list) -region_polys = {} -for i_reg, reg_vert in enumerate(vor.regions): - pt_indices = np.where(vor.point_region == i_reg)[0] - if len(pt_indices) == 0: # skip regions w/o points in them - print(i_reg, 'empty') - continue - - print(i_reg) - print('> vertice indices: ', reg_vert) - print('> point indices: ', pt_indices) - - region_pts[i_reg].extend(pt_indices) - - if np.all(np.array(reg_vert) >= 0): # fully finite-bounded region - print('> finite') - p = Polygon(vor.vertices[reg_vert]) - else: - print('> not finite') - - p_vertices = [] - i_pt = pt_indices[0] # only consider one point, not a duplicate - print('>> point', i_pt) - enclosing_ridge_pts_mask = (vor.ridge_points[:, 0] == i_pt) | (vor.ridge_points[:, 1] == i_pt) - for pointidx, simplex in zip(vor.ridge_points[enclosing_ridge_pts_mask], - ridge_vert[enclosing_ridge_pts_mask]): - print('>> pointidx', pointidx) - print('>> simplex indices', simplex) - - if np.all(simplex >= 0): - p_vertices.extend(vor.vertices[simplex]) - else: - # "loose ridge": contains infinite Voronoi vertex - i = simplex[simplex >= 0][0] # finite end Voronoi vertex - finite_pt = vor.vertices[i] - - 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 - direction = direction / np.linalg.norm(direction) # to unit vector - - far_point = finite_pt + direction * max_dim_extend * farpoints_max_extend_factor - - print(f'>> finite point {i} is {finite_pt}') - print(f'>> far point is {far_point}') - - p_vertices.extend([finite_pt, far_point]) - - print(f'> polygon vertices: {p_vertices}') - 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 - - -#%% - -list(region_polys[1].exterior.coords) - -#%% - -fig, ax = subplot_for_map(show_x_axis=True, show_y_axis=True) - -plot_polygon(ax, geo_shape, facecolor=(0, 0, 0, 0.0), edgecolor='r') -plot_points(ax, points, color='r', markersize=1, labels=list(map(str, range(len(points)))), label_color='r') - -for i_p, p in region_polys.items(): - plot_polygon(ax, p, edgecolor='k', facecolor=(0, 0, 0, 0.1), linestyle='dashed', label=str(i_p)) - - -fig.show() - -region_pts - -#%% - -farpoints_max_extend_factor = 10 -max_dim_extend = vor.points.ptp(axis=0).max() * farpoints_max_extend_factor -center = np.array(MultiPoint(vor.points).convex_hull.centroid) - -# generate lists of full polygon lines, loose ridges and far points of loose ridges from scipy Voronoi result object -poly_lines = [] -loose_ridges = [] -far_points = [] -for pointidx, simplex in zip(vor.ridge_points, vor.ridge_vertices): - simplex = np.asarray(simplex) - print(f'pointidx: {pointidx}, simplex: {simplex}') - if np.all(simplex >= 0): # full finite polygon line - print('> finite') - poly_lines.append(LineString(vor.vertices[simplex])) - else: # "loose ridge": contains infinite Voronoi vertex - print('> loose') - i = simplex[simplex >= 0][0] # finite end Voronoi vertex - 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 - direction = direction / np.linalg.norm(direction) # to unit vector - - far_point = vor.vertices[i] + direction * max_dim_extend * farpoints_max_extend_factor - print(f'> finite vertex {i}: {vor.vertices[i]}') - print(f'> far point: {far_point}') - loose_ridges.append(LineString(np.vstack((vor.vertices[i], far_point)))) - far_points.append(far_point) - -# -# confine the infinite Voronoi regions by constructing a "hull" from loose ridge far points around the centroid of -# the geographic area -# - -# first append the loose ridges themselves -for l in loose_ridges: - poly_lines.append(l) - -# now create the "hull" of far points: `far_points_hull` -far_points = np.array(far_points) -far_points_hull = polygon_around_center(far_points, center) - -if far_points_hull is None: - raise RuntimeError('no polygonal hull of far points could be created') - -# sometimes, this hull does not completely encompass the geographic area `geo_shape` -if not far_points_hull.contains(geo_shape): # if that's the case, merge it by taking the union - far_points_hull = far_points_hull.union(geo_shape) - -if not isinstance(far_points_hull, Polygon): - raise RuntimeError('hull of far points is not Polygon as it should be; try increasing ' - '`farpoints_max_extend_factor`') - -# now add the lines that make up `far_points_hull` to the final `poly_lines` list -far_points_hull_coords = far_points_hull.exterior.coords -for i, pt in list(enumerate(far_points_hull_coords))[1:]: - poly_lines.append(LineString((far_points_hull_coords[i - 1], pt))) -poly_lines.append(LineString((far_points_hull_coords[-1], far_points_hull_coords[0]))) - -#%% - -poly_shapes = [] -for p in polygonize(poly_lines): - if geo_shape is not None and not geo_shape.contains(p): # if `geo_shape` contains polygon `p`, - p = p.intersection(geo_shape) # intersect it with `geo_shape` (i.e. "cut" it) - - if not p.is_empty: - poly_shapes.append(p) - -for p in poly_shapes: - print(list(p.exterior.coords)) - -len(poly_shapes) - -#%% - -for p in poly_shapes: - plot_polygon(ax, p, edgecolor='k', facecolor=(0, 0, 0, 0.1), linestyle='dashed') - -fig.show() \ No newline at end of file diff --git a/geovoronoi/_geom.py b/geovoronoi/_geom.py index e263bb5..ea0c92c 100644 --- a/geovoronoi/_geom.py +++ b/geovoronoi/_geom.py @@ -76,8 +76,6 @@ def calculate_polygon_areas(poly_shapes, m2_to_km2=False): Return the area of the respective polygons in `poly_shapes`. Returns a NumPy array of areas in m² (if `m2_to_km2` is False) or km² (otherwise). """ - areas = np.array([p.area for p in poly_shapes]) - if m2_to_km2: - return areas / 1000000 # = 1000² - else: - return areas + + 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/plotting.py b/geovoronoi/plotting.py index d45d0ef..7cc26dd 100644 --- a/geovoronoi/plotting.py +++ b/geovoronoi/plotting.py @@ -96,9 +96,9 @@ def plot_voronoi_polys(ax, poly_shapes, color=None, edgecolor=None, labels=None, raise ValueError('number of labels (%d) must match number of Voronoi polygons (%d)' % (n_labels, n_features)) - for (i, p), lbl in zip(poly_shapes.items(), labels): + for i, p in poly_shapes.items(): tx, ty = p.centroid.coords[0] - ax.text(tx, ty, lbl, fontsize=label_fontsize, color=_color_for_labels(label_color, color, i)) + ax.text(tx, ty, labels[i], fontsize=label_fontsize, color=_color_for_labels(label_color, color, i)) def plot_points(ax, points, markersize=1, marker='o', color=None, labels=None, label_fontsize=7, label_color=None, @@ -123,10 +123,10 @@ def plot_points(ax, points, markersize=1, marker='o', color=None, labels=None, l % (n_labels, n_features)) drawn_coords = set() - for i, (x_i, y_i, lbl) in enumerate(zip(x, y, labels)): + for i, (x_i, y_i) in enumerate(zip(x, y)): pos = (x_i, y_i) # make hashable if label_draw_duplicates or pos not in drawn_coords: - ax.text(x_i, y_i, lbl, fontsize=label_fontsize, color=_color_for_labels(label_color, color, i)) + ax.text(x_i, y_i, labels[i], fontsize=label_fontsize, color=_color_for_labels(label_color, color, i)) drawn_coords.add(pos)