diff --git a/examples/brandenburg.json b/examples/brandenburg.json new file mode 100644 index 0000000..b386131 --- /dev/null +++ b/examples/brandenburg.json @@ -0,0 +1,7 @@ +{ +"type": "FeatureCollection", +"name": "Germany_AL4", +"features": [ +{ "type": "Feature", "properties": { "srid": "4326", "id": "62504", "name": "Brandenburg", "localname": "Brandenburg", "official_name": "", "boundary": "adminstrative", "admin_level": "4", "note": "", "wikidata": "Q1208", "wikipedia": "de:Brandenburg", "timestamp": "2018\/07\/19 01:55:20", "rpath": "62504,51477,0", "alltags": { "ref": "BB", "name": "Brandenburg", "name:am": "ብራንደንቡርግ", "name:an": "Brandemburgo", "name:ar": "براندنبورغ", "name:ay": "Brandenburg suyu", "name:az": "Brandenburq", "name:ba": "Бранденбург", "name:be": "Брандэнбург", "name:bg": "Бранденбург", "name:bn": "ব্রান্ডেনবুর্গ", "name:ca": "Brandenburg", "name:ce": "Бранденбург", "name:cs": "Braniborsko", "name:cv": "Бранденбург", "name:de": "Brandenburg", "name:el": "Βρανδεμβούργο", "name:eo": "Brandenburgio", "name:es": "Brandemburgo", "name:eu": "Brandenburgo", "name:fa": "براندنبورگ", "name:fr": "Brandebourg", "name:fy": "Brandenburch", "name:gl": "Brandeburgo", "name:gn": "Brandeburgo", "name:he": "ברנדנבורג", "name:hi": "ब्रैंडेनबर्ग", "name:hy": "Բրանդենբուրգ", "name:ia": "Brandeburgo", "name:ie": "Brandenburgia", "name:is": "Brandenborg", "name:it": "Brandeburgo", "name:ja": "ブランデンブルク州", "name:ka": "ბრანდენბურგი", "name:kk": "Бранденбург", "name:ko": "브란덴부르크 주", "name:la": "Brandenburgum", "name:lb": "Brandenburg", "name:li": "Brandeburg", "name:lt": "Brandenburgas", "name:lv": "Brandenburga", "name:mk": "Бранденбург", "name:mn": "Бранденбург", "name:mr": "ब्रांडेनबुर्ग", "name:ne": "ब्रान्डेनबर्ग", "name:nl": "Brandenburg", "name:oc": "Brandeborg", "name:os": "Бранденбург", "name:pa": "ਬ੍ਰਾਂਡਨਬੁਰਕ", "name:pl": "Brandenburgia", "name:ps": "براندنبورگ", "name:pt": "Brandemburgo", "name:ru": "Бранденбург", "name:sk": "Brandenbursko", "name:sq": "Brandenburgu", "name:sr": "Бранденбург", "name:th": "รัฐบรันเดินบวร์ค", "name:tl": "Brandeburgo", "name:uk": "Бранденбург", "name:ur": "برندنبرگ", "name:vo": "Brandänburgän", "name:yi": "בראנדנבורג", "name:zh": "勃兰登堡", "boundary": "administrative", "name:hsb": "Braniborska", "wikidata": "Q1208", "ISO3166-2": "DE-BB", "wikipedia": "de:Brandenburg", "admin_level": "4", "alt_name:cs": "Země Braniborsko", "border_type": "state", "name:prefix": "Bundesland", "de:regionalschluessel": "12", "TMC:cid_58:tabcd_1:Class": "Area", "TMC:cid_58:tabcd_1:LCLversion": "12.0", "TMC:cid_58:tabcd_1:LocationCode": "267", "de:amtlicher_gemeindeschluessel": "12" } }, "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 11.267612, 53.121972 ], [ 11.393711, 53.10984 ], [ 11.395417, 53.140179 ], [ 11.42153, 53.141106 ], [ 11.53464, 53.125218 ], [ 11.569435, 53.165502 ], [ 11.550948, 53.20823 ], [ 11.629009, 53.241997 ], [ 11.699607, 53.244091 ], [ 11.728497, 53.21695 ], [ 11.827851, 53.227466 ], [ 11.795949, 53.253312 ], [ 11.861657, 53.248446 ], [ 11.892566, 53.278955 ], [ 11.948175, 53.271419 ], [ 11.973662, 53.297528 ], [ 12.018622, 53.299422 ], [ 12.014026, 53.33423 ], [ 12.055753, 53.349153 ], [ 12.045029, 53.370756 ], [ 12.061792, 53.371256 ], [ 12.079668, 53.366857 ], [ 12.081924, 53.346194 ], [ 12.106835, 53.34353 ], [ 12.139637, 53.361395 ], [ 12.168702, 53.339156 ], [ 12.229785, 53.358135 ], [ 12.260042, 53.323897 ], [ 12.310277, 53.327242 ], [ 12.393245, 53.302036 ], [ 12.392042, 53.284889 ], [ 12.428612, 53.274005 ], [ 12.445349, 53.249652 ], [ 12.529791, 53.264621 ], [ 12.605645, 53.244335 ], [ 12.660175, 53.254649 ], [ 12.674905, 53.248638 ], [ 12.67066, 53.229694 ], [ 12.757725, 53.224094 ], [ 12.739304, 53.199295 ], [ 12.764253, 53.188972 ], [ 12.848072, 53.200874 ], [ 12.879718, 53.179359 ], [ 12.965443, 53.200246 ], [ 12.979795, 53.19081 ], [ 12.942102, 53.174467 ], [ 12.980444, 53.162734 ], [ 13.078718, 53.200462 ], [ 13.10653, 53.214958 ], [ 13.104125, 53.236278 ], [ 13.149808, 53.25066 ], [ 13.180433, 53.251228 ], [ 13.203488, 53.22126 ], [ 13.229785, 53.216611 ], [ 13.252443, 53.262103 ], [ 13.298454, 53.28097 ], [ 13.360916, 53.277759 ], [ 13.386071, 53.244784 ], [ 13.405712, 53.246824 ], [ 13.44252, 53.286519 ], [ 13.43781, 53.299583 ], [ 13.48359, 53.291228 ], [ 13.502451, 53.327389 ], [ 13.523603, 53.319917 ], [ 13.516355, 53.352135 ], [ 13.556168, 53.376039 ], [ 13.550068, 53.399166 ], [ 13.624767, 53.408976 ], [ 13.64059, 53.447112 ], [ 13.711784, 53.48137 ], [ 13.812113, 53.484135 ], [ 13.826096, 53.498866 ], [ 13.778005, 53.513742 ], [ 13.791149, 53.558362 ], [ 13.82185, 53.518072 ], [ 13.87974, 53.504012 ], [ 13.875533, 53.473292 ], [ 13.918138, 53.45541 ], [ 13.902906, 53.43181 ], [ 13.918155, 53.421914 ], [ 14.001263, 53.434914 ], [ 14.081912, 53.410847 ], [ 14.121929, 53.441704 ], [ 14.179376, 53.422526 ], [ 14.224816, 53.434292 ], [ 14.242334, 53.421925 ], [ 14.230912, 53.367779 ], [ 14.099028, 53.26147 ], [ 14.264627, 53.259103 ], [ 14.261981, 53.276502 ], [ 14.302061, 53.286208 ], [ 14.316346, 53.312279 ], [ 14.402259, 53.330064 ], [ 14.415453, 53.324377 ], [ 14.406803, 53.308729 ], [ 14.421205, 53.276137 ], [ 14.449841, 53.260722 ], [ 14.404364, 53.209488 ], [ 14.376247, 53.199131 ], [ 14.366777, 53.168451 ], [ 14.38642, 53.140675 ], [ 14.349985, 53.057392 ], [ 14.257073, 53.001979 ], [ 14.143936, 52.961109 ], [ 14.160814, 52.886771 ], [ 14.123005, 52.837765 ], [ 14.216389, 52.816773 ], [ 14.28029, 52.774268 ], [ 14.350431, 52.751291 ], [ 14.46523, 52.661411 ], [ 14.596211, 52.610643 ], [ 14.639014, 52.573317 ], [ 14.603838, 52.531031 ], [ 14.633926, 52.491484 ], [ 14.548115, 52.432337 ], [ 14.534376, 52.395 ], [ 14.584829, 52.306347 ], [ 14.575288, 52.289048 ], [ 14.689475, 52.256758 ], [ 14.71555, 52.236128 ], [ 14.685814, 52.193906 ], [ 14.705567, 52.168914 ], [ 14.680016, 52.143383 ], [ 14.681812, 52.11564 ], [ 14.759094, 52.065383 ], [ 14.748422, 52.031928 ], [ 14.714029, 52.003683 ], [ 14.721659, 51.994475 ], [ 14.704855, 51.976008 ], [ 14.721287, 51.951315 ], [ 14.694184, 51.901895 ], [ 14.611105, 51.857131 ], [ 14.590107, 51.821362 ], [ 14.645475, 51.795528 ], [ 14.668163, 51.725857 ], [ 14.7474, 51.675669 ], [ 14.76512, 51.60747 ], [ 14.729824, 51.581589 ], [ 14.696636, 51.596762 ], [ 14.70554, 51.576932 ], [ 14.674128, 51.550852 ], [ 14.609559, 51.550107 ], [ 14.599006, 51.571187 ], [ 14.567686, 51.581087 ], [ 14.516481, 51.55428 ], [ 14.387123, 51.541567 ], [ 14.332329, 51.504899 ], [ 14.273278, 51.532376 ], [ 14.136602, 51.543333 ], [ 14.142298, 51.522966 ], [ 14.108512, 51.521879 ], [ 14.073806, 51.492216 ], [ 14.08894, 51.477779 ], [ 14.032744, 51.475003 ], [ 14.062587, 51.445421 ], [ 14.036877, 51.434537 ], [ 14.044861, 51.419073 ], [ 14.015904, 51.404489 ], [ 14.028122, 51.396264 ], [ 14.003784, 51.395536 ], [ 14.000596, 51.372537 ], [ 13.970463, 51.375699 ], [ 13.972375, 51.393434 ], [ 13.955527, 51.397329 ], [ 13.882903, 51.374429 ], [ 13.83426, 51.383752 ], [ 13.763328, 51.359064 ], [ 13.767724, 51.371696 ], [ 13.599149, 51.368673 ], [ 13.570333, 51.385659 ], [ 13.542538, 51.369405 ], [ 13.523314, 51.380931 ], [ 13.520726, 51.403668 ], [ 13.462982, 51.412214 ], [ 13.475845, 51.418519 ], [ 13.446034, 51.430853 ], [ 13.421893, 51.421568 ], [ 13.428913, 51.431774 ], [ 13.400668, 51.454248 ], [ 13.374004, 51.438229 ], [ 13.378692, 51.425301 ], [ 13.320841, 51.437384 ], [ 13.328584, 51.426473 ], [ 13.285623, 51.410778 ], [ 13.286452, 51.399297 ], [ 13.260211, 51.401434 ], [ 13.268686, 51.384585 ], [ 13.216339, 51.395353 ], [ 13.201536, 51.431875 ], [ 13.174032, 51.428282 ], [ 13.20328, 51.451605 ], [ 13.182092, 51.49191 ], [ 13.201199, 51.491789 ], [ 13.208042, 51.524272 ], [ 13.185566, 51.557986 ], [ 13.142584, 51.567988 ], [ 13.154377, 51.600253 ], [ 13.120197, 51.619986 ], [ 13.08593, 51.608095 ], [ 13.050548, 51.647505 ], [ 13.154312, 51.686403 ], [ 13.154189, 51.710474 ], [ 13.186882, 51.715693 ], [ 13.152096, 51.744237 ], [ 13.163733, 51.75484 ], [ 13.14744, 51.766448 ], [ 13.17055, 51.787018 ], [ 13.123927, 51.85429 ], [ 13.150508, 51.859678 ], [ 13.149145, 51.872145 ], [ 13.118628, 51.883107 ], [ 13.040988, 51.870573 ], [ 13.026893, 51.880785 ], [ 13.045574, 51.900435 ], [ 12.973468, 51.900767 ], [ 12.976643, 51.920968 ], [ 12.955935, 51.922641 ], [ 12.960558, 51.934642 ], [ 12.852415, 51.935165 ], [ 12.844325, 51.967582 ], [ 12.776118, 51.964935 ], [ 12.756244, 51.986499 ], [ 12.668292, 52.01306 ], [ 12.597264, 51.98152 ], [ 12.539432, 51.984922 ], [ 12.535818, 52.002928 ], [ 12.494043, 52.011615 ], [ 12.480728, 52.033175 ], [ 12.429478, 52.018333 ], [ 12.276527, 52.103891 ], [ 12.216113, 52.170543 ], [ 12.247333, 52.183879 ], [ 12.248787, 52.211782 ], [ 12.282094, 52.216715 ], [ 12.296648, 52.22846 ], [ 12.24541, 52.249863 ], [ 12.262923, 52.295045 ], [ 12.30809, 52.344632 ], [ 12.284122, 52.364199 ], [ 12.306555, 52.377714 ], [ 12.291711, 52.386111 ], [ 12.302269, 52.405309 ], [ 12.27461, 52.416183 ], [ 12.297198, 52.423694 ], [ 12.289413, 52.430579 ], [ 12.330787, 52.478274 ], [ 12.308522, 52.479194 ], [ 12.329612, 52.496642 ], [ 12.271845, 52.487501 ], [ 12.257958, 52.518197 ], [ 12.23607, 52.523873 ], [ 12.221768, 52.500231 ], [ 12.184673, 52.495722 ], [ 12.167122, 52.514696 ], [ 12.18874, 52.532226 ], [ 12.142172, 52.528285 ], [ 12.18152, 52.575218 ], [ 12.166939, 52.627464 ], [ 12.23646, 52.628728 ], [ 12.23259, 52.688766 ], [ 12.197363, 52.718902 ], [ 12.222162, 52.739826 ], [ 12.205848, 52.76049 ], [ 12.220406, 52.78989 ], [ 12.244407, 52.786606 ], [ 12.25762, 52.80487 ], [ 12.233428, 52.859922 ], [ 12.197054, 52.877878 ], [ 12.123719, 52.853635 ], [ 12.125539, 52.89439 ], [ 12.112775, 52.875695 ], [ 12.022564, 52.89038 ], [ 11.979908, 52.876112 ], [ 11.83449, 52.910145 ], [ 11.823288, 52.92125 ], [ 11.846216, 52.95132 ], [ 11.788134, 52.960208 ], [ 11.743077, 52.987857 ], [ 11.692056, 52.979161 ], [ 11.676296, 53.008221 ], [ 11.627038, 53.011963 ], [ 11.637715, 53.039368 ], [ 11.510955, 53.047305 ], [ 11.446347, 53.078476 ], [ 11.340537, 53.054726 ], [ 11.273544, 53.098469 ], [ 11.267612, 53.121972 ] ], [ [ 13.08835, 52.41964 ], [ 13.138103, 52.397856 ], [ 13.12722, 52.391666 ], [ 13.13108, 52.387236 ], [ 13.245947, 52.421126 ], [ 13.249725, 52.404957 ], [ 13.296755, 52.41625 ], [ 13.312082, 52.399139 ], [ 13.34328, 52.411273 ], [ 13.388437, 52.377865 ], [ 13.420834, 52.376127 ], [ 13.41878, 52.409995 ], [ 13.463563, 52.421065 ], [ 13.479767, 52.395947 ], [ 13.53843, 52.400603 ], [ 13.5355, 52.38899 ], [ 13.59275, 52.39383 ], [ 13.605798, 52.37361 ], [ 13.642682, 52.377509 ], [ 13.63632, 52.346819 ], [ 13.64843, 52.338276 ], [ 13.699975, 52.375599 ], [ 13.686815, 52.385306 ], [ 13.739034, 52.407335 ], [ 13.740057, 52.432333 ], [ 13.729895, 52.433924 ], [ 13.726043, 52.435298 ], [ 13.723017, 52.436721 ], [ 13.723066, 52.437057 ], [ 13.76116, 52.43771 ], [ 13.701267, 52.468224 ], [ 13.704616, 52.454756 ], [ 13.648568, 52.478744 ], [ 13.611369, 52.470653 ], [ 13.65695, 52.529883 ], [ 13.625691, 52.530182 ], [ 13.63739, 52.54225 ], [ 13.586412, 52.549782 ], [ 13.581446, 52.571081 ], [ 13.508126, 52.592159 ], [ 13.497131, 52.606708 ], [ 13.52287, 52.644812 ], [ 13.489883, 52.655517 ], [ 13.47951, 52.67551 ], [ 13.450794, 52.662706 ], [ 13.469909, 52.651915 ], [ 13.424448, 52.635545 ], [ 13.394554, 52.647554 ], [ 13.357644, 52.622977 ], [ 13.302609, 52.62719 ], [ 13.310177, 52.65737 ], [ 13.282752, 52.660761 ], [ 13.283772, 52.641113 ], [ 13.262141, 52.64072 ], [ 13.264234, 52.626865 ], [ 13.22057, 52.628242 ], [ 13.201643, 52.606479 ], [ 13.217397, 52.587477 ], [ 13.164528, 52.598793 ], [ 13.128966, 52.587302 ], [ 13.153538, 52.57324 ], [ 13.145686, 52.552791 ], [ 13.130503, 52.556023 ], [ 13.117393, 52.517033 ], [ 13.168632, 52.509253 ], [ 13.117693, 52.47732 ], [ 13.109278, 52.450632 ], [ 13.123153, 52.438709 ], [ 13.08835, 52.41964 ] ], [ [ 13.503438, 52.618993 ], [ 13.505442, 52.619687 ], [ 13.505449, 52.619884 ], [ 13.503814, 52.619223 ], [ 13.503438, 52.618993 ] ] ], [ [ [ 12.265677, 52.231337 ], [ 12.272055, 52.231249 ], [ 12.275011, 52.230527 ], [ 12.272864, 52.228934 ], [ 12.265677, 52.231337 ] ] ], [ [ [ 12.282741, 52.227261 ], [ 12.283196, 52.228181 ], [ 12.286261, 52.227816 ], [ 12.285645, 52.227026 ], [ 12.282741, 52.227261 ] ] ], [ [ [ 12.276428, 52.227809 ], [ 12.278846, 52.22944 ], [ 12.282289, 52.228014 ], [ 12.280859, 52.227598 ], [ 12.276428, 52.227809 ] ] ], [ [ [ 12.286549, 52.225852 ], [ 12.286796, 52.226286 ], [ 12.290414, 52.225202 ], [ 12.290055, 52.22474 ], [ 12.286549, 52.225852 ] ] ], [ [ [ 11.267355, 53.12452 ], [ 11.269086, 53.12301 ], [ 11.268632, 53.12211 ], [ 11.267954, 53.122929 ], [ 11.267355, 53.12452 ] ] ], [ [ [ 12.218076, 52.86133 ], [ 12.221282, 52.863301 ], [ 12.222102, 52.862805 ], [ 12.219041, 52.861117 ], [ 12.218076, 52.86133 ] ] ], [ [ [ 12.260021, 52.229894 ], [ 12.260818, 52.230771 ], [ 12.26174, 52.230569 ], [ 12.26097, 52.229695 ], [ 12.260021, 52.229894 ] ] ] ] } } +] +} diff --git a/examples/duplicate_points.png b/examples/duplicate_points.png index e62c854..7feecdb 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..d5bf5d7 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 2021 """ @@ -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 @@ -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,42 @@ # 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) +region_polys, region_pts = 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 region_pts.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 = points_to_region(region_pts) +for i_pt, i_poly in pts_to_poly_assignments.items(): print('> point ', i_pt, '-> voronoi region', i_poly) -# -# plotting -# +#%% 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 region_pts.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: (1,0,0) if len(pt_indices) > 1 else (0,0,1) + for i_poly, pt_indices in region_pts.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, region_polys, 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=(0,0,0,1), 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/examples/random_points_across_italy.png b/examples/random_points_across_italy.png index f653635..d19336e 100644 Binary files a/examples/random_points_across_italy.png and b/examples/random_points_across_italy.png differ diff --git a/examples/random_points_across_italy.py b/examples/random_points_across_italy.py index b098d45..ad8c389 100644 --- a/examples/random_points_across_italy.py +++ b/examples/random_points_across_italy.py @@ -1,28 +1,32 @@ """ Example script that scatters random points across a country and generates the Voronoi regions for them. Both the regions -and their points will be plotted using the `plotting` sub-module of `geovoronoi`. +and their points will be plotted using the `plotting` sub-module of `geovoronoi`. This example will show the effect +of calculating the Voronoi regions separately per country sub-geometry (e.g. separately for islands) vs. calculating +the Voronoi regions for the whole country shape together. Author: Markus Konrad -March 2018 +January 2021 """ import logging +from pprint import pprint 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 +from geovoronoi import coords_to_points, voronoi_regions_from_coords 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 = 100 +#%% + +N_POINTS = 105 COUNTRY = 'Italy' np.random.seed(123) @@ -46,29 +50,68 @@ # 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 +del coords # not used any more 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 +#%% # -# calculate the Voronoi regions, cut them with the geographic area shape and assign the points to them +# Calculate the Voronoi regions, cut them with the geographic area shape and assign the points to them. +# Note how in Sardinia there's only one point which is not assigned to any Voronoi region. +# Since by default all sub-geometries (i.e. islands or other isolated features) in the geographic shape +# are treated separately (set `per_geom=False` to change this) and you need more than one point to generate +# a Voronoi region, this point is left unassigned. By setting `return_unassigned_points=True`, we can get a +# set of unassigned point indices: # -poly_shapes, pts, poly_to_pt_assignments = voronoi_regions_from_coords(coords, area_shape) +region_polys, region_pts, unassigned_pts = voronoi_regions_from_coords(pts, area_shape, + return_unassigned_points=True, + per_geom=True) -# -# plotting -# +print('Voronoi region to point assignments:') +pprint(region_pts) + +print('Unassigned points:') +for i_pt in unassigned_pts: + print('#%d: %.2f, %.2f' % (i_pt, pts[i_pt].x, pts[i_pt].y)) + +#%% plotting fig, ax = subplot_for_map() -plot_voronoi_polys_with_points_in_area(ax, area_shape, poly_shapes, coords, poly_to_pt_assignments) -#plot_voronoi_polys_with_points_in_area(ax, area_shape, poly_shapes, coords) # monocolor +plot_voronoi_polys_with_points_in_area(ax, area_shape, region_polys, pts, region_pts, + point_labels=list(map(str, range(len(pts))))) ax.set_title('%d random points and their Voronoi regions in %s' % (len(pts), COUNTRY)) plt.tight_layout() plt.savefig('random_points_across_italy.png') plt.show() + + +#%% + +# Now we change `per_geom` to False. Since all geometries of the country are treated as one during Voronoi region +# generation, also the single point on Sardinia gets a Voronoi region. Note however, that these regions now can +# span over all geometries, e.g. regions from an island can span over to the mainland or another island and vice versa. +# You can see this effect for point #39 as its region spans from Sicilia to the "tiptoe" of Italy's mainland. + +region_polys2, region_pts2 = voronoi_regions_from_coords(pts, area_shape, per_geom=False) + +print('Voronoi region to point assignments:') +pprint(region_pts) + + +#%% plotting + +fig, ax = subplot_for_map() + +plot_voronoi_polys_with_points_in_area(ax, area_shape, region_polys2, pts, region_pts2, + point_labels=list(map(str, range(len(pts))))) + +ax.set_title('%d random points and their Voronoi regions in %s (with per_geom=False)' % (len(pts), COUNTRY)) + +plt.tight_layout() +plt.savefig('random_points_across_italy_per_geom_false.png') +plt.show() diff --git a/examples/random_points_across_italy_per_geom_false.png b/examples/random_points_across_italy_per_geom_false.png new file mode 100644 index 0000000..52ac13d Binary files /dev/null and b/examples/random_points_across_italy_per_geom_false.png differ diff --git a/examples/random_points_and_area.png b/examples/random_points_and_area.png index f1ad380..7321b49 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..b671362 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,39 +51,40 @@ # 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) +region_polys, region_pts = 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² +poly_areas = calculate_polygon_areas(region_polys, 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 -# +#%% plotting fig, ax = subplot_for_map(show_x_axis=True, show_y_axis=True) -voronoi_labels = ['%d km²' % round(a) for a in poly_areas] -plot_voronoi_polys_with_points_in_area(ax, area_shape, poly_shapes, coords, poly_to_pt_assignments, +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, region_polys, coords, region_pts, 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' % (n_pts, COUNTRY)) plt.tight_layout() plt.savefig('random_points_and_area.png') diff --git a/examples/random_points_brandenburg.png b/examples/random_points_brandenburg.png new file mode 100644 index 0000000..a80cfe5 Binary files /dev/null and b/examples/random_points_brandenburg.png differ diff --git a/examples/random_points_brandenburg.py b/examples/random_points_brandenburg.py new file mode 100644 index 0000000..c12f620 --- /dev/null +++ b/examples/random_points_brandenburg.py @@ -0,0 +1,85 @@ +""" +Example script that scatters random points across Brandenburg and generates the Voronoi regions for them. +The boundary shape of Brandenburg contains a hole (Berlin) and when loaded is regarded as "invalid" shape. + +Author: Markus Konrad +January 2021 +""" + + +import logging +from pprint import pprint + +import matplotlib.pyplot as plt +import numpy as np +import fiona +import pyproj +from shapely.geometry import shape +from shapely.ops import transform + +from geovoronoi import coords_to_points, voronoi_regions_from_coords +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 + +#%% + +INPUT_FILE = 'brandenburg.json' +N_POINTS = 20 + +np.random.seed(20210129) + +print('loading %s' % INPUT_FILE) + +with fiona.open(INPUT_FILE) as f: + src_crs = pyproj.CRS.from_dict(f.meta['crs']) + target_crs = pyproj.CRS.from_epsg(3395) # World Mercator CRS + print('source CRS:', src_crs) + print('target CRS:', target_crs) + crs_transformer = pyproj.Transformer.from_crs(src_crs, target_crs) + brandenburg = shape(f[0]['geometry']) + # note that we also apply ".buffer(0)", otherwise the shape is not valid + # see https://shapely.readthedocs.io/en/stable/manual.html#object.buffer on the "buffer(0)" trick + brandenburg = transform(crs_transformer.transform, brandenburg).buffer(0) + +#%% + +# generate some random points within the bounds +minx, miny, maxx, maxy = brandenburg.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(brandenburg)] # converts to shapely Point +del coords # not used any more + +print('will use %d of %d randomly generated points that are inside geographic area' % (len(pts), N_POINTS)) + +#%% + +# +# calculate the Voronoi regions, cut them with the geographic area shape and assign the points to them +# + +region_polys, region_pts = voronoi_regions_from_coords(pts, brandenburg) + +print('Voronoi region to point assignments:') +pprint(region_pts) + +#%% plotting + +fig, ax = subplot_for_map() + +plot_voronoi_polys_with_points_in_area(ax, brandenburg, region_polys, pts, region_pts, + point_labels=list(map(str, range(len(pts))))) + +ax.set_title('%d random points and their Voronoi regions in Brandenburg' % len(pts)) + +plt.tight_layout() +plt.savefig('random_points_brandenburg.png') +plt.show() diff --git a/examples/toyexample.py b/examples/toyexample.py new file mode 100644 index 0000000..bbfdbf0 --- /dev/null +++ b/examples/toyexample.py @@ -0,0 +1,54 @@ +""" +Artificial mini-example useful for debugging. + +Author: Markus Konrad +January 2021 +""" + + +import logging + +import matplotlib.pyplot as plt +import numpy as np + +from shapely.geometry import Polygon + +from geovoronoi import voronoi_regions_from_coords +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 + +#%% + +# points = np.array([[0, 0], [0, 1], [0, 2], +# [1, 0], [1, 1], [1, 2], +# [2, 0], [2, 1], [2, 2]]) + +points = np.array([[1, 1], [1.1, 0.9], [1.15, 0.8], [2.5, 0]]) +#points = np.array([[1, 1], [1.1, 0.91], [1.2, 0.8]]) + +# 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) + +plot_voronoi_polys_with_points_in_area(ax, shape, region_polys, points, region_pts, + point_labels=list(map(str, range(len(points)))), + voronoi_labels=list(map(str, region_polys.keys()))) + +ax.set_title('toy example') + +plt.tight_layout() +plt.show() diff --git a/examples/using_geopandas.png b/examples/using_geopandas.png index 1272452..8543aeb 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..609b73f 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,16 +47,17 @@ 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) +# we set "per_geom=False" so that the whole continent is treated as one area and Voronoi regions +# span over to Tierra del Fuego +region_polys, region_pts = voronoi_regions_from_coords(coords, south_am_shape, per_geom=False) -# -# Plotting -# +#%% Plotting 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, region_polys, coords, region_pts, + 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/geovoronoi/__init__.py b/geovoronoi/__init__.py index 7d76a66..b4d55c8 100644 --- a/geovoronoi/__init__.py +++ b/geovoronoi/__init__.py @@ -1,20 +1,19 @@ """ geovoronoi – main module -Imports all necessary functions to calculate Voronoi regions from a set of coordinates on a geographic shape. -Addtionally imports some helper funcitons. +Imports all necessary functions to calculate Voronoi regions from a set of coordinates inside a geographic shape. +Addtionally imports some helper functions. Author: Markus Konrad """ -from ._voronoi import (coords_to_points, points_to_coords, voronoi_regions_from_coords, polygon_lines_from_voronoi, - polygon_shapes_from_voronoi_lines, assign_points_to_voronoi_polygons, - get_points_to_poly_assignments) -from ._geom import calculate_polygon_areas +from ._voronoi import (coords_to_points, points_to_coords, voronoi_regions_from_coords, + points_to_region) +from ._geom import calculate_polygon_areas, line_segment_intersection __title__ = 'geovoronoi' -__version__ = '0.2.0' +__version__ = '0.3.0dev' __author__ = 'Markus Konrad' __license__ = 'Apache License 2.0' diff --git a/geovoronoi/_geom.py b/geovoronoi/_geom.py index fd10390..36b4ef3 100644 --- a/geovoronoi/_geom.py +++ b/geovoronoi/_geom.py @@ -1,85 +1,92 @@ """ Geometry helper functions in cartesian 2D space. -"shapely" refers to the [Shapely Python package for computational geometry](http://toblerity.org/shapely/index.html). +"Shapely" refers to the [Shapely Python package for computational geometry](http://toblerity.org/shapely/index.html). Author: Markus Konrad """ import numpy as np -from shapely.geometry import Polygon -def angle_between_pts(a, b, ref_vec=(1.0, 0.0)): +def line_segment_intersection(l_off, l_dir, segm_a, segm_b): """ - Angle *theta* between two points (numpy arrays) `a` and `b` in relation to a reference vector `ref_vec`. - By default, `ref_vec` is the x-axis (i.e. unit vector (1, 0)). - *theta* is in [0, 2Pi] unless `a` and `b` are very close. In this case this function returns NaN. + Check for line - segment intersection between line defined as `l_off + t * l_dir` and line segment between points + `segm_a` and `segm_b`. Hence the line is conceived as line starting at `l_off` and heading in direction `l_dir` + towards infinity. The function returns the intersection point as 1D NumPy array `[x, y]` if the given line hits + the defined segment between points `segm_a` and `segm_b` (endpoints inclusive). If there's no intersection, the + function returns None. + + All arguments must be 1D NumPy arrays of size 2. + + :param l_off: line offset point + :param l_dir: line direction vector + :param segm_a: segment start point + :param segm_b: segment end point + :return: intersection point as 1D NumPy array `[x, y]` or None if there is no intersection. """ - ang = inner_angle_between_vecs(a - b, np.array(ref_vec)) - if not np.isnan(ang) and a[1] < b[1]: - ang = 2 * np.pi - ang - return ang + if not all(isinstance(arg, np.ndarray) and arg.shape == (2,) for arg in (l_off, l_dir, segm_a, segm_b)): + raise ValueError('all arguments must be 1D NumPy arrays of size 2') + if np.isclose(np.linalg.norm(l_dir), 0): + raise ValueError('vector length of `l_dir` must be greater than 0') -def inner_angle_between_vecs(a, b): - """ - Return the inner angle *theta* between numpy vectors `a` and `b`. *theta* is in [0, Pi] if both `a` and `b` are - not at the origin (0, 0), otherwise this function returns NaN. - """ - origin = np.array((0, 0)) - if np.allclose(a, origin) or np.allclose(b, origin): - return np.nan - - au = a / np.linalg.norm(a) - bu = b / np.linalg.norm(b) - ang = np.arccos(np.clip(np.dot(au, bu), -1.0, 1.0)) + if np.isclose(np.linalg.norm(segm_b - segm_a), 0): + raise ValueError('vector length between `segm_a` and `segm_b` must be greater than 0') - return ang + segm_dir = segm_b - segm_a + v = np.array([l_dir, -segm_dir]).T + p = segm_a - l_off + if np.isclose(np.linalg.det(v), 0): # det of direction vector matrix is zero -> parallel direction vectors + if np.isclose(np.linalg.det(np.array([p, l_dir])), 0): # det of vector offset matrix is zero + # -> possible overlap + # order segment end points either horizontally or vertically + if segm_a[0] > segm_b[0] or \ + (np.isclose(segm_a[0], segm_b[0]) and segm_a[1] > segm_b[1]): # if horizontally aligned, + # order on y-axis + segm_b, segm_a = segm_a, segm_b -def polygon_around_center(points, center=None, fix_nan_angles=True): - """ - Order numpy array of coordinates `points` around `center` so that they form a valid polygon. Return that as - shapely `Polygon` object. If no valid polygon can be formed, return `None`. - If `center` is None (default), use midpoint of `points` as center. - """ - if center is None: - center = points.mean(axis=0) - else: - center = np.array(center) - - # angle between each point in `points` and `center` - angles = np.apply_along_axis(angle_between_pts, 1, points, b=center) + nonzero_ind = np.nonzero(l_dir)[0] # norm is > 0 so there must be a nonzero index + t_a = (segm_a - l_off)[nonzero_ind] / l_dir[nonzero_ind] # segm_a = l_off + t_a * l_dir + t_b = (segm_b - l_off)[nonzero_ind] / l_dir[nonzero_ind] # segm_b = l_off + t_b * l_dir - # sort by angles and generate polygon - if fix_nan_angles: - for repl in (0, np.pi): - tmp_angles = angles.copy() - tmp_angles[np.isnan(tmp_angles)] = repl - poly = Polygon(points[np.argsort(tmp_angles)]) - if poly.is_simple and poly.is_valid: - return poly + t = np.array([t_a, t_b]) + t = t[t >= 0] + if len(t) > 0: + return l_off + np.min(t) * l_dir - return None + return None # either parallel directions or line doesn't intersect with any segment endpoint else: - poly = Polygon(points[np.argsort(angles)]) + t = np.matmul(np.linalg.inv(v), p.T) + # intersection at l_off + t_0 * l_dir and segm_a + t_1 * segm_dir + # -> we're only interested if l_dir hits the segment (segm_a,segm_b) when it goes in positive direction, + # hence if t_0 is positive and t_1 is in [0, 1] - if poly.is_simple and poly.is_valid: - return poly + if t[0] >= 0 and 0 <= t[1] <= 1: + return segm_a + t[1] * segm_dir else: return None -def calculate_polygon_areas(poly_shapes, m2_to_km2=False): +def calculate_polygon_areas(region_polys, 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). + Return the area of the respective polygons in `poly_shapes` either in unit square meters (`m2_to_km2` is False) or + in square kilometers (`m2_to_km2` is True). Does not really calculate the area but uses the `area` property of + the Shapely polygons. + + Note: It is important to use an *equal area* projection with meters as units before using the areas of the + Voronoi regions! + + :param region_polys: dict mapping Voronoi region IDs to Shapely Polygon/MultiPolygon objects + :param m2_to_km2: if True, return results as square kilometers, otherwise square meters + :return: dict mapping Voronoi region IDs to area in square meters or square kilometers """ - areas = np.array([p.area for p in poly_shapes]) - if m2_to_km2: - return areas / 1000000 # = 1000² - else: - return areas + + if not isinstance(region_polys, dict): + raise ValueError('`poly_shapes` must be a dict') + + unit_convert = 1000000 if m2_to_km2 else 1 + return {i_poly: p.area / unit_convert for i_poly, p in region_polys.items()} diff --git a/geovoronoi/_voronoi.py b/geovoronoi/_voronoi.py index 5a98679..374632a 100644 --- a/geovoronoi/_voronoi.py +++ b/geovoronoi/_voronoi.py @@ -1,7 +1,7 @@ """ Functions to create Voronoi regions from points inside a geographic area. -"shapely" refers to the [Shapely Python package for computational geometry](http://toblerity.org/shapely/index.html). +"Shapely" refers to the [Shapely Python package for computational geometry](http://toblerity.org/shapely/index.html). Author: Markus Konrad """ @@ -11,12 +11,12 @@ import numpy as np from scipy.spatial import Voronoi -from scipy.spatial.distance import cdist -from shapely.geometry import LineString, asPoint, MultiPoint, Polygon +from scipy.spatial.qhull import QhullError +from shapely.geometry import box, LineString, asPoint, MultiPoint, Polygon, MultiPolygon from shapely.errors import TopologicalError -from shapely.ops import polygonize, cascaded_union +from shapely.ops import cascaded_union -from ._geom import polygon_around_center +from ._geom import line_segment_intersection logger = logging.getLogger('geovoronoi') @@ -24,295 +24,469 @@ def coords_to_points(coords): - """Convert a NumPy array of 2D coordinates `coords` to a list of shapely Point objects""" + """ + Convert a NumPy array of 2D coordinates `coords` to a list of Shapely Point objects. + + This is the inverse of `points_to_coords()`. + + :param coords: NumPy array of shape (N,2) with N coordinates in 2D space + :return: list of length N with Shapely Point objects + """ return list(map(asPoint, coords)) def points_to_coords(pts): - """Convert a list of shapely Point objects to a NumPy array of 2D coordinates `coords`""" + """ + Convert a list of Shapely Point objects to a NumPy array of 2D coordinates `coords`. + + This is the inverse of `coords_to_points()`. + + :param pts: list of length N with Shapely Point objects + :return: NumPy array of shape (N,2) with N coordinates in 2D space + """ return np.array([p.coords[0] for p in pts]) -def voronoi_regions_from_coords(coords, geo_shape, - shapes_from_diff_with_min_area=None, - accept_n_coord_duplicates=None, - return_unassigned_points=False, - farpoints_max_extend_factor=10): +def voronoi_regions_from_coords(coords, geo_shape, per_geom=True, return_unassigned_points=False, + results_per_geom=False, **kwargs): """ - Calculate Voronoi regions from NumPy array of 2D coordinates `coord` that lie within a shape `geo_shape`. Setting - `shapes_from_diff_with_min_area` fixes rare errors where the Voronoi shapes do not fully cover `geo_shape`. Set this - to a small number that indicates the minimum valid area of "fill up" Voronoi region shapes. - Set `accept_n_coord_duplicates` to accept exactly this number of points with exactly the same coordinates. Such - duplicates will belong to the same Voronoi region. If set to `None` then any number of duplicate coordinates is - accepted. Set `return_unassigned_points` to True to additionally return a list of shapely Point objects that could - not be assigned to any Voronoi region. - - This function returns the following values in a tuple: - - 1. `poly_shapes`: a list of shapely Polygon/MultiPolygon objects that represent the generated Voronoi regions - 2. `points`: the input coordinates converted to a list of shapely Point objects - 3. `poly_to_pt_assignments`: a nested list that for each Voronoi region in `poly_shapes` contains a list of indices - into `points` (or `coords`) that represent the points that belong to this Voronoi region. Usually, this is only - a single point. However, in case of duplicate points (e.g. both or more points have exactly the same coordinates) - then all these duplicate points are listed for the respective Voronoi region. - 4. optional if `return_unassigned_points` is True: a list of points that could not be assigned to any Voronoi region - - When calculating the far points of loose ridges for the Voronoi regions, `farpoints_max_extend_factor` is the - factor that is multiplied with the maximum extend per dimension. Increase this number in case the hull of far points - doesn't intersect with `geo_shape`. + Generate Voronoi regions from NumPy array of 2D coordinates or list of Shapely Point objects in `coord`. These + points must lie within a shape `geo_shape` which must be a valid Shapely Polygon or MultiPolygon object. If + `geo_shape` is a MultiPolygon, each of its sub-geometries will be either treated separately during Voronoi + region generation when `per_geom` is True or otherwise the whole MultiPolygon is treated as one object. In the + former case, Voronoi regions may not span from one sub-geometry to another (e.g. from one island to another + island) which also means that sub-geometries may remain empty (e.g. when there are no points on an island). In the + latter case Voronoi regions from one sub-geometry may span to another sub-geometry, hence all sub-geometries + should be covered by a Voronoi region as a result. + + This function returns at least two dicts in a tuple, called `region_polys` and `region_pts`. The first contains a + dict that maps unique Voronoi region IDs to the generated Voronoi region geometries as Shapely Polygon/MultiPolygon + objects. The second contains a dict that maps the region IDs to a list of point indices of `coords`. This dict + describes which Voronoi region contains which points. By definition a Voronoi region surrounds only a single point. + However, if you have duplicate coordinates in `coords`, these duplicate points will be surrounded by the same + Voronoi region. + + The structure of the returned dicts depends on `results_per_geom`. If `results_per_geom` is False, there is a direct + mapping from Voronoi region ID to region geometry or assigned points respectively. If `results_per_geom` is True, + both dicts map a sub-geometry ID from `geo_shape` (the index of the sub-geometry in `geo_shape.geoms`) to the + respective dict which in turn map a Voronoi region ID to its region geometry or assigned points. + + To conclude with N coordinates in `coords` and `results_per_geom` is False (default): + + ``` + region_polys = + { + 0: , + 1: , + ... + N-1: + } + region_pts = + { + 0: [5], + 1: [7, 2], # coords[7] and coords[2] are duplicates + ... + N-1: [3] + } + ``` + + And if `results_per_geom` is True and there are M sub-geometries in `geo_shape`: + + ``` + region_polys = + { + 0: { # first sub-geometry in `geom_shape` with Voronoi regions inside + 4: , + 2: , + ... + }, + ... + M-1: { # last sub-geometry in `geom_shape` with Voronoi regions inside + 9: , + 12: , + ... + }, + } + + region_pts = (similar to above) + ``` + + Setting `results_per_geom` to True only makes sense when `per_geom` is True so that Voronoi region generation is + done separately for each sub-geometry in `geo_shape`. + + :param coords: NumPy array of 2D coordinates as shape (N,2) for N points or list of Shapely Point objects; should + contain at least two points + :param geo_shape: Shapely Polygon or MultiPolygon object that defines the restricting area of the Voronoi regions; + all points in `coords` should be within `geo_shape`; make sure that `geo_shape` is a valid + geometry + :param per_geom: if True, treat sub-geometries in `geo_shape` separately during Voronoi region generation + :param return_unassigned_points: If True, additionally return set of point indices which could not be assigned to + any Voronoi region (usually because there were too few points inside a sub-geometry + to generate Voronoi regions) + :param results_per_geom: partition the result dicts by sub-geometry index + :param kwargs: parameters passed to `region_polygons_from_voronoi()` + :return tuple of length two if return_unassigned_points is False, otherwise of length three with: (1) dict + containing generated Voronoi region geometries as Shapely Polygon/MultiPolygon, (2) dict mapping Voronoi + regions to point indices of `coords`, (3 - optionally) set of point indices which could not be assigned to + any Voronoi region """ - logger.info('running Voronoi tesselation for %d points' % len(coords)) - vor = Voronoi(coords) - logger.info('generated %d Voronoi regions' % (len(vor.regions)-1)) + logger.info('running Voronoi tesselation for %d points / treating geoms separately: %s' % (len(coords), per_geom)) - logger.info('generating Voronoi polygon lines') - poly_lines = polygon_lines_from_voronoi(vor, geo_shape, farpoints_max_extend_factor=farpoints_max_extend_factor) + # check input arguments - logger.info('generating Voronoi polygon shapes') - poly_shapes = polygon_shapes_from_voronoi_lines(poly_lines, geo_shape, - shapes_from_diff_with_min_area=shapes_from_diff_with_min_area) + if len(coords) < 2: + raise ValueError('insufficient number of points provided in `coords`') - logger.info('assigning %d points to %d Voronoi polygons' % (len(coords), len(poly_shapes))) - points = coords_to_points(coords) - poly_to_pt_assignments, unassigned_pts = assign_points_to_voronoi_polygons(points, poly_shapes, - accept_n_coord_duplicates=accept_n_coord_duplicates, - return_unassigned_points=True, - coords=coords) - - if return_unassigned_points: - return poly_shapes, points, poly_to_pt_assignments, list(unassigned_pts) + # 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: - return poly_shapes, points, poly_to_pt_assignments - - -def polygon_lines_from_voronoi(vor, geo_shape, return_only_poly_lines=True, farpoints_max_extend_factor=10): - """ - Takes a scipy Voronoi result object `vor` (see [1]) and a shapely Polygon `geo_shape` the represents the geographic - area in which the Voronoi regions shall be placed. Calculates the following three lists: - - 1. Polygon lines of the Voronoi regions. These can be used to generate all Voronoi region polygons via - `polygon_shapes_from_voronoi_lines`. - 2. Loose ridges of Voronoi regions. - 3. Far points of loose ridges of Voronoi regions. + pts = coords + coords = points_to_coords(pts) - If `return_only_poly_lines` is True, only the first list is returned, otherwise a tuple of all three lists is - returned. + if not isinstance(geo_shape, (Polygon, MultiPolygon)): + raise ValueError('`geo_shape` must be a Polygon or MultiPolygon') - When calculating the far points of loose ridges, `farpoints_max_extend_factor` is the factor that is multiplied - with the maximum extend per dimension. Increase this number in case the hull of far points doesn't intersect - with `geo_shape`. + if not geo_shape.is_valid: + raise ValueError('`geo_shape` is not a valid shape; try applying `.buffer(b)` where `b` is zero ' + 'or a very small number') - Calculation of Voronoi region polygon lines taken and adopted from [2]. - - [1]: https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.Voronoi.html#scipy.spatial.Voronoi - [2]: https://github.com/scipy/scipy/blob/v1.0.0/scipy/spatial/_plotutils.py - """ + # 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 has sub-geometries + geoms = geo_shape.geoms - max_dim_extend = vor.points.ptp(axis=0).max() * farpoints_max_extend_factor - center = np.array(MultiPoint(vor.points).convex_hull.centroid) + # set up data containers + pts_indices = set(range(len(pts))) # point indices to operate on + geom_region_polys = {} + geom_region_pts = {} + unassigned_pts = set() - # 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) - if np.all(simplex >= 0): # full finite polygon line - poly_lines.append(LineString(vor.vertices[simplex])) - else: # "loose ridge": contains infinite Voronoi vertex - 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 - - 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]))) - - if return_only_poly_lines: - return poly_lines - else: - return poly_lines, loose_ridges, far_points + # 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] = {} -def polygon_shapes_from_voronoi_lines(poly_lines, geo_shape=None, shapes_from_diff_with_min_area=None): - """ - Form shapely Polygons objects from a list of shapely LineString objects in `poly_lines` by using - [`polygonize`](http://toblerity.org/shapely/manual.html#shapely.ops.polygonize). If `geo_shape` is not None, then - the intersection between any generated polygon and `geo_shape` is taken in case they overlap (i.e. the Voronoi - regions at the border are "cut" to the `geo_shape` polygon that represents the geographic area holding the - Voronoi regions). Setting `shapes_from_diff_with_min_area` fixes rare errors where the Voronoi shapes do not fully - cover `geo_shape`. Set this to a small number that indicates the minimum valid area of "fill up" Voronoi region - shapes. - Returns a list of shapely Polygons objects. - """ + logger.info('%d of %d points in geometry #%d of %d' % (len(pts_in_geom), len(pts), i_geom + 1, len(geoms))) - # generate shapely Polygon objects from the LineStrings of the Voronoi shapes in `poly_lines` - 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 pts_in_geom: # no points within `geom` -> skip + continue - if not p.is_empty: - poly_shapes.append(p) + # remove the points that we're about to use (point - geometry assignment is bijective) + pts_indices.difference_update(pts_in_geom) - if geo_shape is not None and shapes_from_diff_with_min_area is not None: - # fix rare cases where the generated polygons of the Voronoi regions don't fully cover `geo_shape` - vor_polys_union = cascaded_union(poly_shapes) # union of Voronoi regions + # generate the Voronoi regions using the SciPy Voronoi class + logger.info('generating Voronoi regions') try: - diff = np.array(geo_shape.difference(vor_polys_union), dtype=object) # "gaps" - except TopologicalError: - diff = np.array([], dtype=object) + vor = Voronoi(coords[pts_in_geom]) + 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 # 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) + + # 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: # nothing to do here since we already have results per sub-geometry + region_polys = geom_region_polys + region_pts = geom_region_pts + 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()))) + + region_pts.update({reg_id: geom_pts[old_reg_id] for reg_id, old_reg_id in region_ids_mapping.items()}) - if diff.shape and len(diff) > 0: - diff_areas = np.array([p.area for p in diff]) # areas of "gaps" - # use only those "gaps" bigger than `shapes_from_diff_with_min_area` because very tiny areas are generated - # at the borders due to floating point errors - poly_shapes.extend(diff[diff_areas >= shapes_from_diff_with_min_area]) - - return poly_shapes + if return_unassigned_points: + return region_polys, region_pts, unassigned_pts + else: + return region_polys, region_pts -def assign_points_to_voronoi_polygons(points, poly_shapes, - accept_n_coord_duplicates=None, - return_unassigned_points=False, - coords=None): +def region_polygons_from_voronoi(vor, geom, return_point_assignments=False, + bounds_buf_factor=0.1, + diff_topo_error_buf_factor=0.00000001): """ - Assign a list/array of shapely Point objects `points` to their respective Voronoi polygons passed as list - `poly_shapes`. Return a list of `assignments` of size `len(poly_shapes)` where ith element in `assignments` - contains the index of the point in `points` that resides in the ith Voronoi region. - Normally, 1-to-1 assignments are expected, i.e. for each Voronoi region in `poly_shapes` there's exactly one point - in `points` belonging to this Voronoi region. However, if there are duplicate coordinates in `points`, then those - duplicates will be assigned together to their Voronoi region and hence there is a 1-to-n relationship between - Voronoi regions and points. If `accept_n_coord_duplicates` is set to None, then an an unspecified number of - duplicates are allowed. If `accept_n_coord_duplicates` is 0, then no point duplicates are allowed, otherwise - exactly `accept_n_coord_duplicates` duplicates are allowed. - Set `return_unassigned_points` to additionally return a list of points that could not be assigned to any Voronoi - region. `coords` can be passed in order to avoid another conversion from Point objects to NumPy coordinate array. + 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 """ - # do some checks first + # check input arguments - n_polys = len(poly_shapes) - n_points = len(points) + if not isinstance(vor, Voronoi): + raise ValueError('`vor` must be SciPy Voronoi object') - if n_polys > n_points: - raise ValueError('The number of voronoi regions must be smaller or equal to the number of points') + if not isinstance(geom, (Polygon, MultiPolygon)): + raise ValueError('`geom` must be a Polygon or MultiPolygon') - if accept_n_coord_duplicates is None: - dupl_restricted = False - accept_n_coord_duplicates = n_points - n_polys - else: - dupl_restricted = True - - expected_n_geocoords = n_points - accept_n_coord_duplicates - - if coords is None: - coords = points_to_coords(points) - elif len(coords) != n_points: - raise ValueError('`coords` and `points` must have the same length') - - if accept_n_coord_duplicates >= 0 and n_polys != expected_n_geocoords: - raise ValueError('Unexpected number of geo-coordinates: %d (got %d polygons and expected %d geo-coordinates)' % - (n_points, n_polys, expected_n_geocoords)) - - # get the Voronoi regions' centroids and calculate the distance between all centroid – input coordinate pairs - poly_centroids = np.array([p.centroid.coords[0] for p in poly_shapes]) - poly_pt_dists = cdist(poly_centroids, coords) - - # generate the assignments - assignments = [] - already_assigned = set() - n_assigned_dupl = 0 - for i_poly, vor_poly in enumerate(poly_shapes): - # indices to points sorted by distance to this Voronoi region's centroid - closest_pt_indices = np.argsort(poly_pt_dists[i_poly]) - assigned_pts = [] - n_assigned = len(assigned_pts) - for i_pt in closest_pt_indices: # check each point with increasing distance - pt = points[i_pt] - - if pt.intersects(vor_poly): # if this point is inside this Voronoi region, assign it to the region - if i_pt in already_assigned: - raise RuntimeError('Point %d cannot be assigned to more than one voronoi region' % i_pt) - assigned_pts.append(i_pt) - already_assigned.add(i_pt) - if n_assigned >= accept_n_coord_duplicates - n_assigned_dupl: - break - - if not assigned_pts: - raise RuntimeError('Polygon %d does not contain any point' % i_poly) - - if accept_n_coord_duplicates == 0 and len(assigned_pts) > 1: - raise RuntimeError('No duplicate points allowed but polygon %d contains several points: %s' - % (i_poly, str(assigned_pts))) - - # add the assignments for this Voronoi region - assignments.append(assigned_pts) - n_assigned_dupl += len(assigned_pts)-1 - - # make some final checks - - assert len(assignments) == len(poly_shapes) - - if dupl_restricted: - assert n_assigned_dupl == accept_n_coord_duplicates - assert sum(map(len, assignments)) == len(poly_shapes) + accept_n_coord_duplicates - assert len(already_assigned) == n_points # make sure all points were assigned - unassigned_pt_indices = set() - else: - unassigned_pt_indices = set(range(n_points)) - already_assigned + # 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) - if return_unassigned_points: - return assignments, unassigned_pt_indices + # center of all points + center = np.array(MultiPoint(vor.points).convex_hull.centroid) + + # ridge vertice's coordinates + ridge_vert = np.asarray(vor.ridge_vertices) + + # 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 + + 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) # 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: + # 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) + + # 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(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 the bounding box + i = vertidx[vertidx >= 0][0] # finite end Voronoi vertex + finite_pt = vor.vertices[i] + p_vert_indices.add(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) # 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, + 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('ridge line must intersect with surrounding geometry from `geom`') + elif len(isects) == 1: # one intersection + far_pt = isects[0] + 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: + p_vert_farpoints.add(tuple(far_pt)) + + # create the Voronoi region polygon as convex hull of the ridge vertices and far points (Voronoi regions + # are convex) + + p_pts = vor.vertices[np.asarray(list(p_vert_indices))] + # adding the point itself prevents generating invalid hull (otherwise sometimes the hull becomes a line + # if the vertices are almost colinear) + p_pts = np.vstack((p_pts, vor.points[pt_indices])) + if p_vert_farpoints: + p_pts = np.vstack((p_pts, np.asarray(list(p_vert_farpoints)))) + + 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) + + if not p.is_valid or p.is_empty: + 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 + + # --- 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) + except StopIteration: + if pass_ == 1: # restart w/ second pass + polys_iter = iter(region_polys.items()) + i_reg, p = next(polys_iter) + pass_ += 1 + 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 + try: + diff = geom.buffer(max_extend * diff_topo_error_buf_factor).difference(union_other_regions) + except TopologicalError: + raise RuntimeError('difference operation failed with TopologicalError; try setting a different ' + '`diff_topo_error_buf_factor`') + + if isinstance(diff, MultiPolygon): + diff = diff.geoms + else: + diff = [diff] + + # list of neighbor region IDs + neighbor_regions = [vor.point_region[pt] for pt in region_neighbor_pts[i_reg]] + 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) + if pass_ == 1: + covered_area += (diff_part.area - p.area) + 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)) + + if return_point_assignments: + return region_polys, region_pts else: - return assignments + return region_polys -def get_points_to_poly_assignments(poly_to_pt_assignments): +def points_to_region(region_pts): """ - Reverse of poly to points assignments: Returns a list of size N, which is the number of unique points in - `poly_to_pt_assignments`. Each list element is an index into the list of Voronoi regions. + 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 region_pts: dict mapping Voronoi region polygon IDs to list of point IDs + :return: dict mapping point ID to Voronoi region polygon ID """ - pt_poly = [(i_pt, i_vor) - for i_vor, pt_indices in enumerate(poly_to_pt_assignments) - for i_pt in pt_indices] - return [i_vor for _, i_vor in sorted(pt_poly, key=lambda x: x[0])] + pts_region = {} + for poly, pts in region_pts.items(): + for pt in pts: + 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 pts_region diff --git a/geovoronoi/plotting.py b/geovoronoi/plotting.py index fe34e07..d00c635 100644 --- a/geovoronoi/plotting.py +++ b/geovoronoi/plotting.py @@ -1,24 +1,29 @@ """ 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 """ import numpy as np import matplotlib.pyplot as plt +from matplotlib.collections import PatchCollection +from descartes.patch import PolygonPatch 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) @@ -36,84 +41,186 @@ 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, 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 = generate_n_colors(len(poly_shapes), cmap_name=cmap_name) + vor_colors = {p_id: col + for p_id, col in zip(region_polys.keys(), generate_n_colors(len(region_polys), cmap_name=cmap_name))} + + pt_to_poly = points_to_region(region_pts) - pt_colors = [vor_colors[i_vor] for i_vor in get_points_to_poly_assignments(poly_to_pt_assignments)] + 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 = [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) return vor_colors, pt_colors -def plot_voronoi_polys(ax, poly_shapes, color=None, edgecolor=None, labels=None, label_fontsize=10, label_color=None, +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: + if not isinstance(points, np.ndarray): + coords = points_to_coords(points) + else: + coords = points + return coords[:, 0], coords[:, 1] + + +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, lbl) in enumerate(zip(poly_shapes, labels)): + for i, p in region_polys.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, marker='o', color=None, labels=None, label_fontsize=7, label_color=None, +def plot_points(ax, points, markersize=1, marker='o', color=None, labels=None, label_fontsize=7, label_color=None, label_draw_duplicates=False, **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`. + + :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 """ - if not isinstance(points, np.ndarray): - coords = points_to_coords(points) - else: - coords = points + x, y = xy_from_points(points) - ax.scatter(coords[:, 0], coords[:, 1], s=markersize, marker=marker, color=color, **kwargs) + ax.scatter(x, y, s=markersize, marker=marker, color=color, **kwargs) if labels: # plot labels using matplotlib's text() n_labels = len(labels) - n_features = len(coords) + n_features = len(x) if n_labels != n_features: raise ValueError('number of labels (%d) must match number of points (%d)' % (n_labels, n_features)) drawn_coords = set() - for i, (pos, lbl) in enumerate(zip(coords, labels)): - pos = tuple(pos) # make hashable + 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(pos[0], pos[1], 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) -def plot_voronoi_polys_with_points_in_area(ax, area_shape, poly_shapes, points, poly_to_pt_assignments=None, - area_color='white', area_edgecolor='black', +def plot_line(ax, points, linewidth=1, color=None, **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) + + ax.plot(x, y, linewidth=linewidth, color=color, **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, + **kwargs)) + if label: + 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, 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, points_color=None, points_markersize=5, points_marker='o', @@ -123,29 +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, + 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 = 'black' # better visible default value + 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) @@ -154,66 +288,71 @@ 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 'black' - - -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` + 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*. - Parameters - ---------- + Plots a collection of Polygon and MultiPolygon geometries to `ax`. - ax : matplotlib.axes.Axes - where shapes will be plotted - - 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 + :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 + """ - **kwargs - Additional keyword arguments passed to the collection + color_values = color + color_indices = None - Returns - ------- + if not isinstance(geoms, GeoSeries): + if isinstance(geoms, dict): + geoms_indices = np.array(list(geoms.keys())) + geoms_values = list(geoms.values()) - collection : matplotlib.collections.Collection that was plotted - """ - from descartes.patch import PolygonPatch - from matplotlib.collections import PatchCollection + if isinstance(color, dict): + color_indices = np.array(list(color.keys())) + color_values = np.array(list(color.values())) + else: + geoms_indices = np.arange(len(geoms)) + geoms_values = geoms + geoms = GeoSeries(geoms_values) + else: + geoms_indices = geoms.index.to_numpy() - if not isinstance(geoms, GeoSeries): - geoms = GeoSeries(geoms) + if isinstance(color, list): + color_values = np.array(color) + color_indices = np.arange(len(color)) + elif isinstance(color, np.ndarray): + color_values = color + color_indices = np.arange(len(color)) - geoms, indices = _flatten_multi_geoms(geoms) + geoms, multi_indices = _flatten_multi_geoms(geoms) - if isinstance(color, (list, tuple)): - color = np.array(color)[indices] + 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. if 'markersize' in kwargs: del kwargs['markersize'] collection = PatchCollection([PolygonPatch(poly) for poly in geoms], - color=color, **kwargs) + color=color_values, **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/requirements.txt b/requirements.txt index 9b24f88..1fb7e26 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,52 +1,6 @@ -appdirs==1.4.3 -attrs==19.3.0 -backcall==0.1.0 -Click==7.0 -click-plugins==1.1.1 -cligj==0.5.0 -cycler==0.10.0 -decorator==4.4.2 -descartes==1.1.0 -distlib==0.3.0 -filelock==3.0.12 -Fiona==1.8.13.post1 -geopandas==0.7.0 -hypothesis==5.6.0 -importlib-metadata==1.5.0 -importlib-resources==1.2.0 -ipython==7.13.0 -ipython-genutils==0.2.0 -jedi==0.16.0 -kiwisolver==1.1.0 -matplotlib==3.2.0 -more-itertools==8.2.0 -munch==2.5.0 -nose==1.3.7 -numpy==1.18.1 -packaging==20.3 -pandas==1.0.1 -parso==0.6.2 -pexpect==4.8.0 -pickleshare==0.7.5 -Pillow==7.0.0 -pluggy==0.13.1 -prompt-toolkit==3.0.3 -ptyprocess==0.6.0 -py==1.8.1 -Pygments==2.5.2 -pyparsing==2.4.6 -pyproj==2.5.0 -pytest==5.3.5 -pytest-mpl==0.11 -python-dateutil==2.8.1 -pytz==2019.3 -scipy==1.4.1 -Shapely==1.7.0 -six==1.14.0 -sortedcontainers==2.1.0 -toml==0.10.0 -tox==3.14.5 -traitlets==4.3.3 -virtualenv==20.0.8 -wcwidth==0.1.8 -zipp==3.1.0 +# requirements.txt +# +# installs dependencies from ./setup.py, and the package itself, +# in editable mode for development + +-e .[all] diff --git a/setup.py b/setup.py index 6df8b8d..d62d27c 100644 --- a/setup.py +++ b/setup.py @@ -8,12 +8,25 @@ GITHUB_URL = 'https://github.com/WZBSocialScienceCenter/geovoronoi' __title__ = 'geovoronoi' -__version__ = '0.2.0' +__version__ = '0.3.0dev' __author__ = 'Markus Konrad' __license__ = 'Apache License 2.0' here = os.path.abspath(os.path.dirname(__file__)) +DEPS_BASE = ['numpy>=1.19.0,<2', 'scipy>=1.5.0,<1.7', 'shapely>=1.7.0,<1.8'] + +DEPS_EXTRA = { + 'plotting': ['matplotlib>=3.3.0,<3.4', 'geopandas>=0.8.0,<0.9', 'descartes>=1.1.0,<1.2'], + 'test': ['pytest>=6.2.0,<6.3', 'pytest-mpl>=0.12,<0.13', 'hypothesis>=6.0.0,<6.1', 'tox>=3.21.0,<3.22'], + 'develop': ['ipython>=7.19.0', 'twine>=3.3.0'], +} + +DEPS_EXTRA['all'] = [] +for k, deps in DEPS_EXTRA.items(): + if k != 'all': + DEPS_EXTRA['all'].extend(deps) + # Get the long description from the README file with open(os.path.join(here, 'README.md'), encoding='utf-8') as f: long_description = f.read() @@ -43,10 +56,10 @@ 'Operating System :: OS Independent', 'Programming Language :: Python', 'Programming Language :: Python :: 3', - 'Programming Language :: Python :: 3.4', - 'Programming Language :: Python :: 3.5', 'Programming Language :: Python :: 3.6', 'Programming Language :: Python :: 3.7', + 'Programming Language :: Python :: 3.8', + 'Programming Language :: Python :: 3.9', 'Topic :: Scientific/Engineering :: Information Analysis', 'Topic :: Software Development :: Libraries :: Python Modules', @@ -58,10 +71,6 @@ packages=['geovoronoi'], # include_package_data=True, python_requires='>=3.6', - install_requires=['numpy>=1.18.0', 'scipy>=1.4.0', 'shapely>=1.7.0'], - extras_require={ - 'plotting': ['matplotlib>=3.2.0', 'geopandas>=0.7.0', 'descartes>=1.1.0'], - 'test': ['pytest>=5.3.0', 'pytest-mpl>=0.11', 'hypothesis>=5.6.0', 'tox>=3.14.0'], - 'develop': ['ipython>=7.0.0', 'twine>=3.1.0'], - } + install_requires=DEPS_BASE, + extras_require=DEPS_EXTRA ) 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/baseline/test_issue_7b.png b/tests/baseline/test_issue_7b.png index 5a7ad20..3080177 100644 Binary files a/tests/baseline/test_issue_7b.png and b/tests/baseline/test_issue_7b.png differ diff --git a/tests/baseline/test_voronoi_geopandas_with_plot.png b/tests/baseline/test_voronoi_geopandas_with_plot.png index a663278..5caa98a 100644 Binary files a/tests/baseline/test_voronoi_geopandas_with_plot.png and b/tests/baseline/test_voronoi_geopandas_with_plot.png differ diff --git a/tests/baseline/test_voronoi_italy_with_plot.png b/tests/baseline/test_voronoi_italy_with_plot.png deleted file mode 100644 index 285c40b..0000000 Binary files a/tests/baseline/test_voronoi_italy_with_plot.png and /dev/null differ diff --git a/tests/baseline/test_voronoi_italy_with_plot_10-False.png b/tests/baseline/test_voronoi_italy_with_plot_10-False.png new file mode 100644 index 0000000..44bb8a0 Binary files /dev/null and b/tests/baseline/test_voronoi_italy_with_plot_10-False.png differ diff --git a/tests/baseline/test_voronoi_italy_with_plot_10-True.png b/tests/baseline/test_voronoi_italy_with_plot_10-True.png new file mode 100644 index 0000000..89923cb Binary files /dev/null and b/tests/baseline/test_voronoi_italy_with_plot_10-True.png differ diff --git a/tests/baseline/test_voronoi_italy_with_plot_100-False.png b/tests/baseline/test_voronoi_italy_with_plot_100-False.png new file mode 100644 index 0000000..dba3c5f Binary files /dev/null and b/tests/baseline/test_voronoi_italy_with_plot_100-False.png differ diff --git a/tests/baseline/test_voronoi_italy_with_plot_100-True.png b/tests/baseline/test_voronoi_italy_with_plot_100-True.png new file mode 100644 index 0000000..c6ed126 Binary files /dev/null and b/tests/baseline/test_voronoi_italy_with_plot_100-True.png differ diff --git a/tests/baseline/test_voronoi_italy_with_plot_1000-False.png b/tests/baseline/test_voronoi_italy_with_plot_1000-False.png new file mode 100644 index 0000000..562abeb Binary files /dev/null and b/tests/baseline/test_voronoi_italy_with_plot_1000-False.png differ diff --git a/tests/baseline/test_voronoi_italy_with_plot_1000-True.png b/tests/baseline/test_voronoi_italy_with_plot_1000-True.png new file mode 100644 index 0000000..bbb9887 Binary files /dev/null and b/tests/baseline/test_voronoi_italy_with_plot_1000-True.png differ diff --git a/tests/baseline/test_voronoi_italy_with_plot_20-False.png b/tests/baseline/test_voronoi_italy_with_plot_20-False.png new file mode 100644 index 0000000..a436b05 Binary files /dev/null and b/tests/baseline/test_voronoi_italy_with_plot_20-False.png differ diff --git a/tests/baseline/test_voronoi_italy_with_plot_20-True.png b/tests/baseline/test_voronoi_italy_with_plot_20-True.png new file mode 100644 index 0000000..6b8579b Binary files /dev/null and b/tests/baseline/test_voronoi_italy_with_plot_20-True.png differ diff --git a/tests/baseline/test_voronoi_italy_with_plot_50-False.png b/tests/baseline/test_voronoi_italy_with_plot_50-False.png new file mode 100644 index 0000000..66104e1 Binary files /dev/null and b/tests/baseline/test_voronoi_italy_with_plot_50-False.png differ diff --git a/tests/baseline/test_voronoi_italy_with_plot_50-True.png b/tests/baseline/test_voronoi_italy_with_plot_50-True.png new file mode 100644 index 0000000..31f6fcf Binary files /dev/null and b/tests/baseline/test_voronoi_italy_with_plot_50-True.png differ diff --git a/tests/baseline/test_voronoi_italy_with_plot_500-False.png b/tests/baseline/test_voronoi_italy_with_plot_500-False.png new file mode 100644 index 0000000..39ec9cf Binary files /dev/null and b/tests/baseline/test_voronoi_italy_with_plot_500-False.png differ diff --git a/tests/baseline/test_voronoi_italy_with_plot_500-True.png b/tests/baseline/test_voronoi_italy_with_plot_500-True.png new file mode 100644 index 0000000..d4277b1 Binary files /dev/null and b/tests/baseline/test_voronoi_italy_with_plot_500-True.png differ diff --git a/tests/baseline/test_voronoi_spain_area_with_plot.png b/tests/baseline/test_voronoi_spain_area_with_plot.png index a5ce409..d344236 100644 Binary files a/tests/baseline/test_voronoi_spain_area_with_plot.png and b/tests/baseline/test_voronoi_spain_area_with_plot.png differ diff --git a/tests/baseline/test_voronoi_sweden_duplicate_points_with_plot.png b/tests/baseline/test_voronoi_sweden_duplicate_points_with_plot.png index ac372b3..84ad277 100644 Binary files a/tests/baseline/test_voronoi_sweden_duplicate_points_with_plot.png and b/tests/baseline/test_voronoi_sweden_duplicate_points_with_plot.png differ diff --git a/tests/test_geom.py b/tests/test_geom.py index 7ec2f4f..7f535fa 100644 --- a/tests/test_geom.py +++ b/tests/test_geom.py @@ -1,80 +1,72 @@ -from math import pi, isclose -from itertools import permutations +""" +Tests for the private _geom module. -from hypothesis import given -import numpy as np - -from ._testtools import real_coords_2d -from geovoronoi._geom import angle_between_pts, inner_angle_between_vecs, polygon_around_center, calculate_polygon_areas - - -@given(a=real_coords_2d(), b=real_coords_2d()) -def test_inner_angle_between_vecs(a, b): - a = np.array(a) - b = np.array(b) - origin = np.array((0, 0)) - - ang = inner_angle_between_vecs(a, b) +Author: Markus Konrad +""" - if np.allclose(a, origin, rtol=0) or np.allclose(b, origin, rtol=0): - assert np.isnan(ang) - else: - assert 0 <= ang <= pi - - -@given(a=real_coords_2d(), b=real_coords_2d()) -def test_angle_between_pts(a, b): - a = np.array(a) - b = np.array(b) - - ang = angle_between_pts(a, b) +import pytest +import numpy as np - if np.allclose(a, b, rtol=0): - assert np.isnan(ang) +from geovoronoi._geom import calculate_polygon_areas, line_segment_intersection + + +@pytest.mark.parametrize( + 'l_off, l_dir, segm_a, segm_b, expected', + [ + ([4, 7], [12, -4], [1, 1], [17, 5], [13, 4]), + # overlapping + ([1, 1], [2, 1], [2, 1.5], [4, 2.5], [2, 1.5]), + ([1, 1], [2, 1], [1, 1], [3, 2], [1, 1]), + ([2, 1.5], [2, 1], [2, 1.5], [4, 2.5], [2, 1.5]), + ([3, 2], [2, 1], [2, 1.5], [4, 2.5], [4, 2.5]), + ([3, 2], [2, 1], [4, 2.5], [2, 1.5], [4, 2.5]), + ([4, 2.5], [2, 1], [2, 1.5], [4, 2.5], [4, 2.5]), + ([5, 3], [2, 1], [2, 1.5], [4, 2.5], None), # "wrong" direction + ([5, 3], [-2, -1], [2, 1.5], [4, 2.5], [4, 2.5]), + ([-1, 0], [-2, -1], [2, 1.5], [4, 2.5], None), # "wrong" direction + # parallel + ([1, 3], [2, 1], [2, 1.5], [4, 2.5], None), + ([1, 3], [-2, -1], [2, 1.5], [4, 2.5], None), + # some edge-cases + ([0, 0], [0, 1], [0, 0], [0, 1], [0, 0]), + ([0, 0], [1, 0], [0, 0], [1, 0], [0, 0]), + ] +) +def test_line_segment_intersection(l_off, l_dir, segm_a, segm_b, expected): + res = line_segment_intersection(np.array(l_off), + np.array(l_dir), + np.array(segm_a), + np.array(segm_b)) + if expected is None: + assert res is None else: - assert 0 <= ang <= 2*pi - - -def test_polygon_around_center(): - points = np.array([[0, 0], [1, 0], [1, 1], [1, 0], [1, -1]]) - - for perm_ind in permutations(range(len(points))): - # many of these permutations do not form a valid polygon in that order - # `polygon_around_center` will make sure that the order of points is correct to create a polygon around the - # center of these points - poly = polygon_around_center(points[perm_ind,:]) - - assert poly.is_simple and poly.is_valid - - -def test_polygon_around_center_given_center_is_one_of_points(): - points = np.array([[0, 0], [1, 0], [1, 1], [1, 0], [0.5, 0.5]]) - - for perm_ind in permutations(range(len(points))): - # many of these permutations do not form a valid polygon in that order - # `polygon_around_center` will make sure that the order of points is correct to create a polygon around the - # center of these points - poly = polygon_around_center(points[perm_ind,:], center=(0.5, 0.5)) - - assert poly.is_simple and poly.is_valid + assert np.array_equal(res, np.array(expected)) def test_calculate_polygon_areas_empty(): - areas = calculate_polygon_areas([]) + areas = calculate_polygon_areas({}) assert len(areas) == 0 +def test_calculate_polygon_areas_nondict(): + with pytest.raises(ValueError): + calculate_polygon_areas([]) + + def test_calculate_polygon_areas_world(): import geopandas as gpd world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres')) world = world[world.continent != 'Antarctica'].to_crs(epsg=3395) # meters as unit! + geoms = {i: geom for i, geom in enumerate(world.geometry)} - areas = calculate_polygon_areas(world.geometry) + areas = calculate_polygon_areas(geoms) + assert isinstance(areas, dict) assert len(areas) == len(world) - assert all(0 <= a < 9e13 for a in areas) + assert all(0 < a < 9e13 for a in areas.values()) - areas_km2 = calculate_polygon_areas(world.geometry, m2_to_km2=True) + areas_km2 = calculate_polygon_areas(geoms, m2_to_km2=True) + assert isinstance(areas_km2, dict) assert len(areas_km2) == len(world) - assert all(isclose(a_m, a_km * 1e6) for a_m, a_km in zip(areas, areas_km2)) + assert all(np.isclose(a_m, a_km * 1e6) for a_m, a_km in zip(areas.values(), areas_km2.values())) diff --git a/tests/test_main.py b/tests/test_main.py index 18d8c75..a50c2d7 100644 --- a/tests/test_main.py +++ b/tests/test_main.py @@ -1,130 +1,224 @@ +""" +Tests for the main module. + +Author: Markus Konrad +""" + import numpy as np import geopandas as gpd -from shapely.geometry import Polygon, MultiPolygon +from shapely.geometry import Polygon, MultiPolygon, asPoint from shapely.ops import cascaded_union import pytest -from hypothesis import given +from hypothesis import given, settings import hypothesis.strategies as st 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 -np.random.seed(123) - #%% tests for individual functions @given(coords=coords_2d_array()) def test_coords_to_points_and_points_to_coords(coords): + # test for bijectivity of points_to_coords and coords_to_points assert np.array_equal(points_to_coords(coords_to_points(coords)), coords) +@pytest.mark.parametrize( + 'poly_to_pts, expected', + [ + ({5: [1], 2: [3]}, {1: 5, 3: 2}), + ({5: [], 2: [3]}, {3: 2}), + ({5: [], 2: []}, {}), + ({1: [1]}, {1: 1}), + ({1: [4, 5]}, {4: 1, 5: 1}), + ({5: [1], 2: [1]}, None) + ] +) +def test_get_points_to_poly_assignments(poly_to_pts, expected): + if expected is None: + with pytest.raises(ValueError): + points_to_region(poly_to_pts) + else: + assert points_to_region(poly_to_pts) == expected + @given(available_points=st.permutations(list(range(10))), n_poly=st.integers(0, 10)) -def test_get_points_to_poly_assignments(available_points, n_poly): +def test_get_points_to_poly_assignments_hypothesis(available_points, n_poly): + # generate poly to point assignments n_pts = len(available_points) if n_poly == 0: - poly_to_pts = [] - elif n_poly == 10: - poly_to_pts = [[x] for x in available_points] - else: + poly_to_pts = {} + elif n_poly == 10: # one to one assignment + poly_to_pts = dict(zip(range(n_poly), [[x] for x in available_points])) + else: # one to N assignment (we have duplicate points) pts_per_poly = n_pts // n_poly - poly_to_pts = [] + poly_to_pts = {} n_assigned = 0 + # try to evenly distribute point IDs to polys for p in range(0, n_poly): - poly_to_pts.append([available_points[i] for i in range(p * pts_per_poly, (p+1) * pts_per_poly)]) + poly_to_pts[p] = [available_points[i] for i in range(p * pts_per_poly, (p+1) * pts_per_poly)] n_assigned += pts_per_poly + # fill up if n_assigned < n_pts: - poly_to_pts[-1].extend([available_points[i] for i in range(n_assigned, n_pts)]) + poly_to_pts[n_poly-1].extend([available_points[i] for i in range(n_assigned, n_pts)]) if n_poly > 0: - assert set(sum(poly_to_pts, [])) == set(available_points) + 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, list) + assert isinstance(pts_to_poly, dict) if n_poly == 0: assert len(pts_to_poly) == 0 else: assert len(pts_to_poly) == n_pts - assert all([0 <= i_poly < n_poly for i_poly in pts_to_poly]) - - -#%% tests against fixed issues - -def test_issue_7a(): - centroids = np.array([[537300, 213400], [538700, 213700], [536100, 213400]]) - polygon = Polygon([[540000, 214100], [535500, 213700], [535500, 213000], [539000, 213200]]) - poly_shapes, pts, poly_to_pt_assignments = voronoi_regions_from_coords(centroids, polygon) - - assert isinstance(poly_shapes, list) - assert 0 < len(poly_shapes) <= len(centroids) - assert all([isinstance(p, (Polygon, MultiPolygon)) for p in poly_shapes]) - - assert np.array_equal(points_to_coords(pts), centroids) - - assert isinstance(poly_to_pt_assignments, list) - assert len(poly_to_pt_assignments) == len(poly_shapes) - assert all([isinstance(assign, list) for assign in poly_to_pt_assignments]) - assert all([len(assign) == 1 for assign in poly_to_pt_assignments]) # in this case there is a 1:1 correspondance - - -@pytest.mark.mpl_image_compare -def test_issue_7b(): - centroids = np.array([[496712, 232672], [497987, 235942], [496425, 230252], [497482, 234933], - [499331, 238351], [496081, 231033], [497090, 233846], [496755, 231645], - [498604, 237018]]) - polygon = Polygon([[495555, 230875], [496938, 235438], [499405, 239403], [499676, 239474], - [499733, 237877], [498863, 237792], [499120, 237335], [498321, 235010], - [497295, 233185], [497237, 231359], [496696, 229620], [495982, 230047], - [496154, 230347], [496154, 230347], [495555, 230875]]) - - poly_shapes, pts, poly_to_pt_assignments = voronoi_regions_from_coords(centroids, polygon) - - assert isinstance(poly_shapes, list) - assert 0 < len(poly_shapes) <= len(centroids) - assert all([isinstance(p, (Polygon, MultiPolygon)) for p in poly_shapes]) - - assert np.array_equal(points_to_coords(pts), centroids) - - assert isinstance(poly_to_pt_assignments, list) - assert len(poly_to_pt_assignments) == len(poly_shapes) - assert all([isinstance(assign, list) for assign in poly_to_pt_assignments]) - assert all([len(assign) == 1 for assign in poly_to_pt_assignments]) # in this case there is a 1:1 correspondance - - fig, ax = subplot_for_map() - plot_voronoi_polys_with_points_in_area(ax, polygon, poly_shapes, centroids, poly_to_pt_assignments) - - return fig - - -#%% realistic full tests with plotting + assert set(list(pts_to_poly.keys())) == set(available_points) + assert set(list(pts_to_poly.values())) == set(list(range(n_poly))) +@settings(deadline=10000) +@given(n_pts=st.integers(0, 200), + per_geom=st.booleans(), + return_unassigned_pts=st.booleans(), + results_per_geom=st.booleans()) +def test_voronoi_regions_from_coords_italy(n_pts, per_geom, return_unassigned_pts, results_per_geom): + area_shape = _get_country_shape('Italy') + n_geoms = len(area_shape.geoms) # number of geometries (is 3 -- main land Italy plus Sardinia and Sicilia) + # put random coordinates inside shape + coords = _rand_coords_in_shape(area_shape, n_pts) + n_pts = len(coords) # number of random points inside shape + + if n_pts < 2: # check ValueError when less than 2 points are submitted + with pytest.raises(ValueError): + voronoi_regions_from_coords(coords, area_shape, + per_geom=per_geom, + return_unassigned_points=return_unassigned_pts, + results_per_geom=results_per_geom) + + return + + # generate Voronoi region polygons + res = voronoi_regions_from_coords(coords, area_shape, + per_geom=per_geom, + return_unassigned_points=return_unassigned_pts, + results_per_geom=results_per_geom) + + # in any case, this must return a tuple of results + assert isinstance(res, tuple) + + if return_unassigned_pts: # additionally expect set of unassigned points + assert len(res) == 3 + region_polys, region_pts, unassigned_pts = res + assert isinstance(unassigned_pts, set) + assert all([pt in range(n_pts) for pt in unassigned_pts]) + else: + assert len(res) == 2 + region_polys, region_pts = res + unassigned_pts = None + + # check general result structure + assert isinstance(region_polys, dict) + assert isinstance(region_pts, dict) + assert list(region_polys.keys()) == list(region_pts.keys()) + + if results_per_geom: # expect a dict that maps geom ID to results + if not per_geom: # if geoms are not treated separately, there's only one geom ID + assert list(region_polys.keys()) == list(region_pts.keys()) == [0] + + # iterate through geoms + for i_geom in region_polys.keys(): + # get Voronoi polygons + assert 0 <= i_geom < n_geoms + region_polys_in_geom = region_polys[i_geom] + assert isinstance(region_polys_in_geom, dict) + + # get Voronoi region -> points assignments + region_pts_in_geom = region_pts[i_geom] + assert isinstance(region_pts_in_geom, dict) + assert list(region_polys_in_geom.keys()) == list(region_pts_in_geom.keys()) + + # check region polygons + if region_pts_in_geom: + if per_geom: + geom_area = area_shape.geoms[i_geom].area + else: + geom_area = sum(g.area for g in area_shape.geoms) + _check_region_polys(region_polys_in_geom.values(), region_pts_in_geom.values(), coords, + expected_sum_area=geom_area) + else: + # no polys generated -> must be insufficient number of points in geom + pass + else: # not results_per_geom + # results are *not* given per geom ID + assert len(region_polys) <= n_pts + assert len(region_pts) == len(region_polys) + + # points to region assignments + 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 + + # check result structure + assert isinstance(region_polys, dict) + assert isinstance(region_pts, dict) + assert list(region_polys.keys()) == list(region_pts.keys()) + + # check region polygons + if region_polys: + if per_geom: + geom_area = None # can't determine this here + else: + geom_area = sum(g.area for g in area_shape.geoms) + _check_region_polys(region_polys.values(), region_pts.values(), coords, expected_sum_area=geom_area) + else: + # no polys generated -> must be insufficient number of points + assert n_pts < 4 + + # fig, ax = subplot_for_map(show_x_axis=True, show_y_axis=True) + # plot_voronoi_polys_with_points_in_area(ax, area_shape, region_polys, coords, region_pts, + # point_labels=list(map(str, range(len(coords)))), + # voronoi_labels=list(map(str, region_polys.keys()))) + # fig.show() + + +# #%% realistic full tests with plotting + + +@pytest.mark.parametrize( + 'n_pts,per_geom', [ + (10, True), (10, False), + (20, True), (20, False), + (50, True), (50, False), + (100, True), (100, False), + (500, True), (500, False), + (1000, True), (1000, False), + ] +) @pytest.mark.mpl_image_compare -def test_voronoi_italy_with_plot(): +def test_voronoi_italy_with_plot(n_pts, per_geom): area_shape = _get_country_shape('Italy') - coords = _rand_coords_in_shape(area_shape, 100) - poly_shapes, pts, poly_to_pt_assignments = voronoi_regions_from_coords(coords, area_shape) + coords = _rand_coords_in_shape(area_shape, n_pts) - assert isinstance(poly_shapes, list) - assert 0 < len(poly_shapes) <= 100 - assert all([isinstance(p, (Polygon, MultiPolygon)) for p in poly_shapes]) + # generate Voronoi regions + region_polys, region_pts = voronoi_regions_from_coords(coords, area_shape, per_geom=per_geom) - assert np.array_equal(points_to_coords(pts), coords) + # full checks for voronoi_regions_from_coords() are done in test_voronoi_regions_from_coords_italy() - assert isinstance(poly_to_pt_assignments, list) - assert len(poly_to_pt_assignments) == len(poly_shapes) - assert all([isinstance(assign, list) for assign in poly_to_pt_assignments]) - assert all([len(assign) == 1 for assign in poly_to_pt_assignments]) # in this case there is a 1:1 correspondance + assert isinstance(region_polys, dict) + assert isinstance(region_pts, dict) + assert len(region_polys) == len(region_pts) + assert 0 < len(region_polys) <= n_pts + # generate plot fig, ax = subplot_for_map() - plot_voronoi_polys_with_points_in_area(ax, area_shape, poly_shapes, coords, poly_to_pt_assignments) + plot_voronoi_polys_with_points_in_area(ax, area_shape, region_polys, coords, region_pts, + point_labels=list(map(str, range(len(coords))))) return fig @@ -133,29 +227,26 @@ def test_voronoi_italy_with_plot(): def test_voronoi_spain_area_with_plot(): area_shape = _get_country_shape('Spain') coords = _rand_coords_in_shape(area_shape, 20) - poly_shapes, pts, poly_to_pt_assignments = voronoi_regions_from_coords(coords, area_shape) - assert isinstance(poly_shapes, list) - assert 0 < len(poly_shapes) <= 20 - assert all([isinstance(p, (Polygon, MultiPolygon)) for p in poly_shapes]) + # generate Voronoi regions + region_polys, region_pts = voronoi_regions_from_coords(coords, area_shape) - assert np.array_equal(points_to_coords(pts), coords) + # full checks for voronoi_regions_from_coords() are done in test_voronoi_regions_from_coords_italy() - assert isinstance(poly_to_pt_assignments, list) - assert len(poly_to_pt_assignments) == len(poly_shapes) - assert all([isinstance(assign, list) for assign in poly_to_pt_assignments]) - assert all([len(assign) == 1 for assign in poly_to_pt_assignments]) # in this case there is a 1:1 correspondance + assert isinstance(region_polys, dict) + assert isinstance(region_pts, dict) + assert len(region_polys) == len(region_pts) + assert 0 < len(region_polys) <= 20 - poly_areas = calculate_polygon_areas(poly_shapes, m2_to_km2=True) # converts m² to km² - assert isinstance(poly_areas, np.ndarray) - assert np.issubdtype(poly_areas.dtype, np.float_) - assert len(poly_areas) == len(poly_shapes) - assert np.all(poly_areas > 0) + # generate covered area + region_areas = calculate_polygon_areas(region_polys, m2_to_km2=True) # converts m² to km² + assert isinstance(region_areas, dict) + assert set(region_areas.keys()) == set(region_polys.keys()) + # generate plot fig, ax = subplot_for_map(show_x_axis=True, show_y_axis=True) - - voronoi_labels = ['%d km²' % round(a) for a in poly_areas] - plot_voronoi_polys_with_points_in_area(ax, area_shape, poly_shapes, coords, poly_to_pt_assignments, + voronoi_labels = {k: '%d km²' % round(a) for k, a in region_areas.items()} + plot_voronoi_polys_with_points_in_area(ax, area_shape, region_polys, coords, region_pts, voronoi_labels=voronoi_labels, voronoi_label_fontsize=7, voronoi_label_color='gray') @@ -179,22 +270,17 @@ def test_voronoi_geopandas_with_plot(): 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) - - assert isinstance(poly_shapes, list) - assert 0 < len(poly_shapes) <= len(coords) - assert all([isinstance(p, (Polygon, MultiPolygon)) for p in poly_shapes]) + region_polys, region_pts = voronoi_regions_from_coords(coords, south_am_shape, per_geom=False) - assert np.array_equal(points_to_coords(pts), coords) + # full checks for voronoi_regions_from_coords() are done in test_voronoi_regions_from_coords_italy() - assert isinstance(poly_to_pt_assignments, list) - assert len(poly_to_pt_assignments) == len(poly_shapes) - assert all([isinstance(assign, list) for assign in poly_to_pt_assignments]) - assert all([len(assign) == 1 for assign in poly_to_pt_assignments]) # in this case there is a 1:1 correspondance + assert isinstance(region_polys, dict) + assert isinstance(region_pts, dict) + assert len(region_polys) == len(region_pts) == len(coords) + # generate plot 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, region_polys, coords, region_pts) return fig @@ -207,41 +293,78 @@ def test_voronoi_sweden_duplicate_points_with_plot(): # duplicate a few points rand_dupl_ind = np.random.randint(len(coords), size=10) coords = np.concatenate((coords, coords[rand_dupl_ind])) + n_pts = len(coords) - poly_shapes, pts, poly_to_pt_assignments = voronoi_regions_from_coords(coords, area_shape, - accept_n_coord_duplicates=10) - - assert isinstance(poly_shapes, list) - assert 0 < len(poly_shapes) <= 20 - assert all([isinstance(p, (Polygon, MultiPolygon)) for p in poly_shapes]) + # generate Voronoi regions + region_polys, region_pts = voronoi_regions_from_coords(coords, area_shape) - assert np.array_equal(points_to_coords(pts), coords) + # full checks for voronoi_regions_from_coords() are done in test_voronoi_regions_from_coords_italy() - assert isinstance(poly_to_pt_assignments, list) - assert len(poly_to_pt_assignments) == len(poly_shapes) - assert all([isinstance(assign, list) for assign in poly_to_pt_assignments]) - assert all([0 < len(assign) <= 10 for assign in poly_to_pt_assignments]) # in this case there is not - # everywhere a 1:1 correspondance + assert isinstance(region_polys, dict) + assert isinstance(region_pts, dict) + assert 0 < len(region_polys) <= n_pts + assert 0 < len(region_pts) <= n_pts - pts_to_poly_assignments = np.array(get_points_to_poly_assignments(poly_to_pt_assignments)) + assert all([0 < len(pts_in_region) <= 10 for pts_in_region in region_pts.values()]) - # 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 region_pts.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: (1, 0, 0) if len(pt_indices) > 1 else (0, 0, 1) + for i_poly, pt_indices in region_pts.items()} + # generate plot 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, region_polys, 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=(0, 0, 0, 1), point_labels=pt_labels, - points_markersize=np.array(count_per_pt)*10) + points_markersize=np.square(np.array(list(count_per_pt.values()))) * 10) + + return fig + + +#%% tests against fixed issues + +def test_issue_7a(): + centroids = np.array([[537300, 213400], [538700, 213700], [536100, 213400]]) + n_pts = len(centroids) + polygon = Polygon([[540000, 214100], [535500, 213700], [535500, 213000], [539000, 213200]]) + region_polys, region_pts = voronoi_regions_from_coords(centroids, polygon) + + assert isinstance(region_polys, dict) + assert isinstance(region_pts, dict) + assert len(region_polys) == len(region_pts) == n_pts + + assert all([len(pts_in_region) == 1 for pts_in_region in region_pts.values()]) # no duplicates + + +@pytest.mark.mpl_image_compare +def test_issue_7b(): + centroids = np.array([[496712, 232672], [497987, 235942], [496425, 230252], [497482, 234933], + [499331, 238351], [496081, 231033], [497090, 233846], [496755, 231645], + [498604, 237018]]) + n_pts = len(centroids) + polygon = Polygon([[495555, 230875], [496938, 235438], [499405, 239403], [499676, 239474], + [499733, 237877], [498863, 237792], [499120, 237335], [498321, 235010], + [497295, 233185], [497237, 231359], [496696, 229620], [495982, 230047], + [496154, 230347], [496154, 230347], [495555, 230875]]) + + region_polys, region_pts = voronoi_regions_from_coords(centroids, polygon) + + assert isinstance(region_polys, dict) + assert isinstance(region_pts, dict) + assert len(region_polys) == len(region_pts) == n_pts + + assert all([len(pts_in_region) == 1 for pts_in_region in region_pts.values()]) # no duplicates + + fig, ax = subplot_for_map() + plot_voronoi_polys_with_points_in_area(ax, polygon, region_polys, centroids, region_pts) return fig @@ -257,6 +380,8 @@ def _get_country_shape(country): def _rand_coords_in_shape(area_shape, n_points): + np.random.seed(123) + # generate some random points within the bounds minx, miny, maxx, maxy = area_shape.bounds @@ -267,3 +392,24 @@ def _rand_coords_in_shape(area_shape, n_points): # 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 return points_to_coords(pts) + + +def _check_region_polys(region_polys, region_pts, coords, expected_sum_area, + contains_check_tol=1, area_check_tol=0.1): + # check validity of each region's polygon, check that all assigned points are inside this polygon and + # check that sum of polygons' area matches `expected_sum_area` + sum_area = 0 + for poly, pt_indices in zip(region_polys, region_pts): + assert isinstance(poly, (Polygon, MultiPolygon)) and poly.is_valid and not poly.is_empty + if contains_check_tol != 0: + polybuf = poly.buffer(contains_check_tol) + if polybuf.is_empty or not polybuf.is_valid: # this may happen due to buffering + polybuf = poly + else: + polybuf = poly + assert all([polybuf.contains(asPoint(coords[i_pt])) for i_pt in pt_indices]) + + sum_area += poly.area + + if expected_sum_area is not None: + assert abs(1 - sum_area / expected_sum_area) <= area_check_tol diff --git a/tox.ini b/tox.ini index 989c1dc..ef0ef24 100644 --- a/tox.ini +++ b/tox.ini @@ -4,12 +4,11 @@ # and then run "tox" from this directory. [tox] -envlist = py36, py37, py38 +envlist = py36, py37, py38, py39 [testenv] deps = .[test] extras = plotting commands = pytest - pytest --mpl-generate-path=.tox/mpl-baseline-images - pytest --mpl + pytest --mpl --mpl-baseline-path=tests/baseline