Dividing North America into pieces

Question

I'm trying to create the voronoi diagram of a map of North America, which means effectively slicing the country into pieces based on the location of its capitals. To do so, I'm using Geopandas to get the geographic data on North America, and then using the GeoVoronoi library to create a Voronoi diagram out of that:


import matplotlib.pyplot as plt
import geopandas as gpd
from shapely.ops import cascaded_union

from geovoronoi.plotting import subplot_for_map, plot_voronoi_polys_with_points_in_area
from geovoronoi import voronoi_regions_from_coords, points_to_coords


logging.basicConfig(level=logging.INFO)
geovoronoi_log = logging.getLogger('geovoronoi')
geovoronoi_log.setLevel(logging.INFO)
geovoronoi_log.propagate = True

#
# load geo data
#

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)
north_am = world[world.continent == 'North America'].to_crs(epsg=3395)
cities = cities.to_crs(north_am.crs)   # convert city coordinates to same CRS!

# create the bounding shape as union of all South American countries' shapes
north_am_shape = cascaded_union(north_am.geometry)
north_am_cities = cities[cities.geometry.within(north_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
#

# convert the pandas Series of Point objects to NumPy array of coordinates
coords = points_to_coords(north_am_cities.geometry)

# calculate the regions
poly_shapes, pts, poly_to_pt_assignments = voronoi_regions_from_coords(coords, north_am_shape)


#
# Plotting
#

fig, ax = subplot_for_map()

plot_voronoi_polys_with_points_in_area(ax, north_am_shape, poly_shapes, pts)

ax.set_title('Cities data for South America from GeoPandas
and Voronoi regions around them')

plt.tight_layout()
plt.savefig('using_geopandas.png')
plt.show()

Most of this code is taken straight out of Geovoronoi documentation. And yet, when I run it, I receive the following error:

user_7384740 · Answer

The error you are receiving arises from the fact that the city info you are pulling contains few cities North America or that they are not properly being recognized as within the boundaries of North America. Your question is about creating a Voronoi diagram based on capitals so I have included a link to a data set for US capitals so that you can test an example with a reliable number of cities:

import matplotlib.pyplot as plt
import numpy as np
import geopandas as gpd
from geovoronoi.plotting import subplot_for_map, plot_voronoi_polys_with_points_in_area
from geovoronoi import voronoi_regions_from_coords

cities = gpd.read_file('us-state-capitals.csv')

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
usa = world[world.name == 'United States of America']
usa = usa.to_crs(epsg=3857)
usa_shape =  usa.iloc[0].geometry

coords = np.array(list(zip(cities.Shape_X,cities.Shape_Y)), dtype='float')

poly_shapes, pts, poly_to_pt_assignments = voronoi_regions_from_coords(coords, usa_shape)

fig, ax = subplot_for_map()
plot_voronoi_polys_with_points_in_area(ax, usa_shape, poly_shapes, coords)
ax.set_title('Cities data for South America from GeoPandas
and Voronoi regions around them')
plt.tight_layout()
plt.savefig('using_geopandas.png')
plt.show()

Producing:

For North America you can download a cities CSV and use the following code:

import matplotlib.pyplot as plt
import geopandas as gpd
from shapely.ops import cascaded_union
from geovoronoi.plotting import subplot_for_map, plot_voronoi_polys_with_points_in_area
from geovoronoi import voronoi_regions_from_coords, points_to_coords

cities = gpd.read_file('world_populated_cities.csv')
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
na = world[world.continent == 'North America']
na = na.to_crs(epsg=3857)
cities.geometry.to_crs(epsg=3857)

na_shape = cascaded_union(na.geometry)
cities = cities.to_crs(na.crs)   # convert city coordinates to same CRS!
cities = cities[cities.geometry.within(na_shape)]

coords = points_to_coords(cities.geometry)
poly_shapes, pts, poly_to_pt_assignments = voronoi_regions_from_coords(coords, na_shape)

fig, ax = subplot_for_map()
plot_voronoi_polys_with_points_in_area(ax, na_shape, poly_shapes, coords)
ax.set_title('Cities data for South America from GeoPandas
and Voronoi regions around them')
plt.tight_layout()
plt.savefig('using_geopandas.png')
plt.show()

Producing:

Dividing North America into pieces

1 Answers