Voronoi Grids: A Practical Application

Quick Success Data Science

Voronoi grids, also called Voronoi diagrams, are used to partition a plane into discrete regions around a given set of seed points. For each seed, there is a corresponding region, called a Voronoi cell, within which all points on the plane are closer to that seed than to any other.

Voronoi diagrams have applications in many fields, including computer science, geography, biology, and urban planning. A particularly important application is to map the nearest airfields for planes needing to make an emergency landing.

The government of Melbourne, Australia uses this tool to produce school catchment maps. "Catchment" refers to those students who reside in a particular area and are guaranteed a position in specific schools. Because students are eligible to attend the nearest primary or high school to where they live – as measured by Euclidian distance – the map of school zones is by default a Voronoi diagram.

In this Quick Success Data Science project, we'll explore the concept of Voronoi diagrams by making our own version of the Melbourne catchment map. We'll use a subset of primary schools in the metropolitan area and grid them with the SciPy library's Voronoi class. We'll then use the Folium library to overlay the Voronoi diagram on a street map of Melbourne.

The Dataset

To produce the dataset, I used the Victoria government's Find My School website to look up the addresses for 109 primary schools in the Melbourne metropolitan area. I then used LatLong.net to convert the addresses to decimal latitude and longitude and stored the results in a CSV file in this Gist.

The SciPy Voronoi Implementation

Python's SciPy scientific library is designed for mathematics, science, and engineering, and addresses many standard problem domains in scientific computing. It's built on and supplements Numerical Python (NumPy) and provides many user-friendly and efficient numerical routines.

For making Voronoi diagrams, SciPy provides the scipy.spatial.voronoi() class, which uses the Qhull library to compute the Voronoi grid. As mentioned previously, all the locations within a grid cell should be closer to the seed point used to generate that cell than to any other seed point.

In the diagram below, taken from the SciPy docs, the seed points are arranged orthogonally, leading to a straightforward grid pattern with the seed points in the centers of the cells.

The seed points are shown in blue. The vertices, where lines in the diagram meet, are colored in orange. The ridges (lines) between the seed points are drawn in black. These form the boundaries of the Voronoi cells.

Because most of the seed points in the previous diagram are along the edge of the map, only the center seed point is associated with a finite cell or region, and thus is surrounded on all sides by finite ridges. The other ridges are dashed as they extend to infinity and thus bound infinite cells that never "close," as there are no outer seed points from which to compute an intervening ridge.

The diagram below, also from the docs, was built from a non-orthogonal arrangement of seed points and consequently has more complex cells. Note how some of the ridges are drawn with solid lines and some with dashed.

The ridges drawn with solid lines converge and meet somewhere, creating a closed and finite region. The dashed lines are divergent and never intersect, so their cells have infinite area. In practice, these infinite cells are generally disregarded or artificially bounded, such as with a map-edge polygon.

We can access lists of the regions (cells), and the coordinates used to define the regions, with attributes such as voronoi.regions, voronoi.vertices, voronoi.ridge_vertices, and voronoi.ridge_points. We'll use these later to create plottable polygons for the Voronoi cells.

---

The Code

The code below, written in JupyterLab, follows these steps:

1. Load the CSV file of school names and locations.
2. Use the school locations to create a Voronoi grid.
3. Transform the Voronoi grid into plottable polygons.
4. Truncate the edge cells with a square bounding box.
5. Plot the schools and Voronoi grid on a street map of Melbourne.

Importing Libraries and Preparing the Data

To accomplish our task we'll need pandas, geopandas, shapely (bundled with geopandas), Folium, and SciPy. The pandas and geopandas libraries permit the loading and processing of the data, shapely creates plottable polygons, Folium permits mapping, and SciPy provides the algorithm for building the Voronoi grid.

If you're not familiar with it, GeoPandas is an open-source, third-party library designed to support Geospatial mapping in Python. It extends the datatypes used by the pandas library and makes working with geospatial vector data similar to working with tabular data. It also enables operations in Python that would otherwise require a dedicated geospatial database, such as Post GIS.

import pandas as pd
import geopandas as gpd
import folium
from shapely.geometry import Polygon
from Scipy.spatial import Voronoi

# Load the school locations CSV file into a DataFrame:
df = pd.read_csv('https://bit.ly/3MYYegT')

