diff --git a/examples/duplicate_points.png b/examples/duplicate_points.png index e62c854..5f41ea5 100644 Binary files a/examples/duplicate_points.png and b/examples/duplicate_points.png differ diff --git a/examples/duplicate_points.py b/examples/duplicate_points.py index f17be39..b7337f9 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 -March 2018 +January 2020 """ @@ -23,6 +23,9 @@ geovoronoi_log.setLevel(logging.INFO) geovoronoi_log.propagate = True + +#%% + N_POINTS = 20 N_DUPL = 10 COUNTRY = 'Sweden' @@ -61,11 +64,8 @@ print('duplicated %d random points -> we have %d coordinates now' % (N_DUPL, len(coords))) -# if we didn't know in advance how many duplicates we have (and which points they are), we could find out like this: -# unique_coords, unique_ind, dupl_counts = np.unique(coords, axis=0, return_index=True, return_counts=True) -# n_dupl = len(coords) - len(unique_ind) -# n_dupl -# >>> 10 + +#%% # # calculate the Voronoi regions, cut them with the geographic area shape and assign the points to them @@ -73,44 +73,46 @@ # the duplicate coordinates will belong to the same voronoi region # -poly_shapes, pts, poly_to_pt_assignments = voronoi_regions_from_coords(coords, area_shape, - accept_n_coord_duplicates=N_DUPL) +poly_shapes, poly_to_pt_assignments = voronoi_regions_from_coords(coords, area_shape) # poly_to_pt_assignments is a nested list because a voronoi region might contain several (duplicate) points print('\n\nvoronoi region to points assignments:') -for i_poly, pt_indices in enumerate(poly_to_pt_assignments): +for i_poly, pt_indices in poly_to_pt_assignments.items(): print('> voronoi region', i_poly, '-> points', str(pt_indices)) print('\n\npoints to voronoi region assignments:') -pts_to_poly_assignments = np.array(get_points_to_poly_assignments(poly_to_pt_assignments)) -for i_pt, i_poly in enumerate(pts_to_poly_assignments): +pts_to_poly_assignments = get_points_to_poly_assignments(poly_to_pt_assignments) +for i_pt, i_poly in pts_to_poly_assignments.items(): print('> point ', i_pt, '-> voronoi region', i_poly) +#%% + # # plotting # -# make point labels: counts of duplicates per points -count_per_pt = [sum(pts_to_poly_assignments == i_poly) for i_poly in pts_to_poly_assignments] -pt_labels = list(map(str, count_per_pt)) +# make point labels: counts of duplicate assignments per points +count_per_pt = {pt_indices[0]: len(pt_indices) for pt_indices in poly_to_pt_assignments.values()} +pt_labels = list(map(str, count_per_pt.values())) +distinct_pt_coords = coords[np.asarray(list(count_per_pt.keys()))] # highlight voronoi regions with point duplicates -count_per_poly = np.array(list(map(len, poly_to_pt_assignments))) -vor_colors = np.repeat('blue', len(poly_shapes)) # default color -vor_colors[count_per_poly > 1] = 'red' # hightlight color +vor_colors = {i_poly: 'red' if len(pt_indices) > 1 else 'blue' + for i_poly, pt_indices in poly_to_pt_assignments.items()} fig, ax = subplot_for_map() -plot_voronoi_polys_with_points_in_area(ax, area_shape, poly_shapes, coords, +plot_voronoi_polys_with_points_in_area(ax, area_shape, poly_shapes, distinct_pt_coords, plot_voronoi_opts={'alpha': 0.2}, plot_points_opts={'alpha': 0.4}, - voronoi_color=list(vor_colors), + voronoi_color=vor_colors, + voronoi_edgecolor='black', point_labels=pt_labels, - points_markersize=np.array(count_per_pt)*10) + points_markersize=np.square(np.array(list(count_per_pt.values())))*10) -ax.set_title('%d random points (incl. %d duplicates)\nand their Voronoi regions in %s' % (len(pts), N_DUPL, COUNTRY)) +ax.set_title('%d random points (incl. %d duplicates)\nand their Voronoi regions in %s' % (len(coords), N_DUPL, COUNTRY)) plt.tight_layout() plt.savefig('duplicate_points.png') diff --git a/geovoronoi/_voronoi.py b/geovoronoi/_voronoi.py index 54e2b00..6e6d816 100644 --- a/geovoronoi/_voronoi.py +++ b/geovoronoi/_voronoi.py @@ -133,46 +133,46 @@ def region_polygons_from_voronoi(vor, geom, return_point_assignments=False): else: p_vert_indices = set() p_vert_farpoints = set() - i_pt = pt_indices[0] # only consider one point, not a duplicate - 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]): - - 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) - 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 - finite_pt = vor.vertices[i] - p_vert_indices.add(i) - - # only add a far point if the finite end is not outside the bounding box already - if geom_bb.contains(asPoint(finite_pt)): - 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 - - isects = [] - for i_ext_coord in range(len(geom_bb.exterior.coords) - 1): - isect = line_segment_intersection(midpoint, direction, - np.array(geom_bb.exterior.coords[i_ext_coord]), - np.array(geom_bb.exterior.coords[i_ext_coord+1])) - if isect is not None: - isects.append(isect) - - if len(isects) == 0: - raise RuntimeError('far point must intersect with surrounding geometry from `geom`') - elif len(isects) == 1: - p_vert_farpoints.add(tuple(isects[0])) - else: - closest_isect_idx = np.argmin(np.linalg.norm(midpoint - isects, axis=1)) - p_vert_farpoints.add(tuple(isects[closest_isect_idx])) + for i_pt in pt_indices: # also consider duplicates + 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]): + + 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) + 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 + finite_pt = vor.vertices[i] + p_vert_indices.add(i) + + # only add a far point if the finite end is not outside the bounding box already + if geom_bb.contains(asPoint(finite_pt)): + 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 + + isects = [] + for i_ext_coord in range(len(geom_bb.exterior.coords) - 1): + isect = line_segment_intersection(midpoint, direction, + np.array(geom_bb.exterior.coords[i_ext_coord]), + np.array(geom_bb.exterior.coords[i_ext_coord+1])) + if isect is not None: + isects.append(isect) + + if len(isects) == 0: + raise RuntimeError('far point must intersect with surrounding geometry from `geom`') + elif len(isects) == 1: + p_vert_farpoints.add(tuple(isects[0])) + else: + closest_isect_idx = np.argmin(np.linalg.norm(midpoint - isects, axis=1)) + p_vert_farpoints.add(tuple(isects[closest_isect_idx])) # create the Voronoi region polygon as convex hull of the ridge vertices and far points (Voronoi regions # are convex)