Skip to content

Commit

Permalink
adapted duplicate_points example for new version
Browse files Browse the repository at this point in the history
  • Loading branch information
@internaut
internaut committed Jan 28, 2021
1 parent 7fc4169 commit cfe9e94
Show file tree
Hide file tree
Showing 3 changed files with 63 additions and 61 deletions.
Binary file modified examples/duplicate_points.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
44 changes: 23 additions & 21 deletions examples/duplicate_points.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
Duplicate points, i.e. points with exactly the same coordinates will belong to the same Voronoi region.
Author: Markus Konrad <markus.konrad@wzb.eu>
March 2018
January 2020
"""


Expand All @@ -23,6 +23,9 @@
geovoronoi_log.setLevel(logging.INFO)
geovoronoi_log.propagate = True


#%%

N_POINTS = 20
N_DUPL = 10
COUNTRY = 'Sweden'
Expand Down Expand Up @@ -61,56 +64,55 @@

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
#
# 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')
Expand Down
80 changes: 40 additions & 40 deletions geovoronoi/_voronoi.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down

0 comments on commit cfe9e94

Please sign in to comment.