diff --git a/geovoronoi/_voronoi.py b/geovoronoi/_voronoi.py index 36f101f..eee9a2b 100644 --- a/geovoronoi/_voronoi.py +++ b/geovoronoi/_voronoi.py @@ -32,7 +32,8 @@ def points_to_coords(pts): return np.array([p.coords[0] for p in pts]) -def voronoi_regions_from_coords(coords, geo_shape, per_geom=True, return_unassigned_points=False): +def voronoi_regions_from_coords(coords, geo_shape, per_geom=True, return_unassigned_points=False, + results_per_geom=False): """ 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 @@ -59,6 +60,9 @@ def voronoi_regions_from_coords(coords, geo_shape, per_geom=True, return_unassig logger.info('running Voronoi tesselation for %d points / treating geoms separately: %s' % (len(coords), per_geom)) + if len(coords) < 2: + raise ValueError('insufficient number of points provided in `coords`') + if isinstance(coords, np.ndarray): pts = coords_to_points(coords) else: @@ -78,12 +82,18 @@ def voronoi_regions_from_coords(coords, geo_shape, per_geom=True, return_unassig geoms = geo_shape.geoms pts_indices = set(range(len(pts))) - region_polys = {} - region_pts = {} + geom_region_polys = {} + geom_region_pts = {} unassigned_pts = set() + for i_geom, geom in enumerate(geoms): pts_in_geom = [i for i in pts_indices if geom.contains(pts[i])] + + geom_region_polys[i_geom] = {} + geom_region_pts[i_geom] = {} + logger.info('%d of %d points in geometry #%d of %d' % (len(pts_in_geom), len(pts), i_geom + 1, len(geoms))) + if not pts_in_geom: continue @@ -103,15 +113,29 @@ def voronoi_regions_from_coords(coords, geo_shape, per_geom=True, return_unassig logger.info('generating Voronoi region polygons') geom_polys, geom_pts = region_polygons_from_voronoi(vor, geom, return_point_assignments=True) - 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()))) + # map 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() + + geom_region_polys[i_geom] = geom_polys + geom_region_pts[i_geom] = geom_pts + + if results_per_geom: + region_polys = geom_region_polys + region_pts = geom_region_pts + else: + region_polys = {} + region_pts = {} + + for i_geom, geom_polys in geom_region_polys.items(): + geom_pts = geom_region_pts[i_geom] + + 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()))) - pt_ids_mapping = dict(zip(range(len(pts_in_geom)), pts_in_geom)) - region_pts.update( - {reg_id: [pt_ids_mapping[pt_id] for pt_id in geom_pts[old_reg_id]] - for reg_id, old_reg_id in region_ids_mapping.items()} - ) + region_pts.update({reg_id: geom_pts[old_reg_id] for reg_id, old_reg_id in region_ids_mapping.items()}) if return_unassigned_points: return region_polys, region_pts, unassigned_pts diff --git a/tests/test_main.py b/tests/test_main.py index 4fe6c25..79992e7 100644 --- a/tests/test_main.py +++ b/tests/test_main.py @@ -1,9 +1,9 @@ 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 @@ -78,232 +78,271 @@ def test_get_points_to_poly_assignments_hypothesis(available_points, n_poly): assert set(list(pts_to_poly.values())) == set(list(range(n_poly))) -@pytest.mark.parametrize( - 'n_pts, per_geom, return_unassigned_pts', - [ - (100, True, False), - ] -) -@pytest.mark.mpl_image_compare -def test_voronoi_regions_from_coords_italy(n_pts, per_geom, return_unassigned_pts): +@settings(deadline=10000) +@given(n_pts=st.integers(0, 20), + 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) coords = _rand_coords_in_shape(area_shape, n_pts) n_pts = len(coords) # number of random points inside shape - poly_shapes, poly_to_pt_assignments = voronoi_regions_from_coords(coords, area_shape, - per_geom=per_geom, - return_unassigned_points=return_unassigned_pts) - - assert isinstance(poly_shapes, dict) - assert len(poly_shapes) == n_pts - assert all([isinstance(p, (Polygon, MultiPolygon)) for p in poly_shapes]) - - assert np.array_equal(points_to_coords(pts), coords) - - assert isinstance(poly_to_pt_assignments, dict) - 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, area_shape, poly_shapes, coords, poly_to_pt_assignments) - - return fig - - - -#%% realistic full tests with plotting - - -@pytest.mark.parametrize( - 'n_pts', [5, 10, 20, 50, 100, 105, 1000] -) -@pytest.mark.mpl_image_compare -def test_voronoi_italy_with_plot(n_pts): - area_shape = _get_country_shape('Italy') - coords = _rand_coords_in_shape(area_shape, n_pts) - poly_shapes, poly_to_pt_assignments = voronoi_regions_from_coords(coords, area_shape) - - assert isinstance(poly_shapes, list) - assert 0 < len(poly_shapes) <= 100 - assert all([isinstance(p, (Polygon, MultiPolygon)) for p in poly_shapes]) - - assert np.array_equal(points_to_coords(pts), coords) - - 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, area_shape, poly_shapes, coords, poly_to_pt_assignments) - - return fig - - -@pytest.mark.mpl_image_compare -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]) - - assert np.array_equal(points_to_coords(pts), coords) - - 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 - - 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) - - 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=voronoi_labels, voronoi_label_fontsize=7, - voronoi_label_color='gray') - - return fig - - -@pytest.mark.mpl_image_compare -def test_voronoi_geopandas_with_plot(): - world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres')) - cities = gpd.read_file(gpd.datasets.get_path('naturalearth_cities')) - - # focus on South America, convert to World Mercator (unit: meters) - south_am = world[world.continent == 'South America'].to_crs(epsg=3395) - cities = cities.to_crs(south_am.crs) # convert city coordinates to same CRS! - - # create the bounding shape as union of all South American countries' shapes - south_am_shape = cascaded_union(south_am.geometry) - south_am_cities = cities[cities.geometry.within(south_am_shape)] # reduce to cities in South America - # convert the pandas Series of Point objects to NumPy array of coordinates - 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]) - - assert np.array_equal(points_to_coords(pts), coords) - - 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, south_am_shape, poly_shapes, pts, poly_to_pt_assignments) - - return fig - - -@pytest.mark.mpl_image_compare -def test_voronoi_sweden_duplicate_points_with_plot(): - area_shape = _get_country_shape('Sweden') - coords = _rand_coords_in_shape(area_shape, 20) - - # duplicate a few points - rand_dupl_ind = np.random.randint(len(coords), size=10) - coords = np.concatenate((coords, coords[rand_dupl_ind])) - - 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]) - - assert np.array_equal(points_to_coords(pts), coords) - - 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 - - pts_to_poly_assignments = np.array(get_points_to_poly_assignments(poly_to_pt_assignments)) - - # make point labels: counts of duplicates per points - count_per_pt = [sum(pts_to_poly_assignments == i_poly) for i_poly in pts_to_poly_assignments] - pt_labels = list(map(str, count_per_pt)) - - # highlight voronoi regions with point duplicates - count_per_poly = np.array(list(map(len, poly_to_pt_assignments))) - vor_colors = np.repeat('blue', len(poly_shapes)) # default color - vor_colors[count_per_poly > 1] = 'red' # hightlight color - - fig, ax = subplot_for_map() - - plot_voronoi_polys_with_points_in_area(ax, area_shape, poly_shapes, coords, - plot_voronoi_opts={'alpha': 0.2}, - plot_points_opts={'alpha': 0.4}, - voronoi_color=list(vor_colors), - point_labels=pt_labels, - points_markersize=np.array(count_per_pt)*10) - - return fig - - -#%% 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) + if n_pts < 2: + 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) - 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 + return - fig, ax = subplot_for_map() - plot_voronoi_polys_with_points_in_area(ax, polygon, poly_shapes, centroids, poly_to_pt_assignments) + res = voronoi_regions_from_coords(coords, area_shape, + per_geom=per_geom, + return_unassigned_points=return_unassigned_pts, + results_per_geom=results_per_geom) - return fig + assert isinstance(res, tuple) + if return_unassigned_pts: + 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 + + assert isinstance(region_polys, dict) + assert isinstance(region_pts, dict) + assert list(region_polys.keys()) == list(region_pts.keys()) + + if results_per_geom: + if not per_geom: + assert list(region_polys.keys()) == list(region_pts.keys()) == [0] + + for i_geom in region_polys.keys(): + assert 0 <= i_geom < n_geoms + region_polys_in_geom = region_polys[i_geom] + assert isinstance(region_polys_in_geom, dict) + 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()) + + 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 + pass + else: # not results_per_geom + assert len(region_polys) <= n_pts + assert len(region_pts) == len(region_polys) + pts_region = get_points_to_poly_assignments(region_pts) + if unassigned_pts is not None: + assert set(range(n_pts)) - set(pts_region.keys()) == unassigned_pts + + + +# #%% realistic full tests with plotting +# +# +# @pytest.mark.parametrize( +# 'n_pts', [5, 10, 20, 50, 100, 105, 1000] +# ) +# @pytest.mark.mpl_image_compare +# def test_voronoi_italy_with_plot(n_pts): +# area_shape = _get_country_shape('Italy') +# coords = _rand_coords_in_shape(area_shape, n_pts) +# poly_shapes, poly_to_pt_assignments = voronoi_regions_from_coords(coords, area_shape) +# +# assert isinstance(poly_shapes, list) +# assert 0 < len(poly_shapes) <= 100 +# assert all([isinstance(p, (Polygon, MultiPolygon)) for p in poly_shapes]) +# +# assert np.array_equal(points_to_coords(pts), coords) +# +# 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, area_shape, poly_shapes, coords, poly_to_pt_assignments) +# +# return fig +# +# +# @pytest.mark.mpl_image_compare +# 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]) +# +# assert np.array_equal(points_to_coords(pts), coords) +# +# 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 +# +# 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) +# +# 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=voronoi_labels, voronoi_label_fontsize=7, +# voronoi_label_color='gray') +# +# return fig +# +# +# @pytest.mark.mpl_image_compare +# def test_voronoi_geopandas_with_plot(): +# world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres')) +# cities = gpd.read_file(gpd.datasets.get_path('naturalearth_cities')) +# +# # focus on South America, convert to World Mercator (unit: meters) +# south_am = world[world.continent == 'South America'].to_crs(epsg=3395) +# cities = cities.to_crs(south_am.crs) # convert city coordinates to same CRS! +# +# # create the bounding shape as union of all South American countries' shapes +# south_am_shape = cascaded_union(south_am.geometry) +# south_am_cities = cities[cities.geometry.within(south_am_shape)] # reduce to cities in South America +# +# # convert the pandas Series of Point objects to NumPy array of coordinates +# 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]) +# +# assert np.array_equal(points_to_coords(pts), coords) +# +# 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, south_am_shape, poly_shapes, pts, poly_to_pt_assignments) +# +# return fig +# +# +# @pytest.mark.mpl_image_compare +# def test_voronoi_sweden_duplicate_points_with_plot(): +# area_shape = _get_country_shape('Sweden') +# coords = _rand_coords_in_shape(area_shape, 20) +# +# # duplicate a few points +# rand_dupl_ind = np.random.randint(len(coords), size=10) +# coords = np.concatenate((coords, coords[rand_dupl_ind])) +# +# 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]) +# +# assert np.array_equal(points_to_coords(pts), coords) +# +# 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 +# +# pts_to_poly_assignments = np.array(get_points_to_poly_assignments(poly_to_pt_assignments)) +# +# # make point labels: counts of duplicates per points +# count_per_pt = [sum(pts_to_poly_assignments == i_poly) for i_poly in pts_to_poly_assignments] +# pt_labels = list(map(str, count_per_pt)) +# +# # highlight voronoi regions with point duplicates +# count_per_poly = np.array(list(map(len, poly_to_pt_assignments))) +# vor_colors = np.repeat('blue', len(poly_shapes)) # default color +# vor_colors[count_per_poly > 1] = 'red' # hightlight color +# +# fig, ax = subplot_for_map() +# +# plot_voronoi_polys_with_points_in_area(ax, area_shape, poly_shapes, coords, +# plot_voronoi_opts={'alpha': 0.2}, +# plot_points_opts={'alpha': 0.4}, +# voronoi_color=list(vor_colors), +# point_labels=pt_labels, +# points_markersize=np.array(count_per_pt)*10) +# +# return fig +# +# +# #%% 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 +# #%% a few helper functions @@ -326,3 +365,18 @@ 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): + 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) + else: + polybuf = poly + assert all([polybuf.contains(asPoint(coords[i_pt])) for i_pt in pt_indices]) + + sum_area += poly.area + + #assert np.isclose(sum_area, expected_sum_area)