How to do an equal segment map (tessellation) analysis in Python?
This reading will summarize how to use Tesspy, the Python library used for geographical tessellation to analyze maps.
Sometimes while working on map analysis given a specific point of interest (POIs) (for instance, hospitals, restaurants, and schools), a common approach is to do a radius analysis with the POIs. This can be particularly helpful if we want to detect empty zones in the map of the POIs, or any zone with POIs that are too close to each other (see Figure 1).
A subtle distinction representing a completely different analysis is when we want to divide the maps into geometric shapes with no overlaps and gaps to compute how many POIs are in a segment (see Figure 2).
This type of segmentation is called in mathematics a tessellation, and there are regular and irregular tessellations. Recently the Journal of Open Source Software published an article about a new Python package for geographical tessellation that allows straightforward tessellation analysis over maps in only a few lines of code called TessPy.
This post presents an overview of the principal features that TessPy provides for regular and irregular tessellations and the type of exploratory data analysis we can perform.
Before continuing, we need to install the necessary libraries.
pip install tesspy, contextily, mapclassifyThese are the required packages for the analysis.
from tesspy import Tessellation
import matplotlib.pyplot as plt
import contextily as ctx
import geopandas as gpd
from shapely.geometry import Point
import warnings
warnings.simplefilter("ignore")For this display, we will be tessellating the map of Guadalajara, a city in the western part of Mexico.
The following code initializes a tessellation object of Guadalajara and the polygon object. We only need to type the street, city name, county, or region as it appears on Open Street Map.
guadalajara = Tessellation("Guadalajara")
guadalajara_polygon = guadalajara.get_polygon()Regular and irregular tessellations
Regular tessellations
Two common geometric shapes for tessellation are hexagons and squares. The next block creates and plots a tessellation with zoom level 8.
guadalajara_hexagons = guadalajara.hexagons(8)
fig, ax = plt.subplots(figsize=(10,15))
guadalajara_hexagons.to_crs("EPSG:3857").plot(ax=ax, facecolor="none", edgecolor="k", lw=1)
guadalajara_polygon.to_crs("EPSG:3857").plot(ax=ax, facecolor="none", edgecolor="k", lw=3)
ctx.add_basemap(ax=ax, source=ctx.providers.CartoDB.Voyager)
ax.set_axis_off()
ax.set_title("Guadalajara", fontsize=20)
In general, if we increase the zoom level, the size of the geometric shape will be reduced due to its reference to a special zoom of the Mercator projection.
guadalajara_hexagons = guadalajara.hexagons(9)
fig, ax = plt.subplots(figsize=(10,15))
guadalajara_hexagons.to_crs("EPSG:3857").plot(ax=ax, facecolor="none", edgecolor="k", lw=1)
guadalajara_polygon.to_crs("EPSG:3857").plot(ax=ax, facecolor="none", edgecolor="k", lw=3)
ctx.add_basemap(ax=ax, source=ctx.providers.CartoDB.Voyager)
ax.set_axis_off()
ax.set_title("Guadalajara", fontsize=20)
Increasing zoom levels allows for analyzing smaller map zones (see Figure 3 and Figure 4).
Similarly, we can produce a tessellation with squares with the following code block.
guadalajara_squares = guadalajara.squares(16)
fig, ax = plt.subplots(figsize=(10,15))
guadalajara_squares.to_crs("EPSG:3857").plot(ax=ax, facecolor="none", edgecolor="k", lw=1)
guadalajara_polygon.to_crs("EPSG:3857").plot(ax=ax, facecolor="none", edgecolor="k", lw=3)
ctx.add_basemap(ax=ax, source=ctx.providers.CartoDB.Voyager)
ax.set_axis_off()
ax.set_title("Guadalajara", fontsize=20)
plt.tight_layout()
Irregular tessellations
Another interesting feature that allows TessPy is the tessellation of a region using city blocks.
guadalajara_vor = guadalajara.city_blocks(n_polygons=100)
fig, ax = plt.subplots(figsize=(10,15))
guadalajara_vor.to_crs("EPSG:3857").plot(ax=ax, facecolor="none", edgecolor="k", lw=1)
guadalajara_polygon.to_crs("EPSG:3857").plot(ax=ax, facecolor="none", edgecolor="k", lw=3)
ctx.add_basemap(ax=ax, source=ctx.providers.CartoDB.Voyager)
ax.set_axis_off()
ax.set_title("Guadalajara", fontsize=20)
The n_polygons variable is used as the desired number of polygons, but it might not be achieved due to the size and number of blocks the city or region actually has.
Finally, the last tessellation method we will explore is the Voronoi Diagrams. In cases when we have a POI to study but we won’t be using a radius analysis, Voronoi Diagrams are helpful.
In the following block of code, we construct a tessellation with the POI key emergency, which includes, among others, ambulance stations, landing sites, and fire hydrants registered on Open Street Map.
guadalajara_vor = guadalajara.voronoi(n_polygons=16, poi_categories=["emergency"])
fig, ax = plt.subplots(figsize=(10,15))
guadalajara_vor.to_crs("EPSG:3857").plot(ax=ax, facecolor="none", edgecolor="k", lw=1)
guadalajara_polygon.to_crs("EPSG:3857").plot(ax=ax, facecolor="none", edgecolor="k", lw=3)
ctx.add_basemap(ax=ax, source=ctx.providers.CartoDB.Voyager)
ax.set_axis_off()
ax.set_title("Guadalajara", fontsize=20)
The “Map features” post of the Open Street Map’s wiki shows the complete POI list.
A different type of irregular tessellation available on TessPy is the Adaptive square method which you can explore in the official documentation.
A quick view of Exploratory Data Analysis over tessellations with TessPy
Let’s restart the exercise and suppose we want to analyze the quantity of POI in the keys shop, building, amenity, office, and public transport in Guadalajara City using Voronoi Tessellation. First, we initialize the object Tessellation and the polygon object.
guadalajara = Tessellation("Guadalajara")
guadalajara_polygon = guadalajara.get_polygon()The following line of code defines the map segmentation.
guadalajara_vor = guadalajara.voronoi(n_polygons=16, poi_categories=["shop", "building", "amenity", "office", "public_transport"])It’s possible to extract the POI data from the guadalajara object and compute how many POI are in each key.
poi_df = guadalajara.get_poi_data()
poi_df.head()
poi_df.iloc[:, -5:].sum()These are the results.
+------------------+------+
| amenity | 4699 |
| building | 2458 |
| office | 107 |
| public_transport | 3126 |
| shop | 1150 |
+------------------+------+We extract the amenity_data to perform further analysis.
poi_geodata = gpd.GeoDataFrame(
data=poi_df.drop(columns=["geometry"]),
geometry=poi_df[["center_longitude", "center_latitude"]]
.apply(Point, axis=1)
.values,
crs="EPSG:4326",
)
amenity_data = poi_geodata[poi_geodata["amenity"]].copy()Next, we plot the amenities on the map.
fig, ax = plt.subplots(1, 1, figsize=(15,15))
ax = amenity_data.to_crs("EPSG:3857").plot(ax=ax,color="tab:blue", markersize=0.5, alpha=0.5)
guadalajara_polygon.to_crs("EPSG:3857").plot(ax=ax,facecolor="none", edgecolor="tab:red", lw=1)
ctx.add_basemap(ax=ax,source=ctx.providers.CartoDB.Positron)
ax.set_axis_off()
ax.set_title("Guadalajara - Amenities", fontsize=12)
plt.show()
Finally, using the Voronoi tessellation, we perform a heatmap analysis of amenities data and plot the results into the map.
heat_map = gpd.sjoin(
guadalajara_vor, amenity_data, how="left", predicate="contains"
)
fig, ax = plt.subplots(1, 1, figsize=(15,15))
heat_map.plot(
ax=ax,
column="index_right",
lw=0.1,
alpha=1,
scheme="fisherjenks",
legend=False,
cmap="OrRd",
edgecolor="k"
)
ax.set_title("Heatmap of Guadalajara's Amenities", fontsize=12)
ax.set_axis_off()
plt.show()
Other exciting analyses can be performed with TessPy, like spatial autocorrelations and spatial lag. I encourage the reader to check the official TessPy documentation to find out all the features available in this package.