# Create a GeoDataFrame with Point geometries:
gdf = gpd.GeoDataFrame(df, 
                       geometry=gpd.points_from_xy(df['Longitude'], 
                                                   df['Latitude']), 
                       crs='EPSG:4326')

After making the imports, we use pandas to read the CSV file and geopandas to turn the schools' coordinates into point geometries.

A GeoDataFrame is a pandas DataFrame with a special "geometry" column for location data. This column bundles together the type of geometric object (such as a point, line string, polygon, etc.) and the coordinates (longitude and latitude) needed to draw it.

Notice that when making the GeoDataFrame we have to provide a coordinate reference system using the argument crs=EPSG:4326. "EPSG" stands for "European Petroleum Survey Group" dataset and the 4326 code is used by the Global Positioning System (GPS). This system projects the latitude and longitude coordinates from our 3D planet onto the flat surface of a map. You can read more about it here.

Creating the Voronoi Grid

The next step is to create the Voronoi grid by passing the GeoDataFrame's "Longitude" and "Latitude" columns to the SciPy library's Voronoi() class.

# Create a Voronoi diagram using the GeoDataFrame.
# Use -1 to check for and exclude regions that extend to infinity:
vor = Voronoi(gdf[['Longitude', 'Latitude']])
voronoi_polygons = [Polygon(vor.vertices[region]) 
                    for region in vor.regions 
                    if region and -1 not in region]

# Create a GeoDataFrame with the Voronoi polygons:
gdf_voronoi = gpd.GeoDataFrame(geometry=voronoi_polygons, 
                               crs='EPSG:4326')

Next, we make a list of the Voronoi polygons and remove empty regions and regions that extend to infinity. This requires using the regions attribute of the vor object returned by the Voronoi() class. If this attribute starts with -1, then we know that the region never closes. As part of the process, we pass the vertices attribute to shapely's Polygon() class, which generates plottable polygons that we add to a new GeoDataFrame.

Truncating the Voronoi Cells

Because we're not using all the primary schools in the area, many edge cells will be unconstrained and may extend for ridiculously long distances. To handle these, we'll create a bounding box and use it to truncate any Voronoi cells (polygons) that extend beyond it.

# Define the bounding box lat-lon limits:
max_lat, min_lat, max_lon, min_lon = (-37.75, -37.9, 
                                      145.18, 144.84)

# Create the bounding box as a Shapely Polygon
bounding_box = Polygon.from_bounds(min_lon, min_lat, 
                                   max_lon, max_lat)

# Truncate each Voronoi polygon with the bounding box:
truncated_polygons = [polygon.intersection(bounding_box) for 
                      polygon in gdf_voronoi.geometry]

# Create a GeoDataFrame with the truncated polygons:
gdf_truncated = gpd.GeoDataFrame(geometry=truncated_polygons, 
                                 crs='EPSG:4326')

We first define the lat-lon limits of the bounding box, using values slightly larger than those found in the CSV file. We then call Polygon() again to create the plottable polygon, then use list comprehension and the shapely Polygon intersection() method to truncate each polygon in the GeoDataFrame. Note how we use the special geometry column (gdf_voronoi.geometry) for this. We finish by creating a new GeoDataFrame to hold the truncated polygons.

Plotting the Map

The Voronoi diagram is useless without a way to relate it to the real world. So, we'll use Folium to draw it over a street map of the city (by choosing OpenStreetMap as the tiles argument).

We'll also need the school locations and related information. This requires looping through the GeoDataFrame and adding each marker to the map. We can control the icon inside the marker by setting the icon argument to the "home" symbol. You can find a list of icons like this on the Glyphicons page. You can also create custom icons, as described in this article.

# Create a Folium map centered on the average coordinates of the schools:
map_center = [gdf['Latitude'].mean(), gdf['Longitude'].mean()]
school_map = folium.Map(location=map_center, 
                        zoom_start=12, 
                        tiles='OpenStreetMap')

# Plot the truncated Voronoi polygons on the map:
folium.GeoJson(gdf_truncated).add_to(school_map)

# Add markers for each school:
for index, school in gdf.iterrows():
    folium.Marker(
        location=[school['Latitude'], school['Longitude']],
        popup=f"{school['School']}n
        {school['Street Address']}n{school['Town']}",
        icon=folium.Icon(color='blue', icon='home')
    ).add_to(school_map)

# Save the map as an HTML file (optional):
# school_map.save('school_voronoi_map_truncated.html')

# Display the map in the notebook:
school_map

