Street Population Density#
In the last two notebooks, we worked on how to determine the population of any given area using the raster dataset GHS-POP (Global Human Settlement Population Layer) for the P2023 release[1]. As our goal is to repeatedly determine the population density of Superblocks, we want to have an effective way to do this. As the Superblocks are made up of subgraphs of our road network, we showed a way to tessellate using roads.
Show code cell source
import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
import osmnx as ox
from networkx import ego_graph
from matplotlib._color_data import CSS4_COLORS
from rasterio.features import shapes
from rasterstats import zonal_stats
from scipy.spatial import Voronoi, voronoi_plot_2d
from scipy.odr import ODR, Model, RealData
from shapely import STRtree
from shapely.geometry import shape
from superblockify import get_edge_cells, get_ghsl, resample_load_window, \
add_edge_population, add_edge_cells, get_population_area
2024-09-21 05:45:10,503 | INFO | __init__.py:11 | superblockify version 1.0.0
Voronoi Tessellation#
We can use a Voronoi tessellation to divide up the area around a given node into cells. This looks like the following:
Show code cell source
# Generate Example Voronoi Tessellation
np.random.seed(1671593280)
points = np.random.multivariate_normal(mean=[0, 0], cov=[[1, 0], [0, 0.6]], size=40)
vor = Voronoi(points)
fig = voronoi_plot_2d(vor, figsize=(5, 4), show_vertices=False, line_colors='green',
line_width=1.5, line_alpha=0.6, point_size=4)
plt.tick_params(axis='both', which='both', bottom=False, top=False, labelbottom=False,
right=False, left=False, labelleft=False)
plt.show()
To keep the code effective, we renounce calculating the Voronoi tessellation for each partition of the road network. Instead, tessellate the whole road network once using roads and calculate the population and area of each road cell. Then, we can use this information to determine the population density of each Superblock repeatedly. More on the edge tesselation in the previous chapter Street Tessellation.
There are two approaches we want to consider integrating the population raster data over the road cells. The first approach is to use up-sampled raster data, as we did before. The second approach is to polygonize the raster data and intersect the road cells with it. For both, we need to consider used memory and computation time. Loading in a big city with an up-sampled raster dataset can take a lot of time and memory, but it needs to be up-sampled to reach better population approximations. Calculating intersections between polygons could be computationally heavy.
Data Preparation#
First, let’s get an example road network and tessellate it. We are Using La Crosse, Wisconsin as an example. By the Native American Neshnabé/Potawatomi, it is called ensewaksak. The working CRS is the World Mollweide projection, as it is an equal-area projection, and the GHS-POP raster data is in this projection. But the edge cells need to be calculated in the local CRS, to avoid distortion.
lc_graph = ox.graph_from_place("La Crosse, La Crosse County, Wisconsin, United "
"States",
network_type="drive")
lc_graph = ox.project_graph(lc_graph) # to local UTM zone
# Get Road Cells
lc_road_cells = get_edge_cells(lc_graph)
# Project to World Mollweide
lc_road_cells = lc_road_cells.to_crs("World Mollweide")
2024-09-21 05:45:15,974 | INFO | tessellation.py:101 | Calculating edge cells for graph with 4799 edges.
2024-09-21 05:45:21,915 | INFO | tessellation.py:155 | Tessellated 2606 edge cells in 0:00:05.
How does the tessellation look like?
# random color of CSS4_COLORS
lc_road_cells["color_name"] = np.random.choice(list(CSS4_COLORS.keys()),
size=len(lc_road_cells))
lc_road_cells.explore(color="color_name", style_kwds={"color": "black", "opacity": 0.5,
"weight": 1}, zoom_start=14)