Geographic Data Science

How to create Voronoi regions with Geospatial data in Python

A step by step guide on creating and plotting Voronoi diagrams in Python

Photo by Andy Holmes on Unsplash

Assume you are planning to walk to a station to pick up a scooter. There are numerous stations available in nearby. Which one should you go to pick up your scooter and ride?

The closest station, right!

But, how do you know the closest station from your place?

Enter Voronoi diagrams.

A Voronoi diagram is a collection of polygons with all the points on a plane that is closest to the single object.

In other words, each polygon division correlates with a single object and contains all points which are closest to this single object.

Applications

Applications of Voronoi diagrams are many and often include determining which feature is closest to any given point. For example, to determine which school is nearest at a given point in a neighbourhood. Or which cell tower is closest to my cell phone to make a phone call.

The famous John Snow’s cholera outbreak mapping in 1854 that killed 500 people in five days is said to have plotted water pumps data on a chart and effectively constructing a Voronoi diagram for inner Neighbourhoods in London.

We Can not list out all possible use case for Voronoi diagrams here. Still, its applications are far-reaching including Anthropology and Archeology, Statistics and Data Analysis, to Marketing and Meteorology. Here is a long list of some of its applications.

Create Voronoi Regions with Python

In this tutorial, we create a Voronoi diagram using preschools dataset. The dataset contains all preschools in Uppsala county, Sweden. To create Voronoi diagrams, we are going to use several libraries, including Geopandas and Geovoronoi.

Let us import the libraries and read the data first.

import numpy as np
import geopandas as gpd
import contextily as ctx
import matplotlib.pyplot as plt
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
gdf = gpd.read_file("data/preschools.shp")
gdf.head()

As shown here in the first few rows of the data, we have a geometry column with Points and some other attributes like the county, post address, etc..

Photo by Andy Holmes on Unsplash

We can plot the points data. Just to give some context, we also read the boundary of the county and plot the points on top of it.

boundary = gpd.read_file(“data/uppsala.shp”)
fig, ax = plt.subplots(figsize=(12, 10))
boundary.plot(ax=ax, color=”gray”)
gdf.plot(ax=ax, markersize=3.5, color=”black”)
ax.axis(“off”)
plt.axis(‘equal’)
plt.show()

The output is a map showing all preschools (black dots) in Uppsala county, Sweden.

Photo by Andy Holmes on Unsplash

Before we calculate the Voronoi regions, we need to make sure two things. First, we need to check the projection of the data and then convert it to Web Mercator projection (epsg=339599.

boundary = boundary.to_crs(epsg=3395)
gdf_proj = gdf.to_crs(boundary.crs)

Second, we need to prepare the data to a format that Geovoronoi library can use. Here, we convert the boundary geometry into a union of the polygon. We also convert the Geopandas GeoSeries of Point objects to NumPy array of coordinates.

boundary_shape = cascaded_union(boundary.geometry)
coords = points_to_coords(gdf_proj.geometry)

Calculate Voronoi Regions

Now, that we have prepared the data, we can calculate Voronoi regions simply using Geovoronoi’s method voronoi_regions_from_coords().

# Calculate Voronoi Regions
poly_shapes, pts, poly_to_pt_assignments = voronoi_regions_from_coords(coords, boundary_shape)

The output holds the shapes, the points and identification link between the two.

Plot Voronoi Diagram

To plot Voronoi Diagrams, we also use the functionality of Geovoronoi — plot_voronoi_polys_with_points_in_area(). Here, we provide all the output from the above Voronoi calculations and the boundary shape.

fig, ax = subplot_for_map()
plot_voronoi_polys_with_points_in_area(ax, boundary_shape, poly_shapes, pts, poly_to_pt_assignments)
ax.set_title('Voronoi regions of Schools in Uppsala')
plt.tight_layout()
plt.show()

The output is a Voronoi diagram which delineates school points into regions that are closest to each location.

Photo by Andy Holmes on Unsplash

The above map is small, but the good thing is that we can tweak it using the familiar Matplotlib interface. Let us first see the parameters and documentation for the plot_voronoi_polys_with_points_in_area() function.

plot_voronoi_polys_with_points_in_area(
    ax,
    area_shape,
    poly_shapes,
    points,
    poly_to_pt_assignments=None,
    area_color='white',
    area_edgecolor='black',
    voronoi_and_points_cmap='tab20',
    voronoi_color=None,
    voronoi_edgecolor=None,
    points_color=None,
    points_markersize=5,
    points_marker='o',
    voronoi_labels=None,
    voronoi_label_fontsize=10,
    voronoi_label_color=None,
    point_labels=None,
    point_label_fontsize=7,
    point_label_color=None,
    plot_area_opts=None,
    plot_voronoi_opts=None,
    plot_points_opts=None,
)
Docstring:
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.

As you can see, there are a lot of parameters we can use to format the visualization. The documentation is also very clear and helpful.

We can increase the figure size and marker size of the points. We can also change the colour using Matplotlib Colormap. This is just one way, and there are other ways you can tweak which I leave it to you.

fig, ax = plt.subplots(figsize=(14,12))
plot_voronoi_polys_with_points_in_area(ax, boundary_shape, poly_shapes, pts, poly_to_pt_assignments,
 voronoi_and_points_cmap=’tab20c’,
 points_markersize=20)
ax.set_title(‘Upssalla Preschools — Voronoi Regions’)
ax.axis(“off”)
plt.tight_layout()
plt.show()

The output map is much better now compared to the previous with the defaults.

Photo by Andy Holmes on Unsplash

Conclusion

Voronoi Diagrams are useful and are used widely in many Geographic applications. In this tutorial, we have covered how to create Voronoi diagrams using Python.

The code and the data for this article are available in this Github repository.

Leave a Reply

Your email address will not be published. Required fields are marked *