From fad86de8e6c4ffacac41929bd1bb38982a48d1cf Mon Sep 17 00:00:00 2001 From: Markus Konrad Date: Wed, 3 Feb 2021 17:08:28 +0100 Subject: [PATCH] updated docs --- examples/duplicate_points.py | 4 +- geovoronoi/__init__.py | 2 +- geovoronoi/_voronoi.py | 177 ++++++++++++++++++++------- geovoronoi/plotting.py | 227 +++++++++++++++++++++++------------ tests/_testtools.py | 6 + tests/test_geom.py | 6 + tests/test_main.py | 16 ++- 7 files changed, 312 insertions(+), 126 deletions(-) diff --git a/examples/duplicate_points.py b/examples/duplicate_points.py index 5300a72..a111757 100644 --- a/examples/duplicate_points.py +++ b/examples/duplicate_points.py @@ -14,7 +14,7 @@ import numpy as np import geopandas as gpd -from geovoronoi import coords_to_points, points_to_coords, voronoi_regions_from_coords, get_points_to_poly_assignments +from geovoronoi import coords_to_points, points_to_coords, voronoi_regions_from_coords, points_to_region from geovoronoi.plotting import subplot_for_map, plot_voronoi_polys_with_points_in_area @@ -82,7 +82,7 @@ print('> voronoi region', i_poly, '-> points', str(pt_indices)) print('\n\npoints to voronoi region assignments:') -pts_to_poly_assignments = get_points_to_poly_assignments(poly_to_pt_assignments) +pts_to_poly_assignments = points_to_region(poly_to_pt_assignments) for i_pt, i_poly in pts_to_poly_assignments.items(): print('> point ', i_pt, '-> voronoi region', i_poly) diff --git a/geovoronoi/__init__.py b/geovoronoi/__init__.py index 18a209b..b4d55c8 100644 --- a/geovoronoi/__init__.py +++ b/geovoronoi/__init__.py @@ -9,7 +9,7 @@ from ._voronoi import (coords_to_points, points_to_coords, voronoi_regions_from_coords, - get_points_to_poly_assignments) + points_to_region) from ._geom import calculate_polygon_areas, line_segment_intersection diff --git a/geovoronoi/_voronoi.py b/geovoronoi/_voronoi.py index 84c25e6..374632a 100644 --- a/geovoronoi/_voronoi.py +++ b/geovoronoi/_voronoi.py @@ -133,9 +133,12 @@ def voronoi_regions_from_coords(coords, geo_shape, per_geom=True, return_unassig logger.info('running Voronoi tesselation for %d points / treating geoms separately: %s' % (len(coords), per_geom)) + # check input arguments + if len(coords) < 2: raise ValueError('insufficient number of points provided in `coords`') + # we need the points as NumPy coordinates array and as list of Shapely Point objects if isinstance(coords, np.ndarray): pts = coords_to_points(coords) else: @@ -149,63 +152,75 @@ def voronoi_regions_from_coords(coords, geo_shape, per_geom=True, return_unassig raise ValueError('`geo_shape` is not a valid shape; try applying `.buffer(b)` where `b` is zero ' 'or a very small number') + # get sub-geometries of `geo_shape` if they exist and if we want to treat sub-geometries separately if not per_geom or isinstance(geo_shape, Polygon): geoms = [geo_shape] - else: # Multipolygon + else: # Multipolygon has sub-geometries geoms = geo_shape.geoms - pts_indices = set(range(len(pts))) + # set up data containers + pts_indices = set(range(len(pts))) # point indices to operate on geom_region_polys = {} geom_region_pts = {} unassigned_pts = set() + # iterate through sub-geometries in `geo_shape` for i_geom, geom in enumerate(geoms): + # get point indices of points that lie within `geom` pts_in_geom = [i for i in pts_indices if geom.contains(pts[i])] + # start with empty data for this sub-geometry geom_region_polys[i_geom] = {} geom_region_pts[i_geom] = {} logger.info('%d of %d points in geometry #%d of %d' % (len(pts_in_geom), len(pts), i_geom + 1, len(geoms))) - if not pts_in_geom: + if not pts_in_geom: # no points within `geom` -> skip continue + # remove the points that we're about to use (point - geometry assignment is bijective) pts_indices.difference_update(pts_in_geom) + # generate the Voronoi regions using the SciPy Voronoi class logger.info('generating Voronoi regions') try: vor = Voronoi(coords[pts_in_geom]) - except QhullError as exc: + except QhullError as exc: # handle error for insufficient number of input points by skipping this sub-geom. 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 + raise exc # otherwise re-raise the error + # generate Voronoi region polygons and region -> points mapping logger.info('generating Voronoi region polygons') geom_polys, geom_pts = region_polygons_from_voronoi(vor, geom, return_point_assignments=True, **kwargs) - # map back to original point indices + # point indices in `geom_pts` refer to `coords[pts_in_geom]` -> map them back to original point indices pts_in_geom_arr = np.asarray(pts_in_geom) for i_reg, pt_indices in geom_pts.items(): geom_pts[i_reg] = pts_in_geom_arr[pt_indices].tolist() + # save data for this sub-geometry geom_region_polys[i_geom] = geom_polys geom_region_pts[i_geom] = geom_pts + # collect results logger.info('collecting Voronoi region results') - if results_per_geom: + if results_per_geom: # nothing to do here since we already have results per sub-geometry region_polys = geom_region_polys region_pts = geom_region_pts - else: + else: # merge the results from the sub-geometries to dicts that map region IDs to polygons / point indices + # across all sub-geometries region_polys = {} region_pts = {} for i_geom, geom_polys in geom_region_polys.items(): geom_pts = geom_region_pts[i_geom] + # assign new region IDs region_ids = range(len(region_polys), len(region_polys) + len(geom_polys)) region_ids_mapping = dict(zip(region_ids, geom_polys.keys())) region_polys.update(dict(zip(region_ids, geom_polys.values()))) @@ -221,44 +236,95 @@ def voronoi_regions_from_coords(coords, geo_shape, per_geom=True, return_unassig def region_polygons_from_voronoi(vor, geom, return_point_assignments=False, bounds_buf_factor=0.1, diff_topo_error_buf_factor=0.00000001): + """ + Construct Shapely Polygon/Multipolygon objects for each Voronoi region in `vor` so that they cover the geometry + `geom`. + + :param vor: SciPy Voronoi object + :param geom: Shapely Polygon/MultiPolygon object that confines the constructed Voronoi regions into this shape + :param return_point_assignments: if True, also return a dict that maps Voronoi region IDs to list of point indices + inside that Voronoi region + :param bounds_buf_factor: factor multiplied to largest extend of the `geom` bounding box when calculating the buffer + distance to enlarge the `geom` bounding box + :param diff_topo_error_buf_factor: factor multiplied to largest extend of the `geom` bounding box when calculating + the buffer distance to enlarge `geom` in the rare case of a TopologicalError + :return: dict mapping Voronoi region IDs to Shapely Polygon/Multipolygon objects; if `return_point_assignments` is + True, this function additionally returns a dict that maps Voronoi region IDs to list of point indices + inside that Voronoi region + """ + + # check input arguments + + if not isinstance(vor, Voronoi): + raise ValueError('`vor` must be SciPy Voronoi object') + + if not isinstance(geom, (Polygon, MultiPolygon)): + raise ValueError('`geom` must be a Polygon or MultiPolygon') + + # construct geom bounding box with buffer max_extend = max(geom.bounds[2] - geom.bounds[0], geom.bounds[3] - geom.bounds[1]) bounds_buf = max_extend * bounds_buf_factor geom_bb = box(*geom.bounds).buffer(bounds_buf) + + # center of all points center = np.array(MultiPoint(vor.points).convex_hull.centroid) - ridge_vert = np.array(vor.ridge_vertices) - region_pts = defaultdict(list) - region_neighbor_pts = defaultdict(set) - region_polys = {} + # ridge vertice's coordinates + ridge_vert = np.asarray(vor.ridge_vertices) - logger.debug('generating polygons') + # data containers for output + region_pts = defaultdict(list) # maps region ID to list of point indices + region_neighbor_pts = defaultdict(set) # maps region ID to set of direct neighbor point IDs + region_polys = {} # maps region ID to Shapely Polygon/MultiPolygon object - covered_area = 0 + logger.debug('generating preliminary polygons') + + # this function has two phases: + # - phase 1 generates preliminary polygons from the Voronoi region vertices in `vor`; these polygons may not cover + # the full area of `geom` yet + # - phase 2 fills up the preliminary polygons to fully cover the area of the surrounding `geom` + + # --- phase 1: preliminary polygons --- + + covered_area = 0 # keeps track of area so far covered by generated Voronoi polygons + # 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): - pt_indices, = np.nonzero(vor.point_region == i_reg) + pt_indices, = np.nonzero(vor.point_region == i_reg) # points within this region if len(pt_indices) == 0: # skip regions w/o points in them continue + # update the region -> point mappings region_pts[i_reg].extend(pt_indices.tolist()) + # construct the theoretical Polygon `p` of this region if np.all(np.array(reg_vert) >= 0): # fully finite-bounded region p = Polygon(vor.vertices[reg_vert]) else: - p_vert_indices = set() - p_vert_farpoints = set() - for i_pt in pt_indices: # also consider duplicates + # 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 + p_vert_indices = set() # collect unique indices of vertices into `ridge_vert` + p_vert_farpoints = set() # collect unique calculated finite far points + + # iterate through points inside this region (this is usually only one point, but we also consider + # duplicates) + for i_pt in pt_indices: + # create a mask that selects ridge vertices for when this point lies on either side of the ridge 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]): + # iterate through the point-to-point pairs (`ridge_points`) and the line segments that represent the + # ridge given by ridge vertices + for pointidx, vertidx in zip(vor.ridge_points[enclosing_ridge_pts_mask], + ridge_vert[enclosing_ridge_pts_mask]): + # update the set of neighbor point indices for this region region_neighbor_pts[i_reg].update([x for x in pointidx if x != i_pt]) - if np.all(simplex >= 0): # both vertices of the ridge are finite points - p_vert_indices.update(simplex) + if np.all(vertidx >= 0): # both vertices of the ridge are finite points + p_vert_indices.update(vertidx) # we can add both points else: # "loose ridge": contains infinite Voronoi vertex - # we calculate the far point, i.e. the point of intersection with a surrounding polygon boundary - i = simplex[simplex >= 0][0] # finite end Voronoi vertex + # we calculate the far point, i.e. the point of intersection with the bounding box + i = vertidx[vertidx >= 0][0] # finite end Voronoi vertex finite_pt = vor.vertices[i] p_vert_indices.add(i) @@ -266,9 +332,12 @@ def region_polygons_from_voronoi(vor, geom, return_point_assignments=False, 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 + midpoint = vor.points[pointidx].mean(axis=0) # point in the middle should be directly on the + # ridge + direction = np.sign(np.dot(midpoint - center, n)) * n # re-orient according to center + # find intersection(s) from ridge line going from `midpoint` towards `direction` + # eventually hitting a side of the geom bounding box isects = [] for i_ext_coord in range(len(geom_bb.exterior.coords) - 1): isect = line_segment_intersection(midpoint, direction, @@ -278,15 +347,16 @@ def region_polygons_from_voronoi(vor, geom, return_point_assignments=False, isects.append(isect) if len(isects) == 0: - raise RuntimeError('far point must intersect with surrounding geometry from `geom`') - elif len(isects) == 1: + raise RuntimeError('ridge line must intersect with surrounding geometry from `geom`') + elif len(isects) == 1: # one intersection far_pt = isects[0] - else: + else: # multiple intersections - take closest intersection closest_isect_idx = np.argmin(np.linalg.norm(midpoint - isects, axis=1)) far_pt = isects[closest_isect_idx] + # only use this far point if the found ridge points in the 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 + if (far_pt - finite_pt)[nonzero_coord] / direction[nonzero_coord] > 0: p_vert_farpoints.add(tuple(far_pt)) # create the Voronoi region polygon as convex hull of the ridge vertices and far points (Voronoi regions @@ -301,12 +371,15 @@ def region_polygons_from_voronoi(vor, geom, return_point_assignments=False, p = MultiPoint(p_pts).convex_hull + # check the generated Polygon + if not isinstance(p, Polygon): raise RuntimeError('generated convex hull is not a polygon') if not p.is_valid or p.is_empty: raise RuntimeError('generated polygon is not valid or is empty') + # if the generated preliminary polygon `p` is not completely inside `geom`, then cut it (i.e. use intersection) if not geom.contains(p): p = p.intersection(geom) @@ -314,12 +387,21 @@ def region_polygons_from_voronoi(vor, geom, return_point_assignments=False, raise RuntimeError('generated polygon is not valid or is empty after intersection with the' ' surrounding geometry `geom`') + # update covered area and set the preliminary polygon for this region covered_area += p.area region_polys[i_reg] = p - uncovered_area_portion = (geom.area - covered_area) / geom.area + # --- phase 2: filling preliminary polygons --- + + logger.debug('filling preliminary polygons to fully cover the surrounding area') + + uncovered_area_portion = (geom.area - covered_area) / geom.area # initial portion of uncovered area polys_iter = iter(region_polys.items()) + # 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 while not np.isclose(uncovered_area_portion, 0) and 0 < uncovered_area_portion <= 1: try: i_reg, p = next(polys_iter) @@ -328,13 +410,16 @@ def region_polygons_from_voronoi(vor, geom, return_point_assignments=False, polys_iter = iter(region_polys.items()) i_reg, p = next(polys_iter) pass_ += 1 - else: + else: # no more preliminary Voronoi region polygons left after both passes break logger.debug('will fill up %f%% uncovered area' % (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]) + # generate difference between geom and other regions' polygons -- what's left is the current region's area + # plus any so far uncovered area try: diff = geom.difference(union_other_regions) except TopologicalError: # may happen in rare circumstances @@ -349,17 +434,22 @@ def region_polygons_from_voronoi(vor, geom, return_point_assignments=False, else: diff = [diff] + # list of neighbor region IDs neighbor_regions = [vor.point_region[pt] for pt in region_neighbor_pts[i_reg]] - add = [] - for diff_part in diff: + add = [] # will hold shapes to be added to the current region's polygon + for diff_part in diff: # iterate through the geometries in the difference if diff_part.is_valid and not diff_part.is_empty: if pass_ == 1: + # for first pass, we want an intersection between the so far uncovered area and the current + # region's polygon, so that the uncovered area is directly connected with the region's polygon isect = p.intersection(diff_part) has_isect = isinstance(isect, (Polygon, MultiPolygon)) and not isect.is_empty else: has_isect = True # we don't need an intersection on second pass -- handle isolated features if has_isect: + # we make sure that there's no neighboring region's polygon in between the current region's polygon + # and the area we want to add 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) @@ -368,9 +458,11 @@ def region_polygons_from_voronoi(vor, geom, return_point_assignments=False, else: covered_area += diff_part.area + # add new areas as union if add: region_polys[i_reg] = cascaded_union([p] + add) + # update the portion of uncovered area uncovered_area_portion = (geom.area - covered_area) / geom.area logger.debug('%f%% uncovered area left' % (uncovered_area_portion * 100)) @@ -381,21 +473,20 @@ def region_polygons_from_voronoi(vor, geom, return_point_assignments=False, return region_polys -def get_points_to_poly_assignments(poly_to_pt_assignments): +def points_to_region(region_pts): """ Reverse Voronoi region polygon IDs to point ID assignments by returning a dict that maps each point ID to its Voronoi region polygon ID. All IDs should be integers. - :param poly_to_pt_assignments: dict mapping Voronoi region polygon IDs to list of point IDs + :param region_pts: dict mapping Voronoi region polygon IDs to list of point IDs :return: dict mapping point ID to Voronoi region polygon ID """ - pt_to_poly = {} - for poly, pts in poly_to_pt_assignments.items(): + pts_region = {} + for poly, pts in region_pts.items(): for pt in pts: - if pt in pt_to_poly.keys(): - raise ValueError('invalid assignments in `poly_to_pt_assignments`: ' - 'point %d is assigned to several polygons' % pt) - pt_to_poly[pt] = poly + if pt in pts_region.keys(): + raise ValueError('invalid assignments in `region_pts`: point %d is assigned to several polygons' % pt) + pts_region[pt] = poly - return pt_to_poly + return pts_region diff --git a/geovoronoi/plotting.py b/geovoronoi/plotting.py index 6997164..d00c635 100644 --- a/geovoronoi/plotting.py +++ b/geovoronoi/plotting.py @@ -1,9 +1,6 @@ """ Functions for plotting Voronoi regions. -Most functions rely on [geopandas](http://geopandas.org/) plotting functions. So to use them, this package must be -installed. - Author: Markus Konrad """ @@ -14,13 +11,19 @@ from geopandas.plotting import _flatten_multi_geoms from geopandas import GeoSeries -from ._voronoi import points_to_coords, get_points_to_poly_assignments +from ._voronoi import points_to_coords, points_to_region def subplot_for_map(show_x_axis=False, show_y_axis=False, 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 aspect: aspect ratio + :param kwargs: additional parameters passed to `plt.subplots()` + :return: tuple with (matplotlib Figure, matplotlib Axes) """ fig, ax = plt.subplots(**kwargs) ax.set_aspect(aspect) @@ -38,28 +41,40 @@ def generate_n_colors(n, cmap_name='tab20'): """ Get a list of `n` numbers from matplotlib color map `cmap_name`. If `n` is larger than the number of colors in the color map, the colors will be recycled, i.e. they are not unique in this case. + + :param n: number of colors to generate + :param cmap_name: matplotlib color map name + :return: list of `n` colors """ pt_region_colormap = plt.get_cmap(cmap_name) max_i = len(pt_region_colormap.colors) return [pt_region_colormap(i % max_i) for i in range(n)] -def colors_for_voronoi_polys_and_points(poly_shapes, poly_to_pt_assignments, points=None, cmap_name='tab20'): +def colors_for_voronoi_polys_and_points(region_polys, region_pts, point_indices=None, cmap_name='tab20', + unassigned_pts_color=(0,0,0,1)): """ - Generate colors for the shapes and points in `poly_shapes` and `poly_to_pt_assignments` using matplotlib color + Generate colors for the shapes and points in `region_polys` and `region_pts` using matplotlib color map `cmap_name`. + + :param region_polys: dict mapping region IDs to Voronoi region geometries + :param region_pts: dict mapping Voronoi region IDs to point indices + :param point_indices: optional list of point indices; used to set color to `unassigned_pts_color` for unused point + indices + :param cmap_name: matplotlib color map name + :return: tuple with (dict mapping Voronoi region ID to color, list of point colors) """ vor_colors = {p_id: col - for p_id, col in zip(poly_shapes.keys(), generate_n_colors(len(poly_shapes), cmap_name=cmap_name))} + for p_id, col in zip(region_polys.keys(), generate_n_colors(len(region_polys), cmap_name=cmap_name))} - pt_to_poly = get_points_to_poly_assignments(poly_to_pt_assignments) + pt_to_poly = points_to_region(region_pts) - if points is not None: - pt_to_poly.update({i_pt: None for i_pt in range(len(points)) if i_pt not in pt_to_poly.keys()}) + if point_indices is not None: + pt_to_poly.update({i_pt: None for i_pt in point_indices if i_pt not in pt_to_poly.keys()}) pt_to_poly = sorted(pt_to_poly.items(), key=lambda x: x[0]) - pt_colors = [(0,0,0,1) if i_vor is None else vor_colors[i_vor] for _, i_vor in pt_to_poly] + pt_colors = [unassigned_pts_color if i_vor is None else vor_colors[i_vor] for _, i_vor in pt_to_poly] assert len(vor_colors) <= len(pt_colors) @@ -67,6 +82,13 @@ def colors_for_voronoi_polys_and_points(poly_shapes, poly_to_pt_assignments, poi def xy_from_points(points): + """ + Helper function to get x and y coordinates from NumPy array or list of Shapely Points `points` especially used + for plotting. + + :param points: NumPy array or list of Shapely Points with length N + :return: tuple with (NumPy array of N x-coordinates, NumPy array of N y-coordinates) + """ if hasattr(points, 'xy'): return points.xy else: @@ -77,26 +99,38 @@ def xy_from_points(points): return coords[:, 0], coords[:, 1] -def plot_voronoi_polys(ax, poly_shapes, color=None, edgecolor=None, labels=None, label_fontsize=10, label_color=None, +def plot_voronoi_polys(ax, region_polys, color=None, edgecolor=None, labels=None, label_fontsize=10, label_color=None, **kwargs): """ Plot Voronoi region polygons in `poly_shapes` on matplotlib Axes object `ax`. Use fill color `color`, edge color `edgecolor`. Optionally supply a list of labels `labels` that will be displayed on the respective Voronoi region - using the styling options `label_fontsize` and `label_color`. All color parameters can also be a sequence. - Additional parameters can be passed to GeoPandas' `_plot_polygon_collection_with_color` function as `kwargs`. + using the styling options `label_fontsize` and `label_color`. All color parameters can also be a dict or sequence + mapped to the Voronoi region accordingly. + + Additional parameters can be passed to `plot_polygon_collection_with_color()` function as `kwargs`. + + :param ax: matplotlib Axes object to plot on + :param region_polys: dict mapping region IDs to Voronoi region geometries + :param color: region polygons' fill color + :param edgecolor: region polygons' edge color + :param labels: region labels + :param label_fontsize: region labels font size + :param label_color: region labels color + :param kwargs: additional parameters passed to `plot_polygon_collection_with_color()` function + :return: None """ - _plot_polygon_collection_with_color(ax, poly_shapes, color=color, edgecolor=edgecolor, **kwargs) + plot_polygon_collection_with_color(ax, region_polys, color=color, edgecolor=edgecolor, **kwargs) if labels: # plot labels using matplotlib's text() n_labels = len(labels) - n_features = len(poly_shapes) + n_features = len(region_polys) if n_labels != n_features: raise ValueError('number of labels (%d) must match number of Voronoi polygons (%d)' % (n_labels, n_features)) - for i, p in poly_shapes.items(): + for i, p in region_polys.items(): tx, ty = p.centroid.coords[0] ax.text(tx, ty, labels[i], fontsize=label_fontsize, color=_color_for_labels(label_color, color, i)) @@ -108,7 +142,20 @@ def plot_points(ax, points, markersize=1, marker='o', color=None, labels=None, l marker size `markersize`. Define marker style with parameters `marker` and `color`. Optionally supply a list of labels `labels` that will be displayed next to the respective point using the styling options `label_fontsize` and `label_color`. All color parameters can also be a sequence. + Additional parameters can be passed to matplotlib's `scatter` function as `kwargs`. + + :param ax: matplotlib Axes object to plot on + :param points: NumPy array or list of Shapely Points + :param markersize: marker size + :param marker: marker type + :param color: marker color(s) + :param labels: point labels + :param label_fontsize: point labels font size + :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 """ x, y = xy_from_points(points) @@ -130,13 +177,17 @@ def plot_points(ax, points, markersize=1, marker='o', color=None, labels=None, l drawn_coords.add(pos) -def plot_lines(ax, points, linewidth=1, color=None, **kwargs): +def plot_line(ax, points, linewidth=1, color=None, **kwargs): """ - Plot points `points` (either list of Point objects or NumPy coordinate array) on matplotlib Axes object `ax` with - marker size `markersize`. Define marker style with parameters `marker` and `color`. Optionally supply a list of - labels `labels` that will be displayed next to the respective point using the styling options `label_fontsize` and - `label_color`. All color parameters can also be a sequence. - Additional parameters can be passed to matplotlib's `scatter` function as `kwargs`. + Plot sequence of points `points` (either list of Point objects or NumPy coordinate array) on matplotlib Axes object + `ax` as line. + + :param ax: matplotlib Axes object to plot on + :param points: NumPy array or list of Shapely Point objects + :param linewidth: line width + :param color: line color + :param kwargs: additional parameters passed to matplotlib's `plot` function + :return: None """ x, y = xy_from_points(points) @@ -145,6 +196,21 @@ def plot_lines(ax, points, linewidth=1, color=None, **kwargs): def plot_polygon(ax, polygon, facecolor=None, edgecolor=None, linewidth=1, linestyle='solid', label=None, label_fontsize=7, label_color='k', **kwargs): + """ + Plot a Shapely polygon on matplotlib Axes object `ax`. + + :param ax: matplotlib Axes object to plot on + :param polygon: Shapely Polygon/Multipolygon object + :param facecolor: fill color of the polygon + :param edgecolor: edge color of the polygon + :param linewidth: line width + :param linestyle: line style + :param label: optional label to show at polygon's centroid + :param label_fontsize: label font size + :param label_color: label color + :param kwargs: additonal parameters passed to matplotlib `PatchCollection` object + :return: None + """ ax.add_collection(PatchCollection([PolygonPatch(polygon)], facecolor=facecolor, edgecolor=edgecolor, linewidth=linewidth, linestyle=linestyle, @@ -153,7 +219,7 @@ def plot_polygon(ax, polygon, facecolor=None, edgecolor=None, linewidth=1, lines ax.text(polygon.centroid.x, polygon.centroid.y, label, fontsize=label_fontsize, color=label_color) -def plot_voronoi_polys_with_points_in_area(ax, area_shape, poly_shapes, points, poly_to_pt_assignments=None, +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), voronoi_and_points_cmap='tab20', voronoi_color=None, voronoi_edgecolor=None, @@ -164,30 +230,56 @@ def plot_voronoi_polys_with_points_in_area(ax, area_shape, poly_shapes, points, plot_voronoi_opts=None, plot_points_opts=None): """ - All-in-one function to plot Voronoi region polygons `poly_shapes` and the respective points `points` inside a - geographic area `area_shape` on a matplotlib Axes object `ax`. By default, the regions will be blue and the points - black. Optionally pass `poly_to_pt_assignments` to show Voronoi regions and their respective points with the same - color (which is randomly drawn from color map `voronoi_and_points_cmap`). Labels for Voronoi regions can be passed - as `voronoi_labels`. Labels for points can be passed as `point_labels`. Use style options to customize the plot. - Pass additional (matplotlib) parameters to the individual plotting steps as `plot_area_opts`, `plot_voronoi_opts` or - `plot_points_opts` respectively. + All-in-one function to plot Voronoi region polygons `region_polys` and the respective points `points` inside a + geographic area `area_shape` on a matplotlib Axes object `ax`. + + By default, the regions will be blue and the points black. Optionally pass `region_pts` to show Voronoi regions and + their respective points with the same color (which is randomly drawn from color map `voronoi_and_points_cmap`). + Labels for Voronoi regions can be passed as `voronoi_labels`. Labels for points can be passed as `point_labels`. + Use style options to customize the plot. Pass additional (matplotlib) parameters to the individual plotting steps + as `plot_area_opts`, `plot_voronoi_opts` or `plot_points_opts` respectively. + + :param ax: matplotlib Axes object to plot on + :param area_shape: geographic shape surrounding the Voronoi regions + :param region_polys: dict mapping region IDs to Voronoi region geometries + :param points: NumPy array or list of Shapely Point objects + :param region_pts: dict mapping Voronoi region IDs to point indices of `points` + :param area_color: fill color of `area_shape` + :param area_edgecolor: edge color of `area_shape` + :param voronoi_and_points_cmap: matplotlib color map name used for Voronoi regions and points colors when colors are + not given by `voronoi_color` + :param voronoi_color: dict mapping Voronoi region ID to fill color or None to use `voronoi_and_points_cmap` + :param voronoi_edgecolor: Voronoi region polygon edge colors + :param points_color: points color + :param points_markersize: points marker size + :param points_marker: points marker type + :param voronoi_labels: Voronoi region labels displayed at centroid of Voronoi region polygon + :param voronoi_label_fontsize: Voronoi region labels font size + :param voronoi_label_color: Voronoi region labels color + :param point_labels: point labels + :param point_label_fontsize: point labels font size + :param point_label_color: point labels color + :param plot_area_opts: options passed to function for plotting the geographic shape + :param plot_voronoi_opts: options passed to function for plotting the Voronoi regions + :param plot_points_opts: options passed to function for plotting the points + :return: None """ plot_area_opts = plot_area_opts or {} plot_voronoi_opts = plot_voronoi_opts or {'alpha': 0.5} plot_points_opts = plot_points_opts or {} - _plot_polygon_collection_with_color(ax, [area_shape], color=area_color, edgecolor=area_edgecolor, **plot_area_opts) + plot_polygon_collection_with_color(ax, [area_shape], color=area_color, edgecolor=area_edgecolor, **plot_area_opts) - if voronoi_and_points_cmap and poly_to_pt_assignments and \ + if voronoi_and_points_cmap and region_pts and \ not all(map(bool, (voronoi_color, voronoi_edgecolor, points_color))): - voronoi_color, points_color = colors_for_voronoi_polys_and_points(poly_shapes, poly_to_pt_assignments, - points=points, + voronoi_color, points_color = colors_for_voronoi_polys_and_points(region_polys, region_pts, + point_indices=list(range(len(points))), cmap_name=voronoi_and_points_cmap) if voronoi_color is None and voronoi_edgecolor is None: voronoi_edgecolor = (0,0,0,1) # better visible default value - plot_voronoi_polys(ax, poly_shapes, color=voronoi_color, edgecolor=voronoi_edgecolor, + plot_voronoi_polys(ax, region_polys, color=voronoi_color, edgecolor=voronoi_edgecolor, labels=voronoi_labels, label_fontsize=voronoi_label_fontsize, label_color=voronoi_label_color, **plot_voronoi_opts) @@ -196,47 +288,19 @@ def plot_voronoi_polys_with_points_in_area(ax, area_shape, poly_shapes, points, **plot_points_opts) -def _color_for_labels(label_color, default_color, seq_index): - """Helper function to get a color for a label with index `seq_index`.""" - if label_color is None: - if hasattr(default_color, '__getitem__'): - c = default_color[seq_index] - else: - c = default_color - else: - c = label_color - - return c or (0,0,0,1) - - -def _plot_polygon_collection_with_color(ax, geoms, color=None, **kwargs): +def plot_polygon_collection_with_color(ax, geoms, color=None, **kwargs): """ - This is a modified version of geopanda's `plot_polygon_collection` function that also accepts a sequences of colors - passed as `color` for each polygon in `geoms` and *uses them correctly even when `geoms` contains MultiPolygon - objects*. - - Plots a collection of Polygon and MultiPolygon geometries to `ax` - - Parameters - ---------- + This is a modified version of geopanda's `_plot_polygon_collection` function that also accepts a sequence or dict + of colors passed as `color` for each polygon in `geoms` and *uses them correctly even when `geoms` contains + MultiPolygon objects*. - ax : matplotlib.axes.Axes - where shapes will be plotted + Plots a collection of Polygon and MultiPolygon geometries to `ax`. - geoms : a sequence of `N` Polygons and/or MultiPolygons (can be mixed) - - color : a sequence of `N` colors (optional) or None for default color or a single color for all shapes - - edgecolor : single color or sequence of `N` colors - Color for the edge of the polygons - - **kwargs - Additional keyword arguments passed to the collection - - Returns - ------- - - collection : matplotlib.collections.Collection that was plotted + :param ax: matplotlib.axes.Axes object where shapes will be plotted + :param geoms: a sequence or dict of `N` Polygons and/or MultiPolygons (can be mixed) + :param color: a sequence or dict of `N` colors (optional) or None for default color or a single color for all shapes + :param kwargs: additional keyword arguments passed to the collection + :return: matplotlib.collections.Collection that was plotted """ color_values = color @@ -266,7 +330,7 @@ def _plot_polygon_collection_with_color(ax, geoms, color=None, **kwargs): geoms, multi_indices = _flatten_multi_geoms(geoms) - if color_indices is not None: + if color_indices is not None: # retain correct color indices color_values = color_values[np.nonzero(geoms_indices[multi_indices][..., np.newaxis] == color_indices)[1]] # PatchCollection does not accept some kwargs. @@ -279,3 +343,16 @@ def _plot_polygon_collection_with_color(ax, geoms, color=None, **kwargs): ax.add_collection(collection, autolim=True) ax.autoscale_view() return collection + + +def _color_for_labels(label_color, default_color, seq_index): + """Helper function to get a color for a label with index `seq_index`.""" + if label_color is None: + if hasattr(default_color, '__getitem__'): + c = default_color[seq_index] + else: + c = default_color + else: + c = label_color + + return c or (0,0,0,1) diff --git a/tests/_testtools.py b/tests/_testtools.py index 9b11188..227a5af 100644 --- a/tests/_testtools.py +++ b/tests/_testtools.py @@ -1,3 +1,9 @@ +""" +Set of utility functions for testing. + +Author: Markus Konrad +""" + from functools import partial import numpy as np diff --git a/tests/test_geom.py b/tests/test_geom.py index 0be17f2..7f535fa 100644 --- a/tests/test_geom.py +++ b/tests/test_geom.py @@ -1,3 +1,9 @@ +""" +Tests for the private _geom module. + +Author: Markus Konrad +""" + import pytest import numpy as np diff --git a/tests/test_main.py b/tests/test_main.py index eebc93c..a50c2d7 100644 --- a/tests/test_main.py +++ b/tests/test_main.py @@ -1,3 +1,9 @@ +""" +Tests for the main module. + +Author: Markus Konrad +""" + import numpy as np import geopandas as gpd from shapely.geometry import Polygon, MultiPolygon, asPoint @@ -9,7 +15,7 @@ from ._testtools import coords_2d_array from geovoronoi import ( voronoi_regions_from_coords, coords_to_points, points_to_coords, calculate_polygon_areas, - get_points_to_poly_assignments + points_to_region ) from geovoronoi.plotting import subplot_for_map, plot_voronoi_polys_with_points_in_area @@ -35,9 +41,9 @@ def test_coords_to_points_and_points_to_coords(coords): def test_get_points_to_poly_assignments(poly_to_pts, expected): if expected is None: with pytest.raises(ValueError): - get_points_to_poly_assignments(poly_to_pts) + points_to_region(poly_to_pts) else: - assert get_points_to_poly_assignments(poly_to_pts) == expected + assert points_to_region(poly_to_pts) == expected @given(available_points=st.permutations(list(range(10))), n_poly=st.integers(0, 10)) @@ -64,7 +70,7 @@ def test_get_points_to_poly_assignments_hypothesis(available_points, n_poly): if n_poly > 0: assert set(sum(list(poly_to_pts.values()), [])) == set(available_points) - pts_to_poly = get_points_to_poly_assignments(poly_to_pts) + pts_to_poly = points_to_region(poly_to_pts) assert isinstance(pts_to_poly, dict) @@ -154,7 +160,7 @@ def test_voronoi_regions_from_coords_italy(n_pts, per_geom, return_unassigned_pt assert len(region_pts) == len(region_polys) # points to region assignments - pts_region = get_points_to_poly_assignments(region_pts) + pts_region = points_to_region(region_pts) if unassigned_pts is not None: # check that unassigned points are not in the result set assert set(range(n_pts)) - set(pts_region.keys()) == unassigned_pts