-
Notifications
You must be signed in to change notification settings - Fork 27
/
duplicate_points.py
117 lines (82 loc) · 4.2 KB
/
duplicate_points.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
"""
Example script to show how to handle duplicate points for which Voronoi regions should be generated.
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
"""
import logging
import matplotlib.pyplot as plt
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.plotting import subplot_for_map, plot_voronoi_polys_with_points_in_area
logging.basicConfig(level=logging.INFO)
geovoronoi_log = logging.getLogger('geovoronoi')
geovoronoi_log.setLevel(logging.INFO)
geovoronoi_log.propagate = True
N_POINTS = 20
N_DUPL = 10
COUNTRY = 'Sweden'
np.random.seed(123)
print('loading country `%s` from naturalearth_lowres' % COUNTRY)
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
area = world[world.name == COUNTRY]
assert len(area) == 1
print('CRS:', area.crs) # gives epsg:4326 -> WGS 84
area = area.to_crs(epsg=3395) # convert to World Mercator CRS
area_shape = area.iloc[0].geometry # get the Polygon
# generate some random points within the bounds
minx, miny, maxx, maxy = area_shape.bounds
randx = np.random.uniform(minx, maxx, N_POINTS)
randy = np.random.uniform(miny, maxy, N_POINTS)
coords = np.vstack((randx, randy)).T
# 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
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
# duplicate a few points
rand_dupl_ind = np.random.randint(len(coords), size=N_DUPL)
coords = np.concatenate((coords, coords[rand_dupl_ind]))
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_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):
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):
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))
# 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
fig, ax = subplot_for_map()
plot_voronoi_polys_with_points_in_area(ax, area_shape, poly_shapes, coords,
plot_voronoi_opts={'alpha': 0.2},
plot_points_opts={'alpha': 0.4},
voronoi_color=list(vor_colors),
point_labels=pt_labels,
points_markersize=np.array(count_per_pt)*10)
ax.set_title('%d random points (incl. %d duplicates)\nand their Voronoi regions in %s' % (len(pts), N_DUPL, COUNTRY))
plt.tight_layout()
plt.savefig('duplicate_points.png')
plt.show()