Here's the resulting map, which permits zooming and panning.

Some schools at the edge of the diagram are not in colored cells as their "ridges" extend to infinity without converging. These can be "fixed" by merging the Voronoi cell with another polygon. We'll discuss this option shortly.

Clicking on the markers launches a pop-up window with the school's name and address.

To turn off the blue fill color for the polygons, just modify the style_function parameter when creating the folium.GeoJson object, as shown below.

folium.GeoJson(gdf_truncated, 
               style_function=lambda x: {'fillColor': 'none'}).add_to(school_map)

Now all you'll see are the cell boundaries.

Incorporating a Municipal Polygon

You can view the actual Melbourne catchment map at the beginning of this article, at the Find My School site (click on a school icon to activate the grid), or at mangomap.com.

Unlike our map, the real school zones blend the Voronoi polygons with other boundaries, such as the coastline along Port Philip Bay. Incorporating all these additional boundaries is beyond the scope of this article, but we can see how it's done by using the City of Melbourne's municipal boundary polygon.

The city provides this polygon as a shapefile. A shapefile is a geospatial vector data format for geographic information system (GIS) software. Despite its name, a "shapefile" isn't a single file but a collection of files in a single folder. You can read more about them here.

To download the shapefile, first navigate to the Municipal Boundary page of the City of Melbourne's Open Data site and then export the shapefile (as "whole dataset"). This data is shared under a CC Attribution 4.0 International license.

Store the zipped folder in the same folder as your Python script or notebook, so you can use the code provided below. Do not unzip the folder.

Here's the code for truncating and merging the existing polygons with the municipal boundary. It's similar to the truncation process used for the bounding box and assumes that you've already run the previous code (for those running the code in a notebook).

# Read in shapefile of Melbourne City Limits as GeoDataFrame:
city_limits = gpd.read_file('municipal-boundary.zip')

# Truncate Voronoi polygons by the city_limits polygon:
truncated_polygons = [polygon.intersection(city_limits.geometry.iloc[0]) 
                      for polygon in gdf_voronoi.geometry]

# Create a GeoDataFrame with the truncated polygons:
gdf_truncated = gpd.GeoDataFrame(geometry=truncated_polygons, 
                                 crs='EPSG:4326')

# Create a Folium map centered on the average coordinates of the schools:
map_center = [gdf['Latitude'].mean(), gdf['Longitude'].mean()]
school_map = folium.Map(location=map_center, 
                        zoom_start=12, 
                        tiles='OpenStreetMap')

# Plot truncated Voronoi polygons on the map:
folium.GeoJson(gdf_truncated).add_to(school_map)

# Add markers for each school:
for index, school in gdf.iterrows():
    folium.Marker(
        location=[school['Latitude'], school['Longitude']],
        popup=f"{school['School']}n
        {school['Street Address']}n{school['Town']}",
        icon=folium.Icon(color='blue', icon='home')
    ).add_to(school_map)

# Display map in notebook:
school_map

The new map, shown below, contains only the Voronoi cells within the Melbourne city limits. Many cell boundaries now conform to non-Voronoi-related features such as the Yarra River, the docks, and major roads.

Some cells are isolated from their respective schools. This is because the Melbourne municipal polygon isn't actually used to define the school zones. We're using it here to demonstrate how the final map could be created once you have a collection of all the correct polygons. As you saw, it's a simple process.

---

Summary

A Voronoi diagram is a nearest-neighbor mapping that lets you subdivide a plane with n generating points into n convex polygons. Each polygon contains one generating point and every location in a given polygon is closer to its generating point than to any other generating point.

In this example, we used a Voronoi diagram to partition a map of Melbourne, Australia into school catchment areas. Students in each zone are closer to the associated primary school than to any other primary school in the city. This map differs from the real thing as we didn't use all the schools or incorporate all the modifying polygons, such as for rivers and coastlines.

Python's SciPy library contains built-in functionality for generating Voronoi grids. With the help of additional third-party libraries like geopandas and Folium, you can project these grids onto maps for real-world applications.

---

Further Reading

If you want to learn more about Voronoi Diagrams, check out this informative article by Francesco Bellilli.

The fascinating world of Voronoi diagrams

---

Thanks!

Thanks for reading and please follow me for more Quick Success Data Science projects in the future.

Tags: Geospatial Hands On Tutorials Python Programming Scipy Voronoi Diagrams

Comment