Note
Measuring density is a typical exercise in urban analytics. momepy allows to measure different types (see API/Intensity); this notebook will outline the main principles.
momepy
[1]:
import momepy import geopandas as gpd import matplotlib.pyplot as plt
We will again use osmnx to get the data for our example and after preprocessing of building layer will generate tessellation layer.
osmnx
[2]:
import osmnx as ox point = (40.731603, -73.977857) dist = 1000 gdf = ox.geometries.geometries_from_point(point, dist=dist, tags={'building':True}) gdf_projected = ox.projection.project_gdf(gdf) gdf_projected = gdf_projected[gdf_projected.geom_type.isin(['Polygon', 'MultiPolygon'])] buildings = momepy.preprocess(gdf_projected, size=30, compactness=True, islands=True) buildings['uID'] = momepy.unique_id(buildings) limit = momepy.buffered_limit(buildings) tessellation = momepy.Tessellation(buildings, unique_id='uID', limit=limit).tessellation
Loop 1 out of 2.
Identifying changes: 100%|██████████| 3200/3200 [00:02<00:00, 1481.66it/s] Changing geometry: 100%|██████████| 2157/2157 [00:17<00:00, 120.11it/s]
Loop 2 out of 2.
Identifying changes: 100%|██████████| 1459/1459 [00:00<00:00, 2104.02it/s] Changing geometry: 100%|██████████| 863/863 [00:05<00:00, 171.39it/s]
Inward offset... Discretization...
9%|▊ | 71/833 [00:00<00:01, 706.33it/s]
Generating input point array...
100%|██████████| 833/833 [00:01<00:00, 587.68it/s]
Generating Voronoi diagram... Generating GeoDataFrame...
Vertices to Polygons: 100%|██████████| 356220/356220 [00:06<00:00, 51892.45it/s]
Dissolving Voronoi polygons...
100%|██████████| 3/3 [00:00<00:00, 716.73it/s]
Preparing limit for edge resolving... Building R-tree... Identifying edge cells... Cutting...
[3]:
f, ax = plt.subplots(figsize=(10, 10)) tessellation.plot(ax=ax) buildings.plot(ax=ax, color='white', alpha=.5) ax.set_axis_off() plt.show()
We have some edge effect here as we are using the buffer as a limit for tessellation in the middle of the urban fabric, but for these examples, we can work with it anyway. Keep in mind that values on the edge of this area will be skewed.
Covered area ratio, in our case measured on tessellation cells, requires GeoDataFrame containing spatial unit (cell, plot), and GeoDataFrame containing covering objects (buildings). On top of that, it currently requires passed areas for both gdfs and unique ID which links together spatial units and objects on them. We can either calculate areas before:
GeoDataFrame
[4]:
tessellation['area'] = momepy.Area(tessellation).series buildings['area'] = momepy.Area(buildings).series tess_car = momepy.AreaRatio(tessellation, buildings, 'area', 'area', 'uID') tessellation['CAR'] = tess_car.series
Or we can pass momepy.Area().series directly:
momepy.Area().series
[5]:
tess_car = momepy.AreaRatio(tessellation, buildings, momepy.Area(tessellation).series, momepy.Area(buildings).series, 'uID') tessellation['CAR'] = tess_car.series
[6]:
f, ax = plt.subplots(figsize=(10, 10)) tessellation.plot(ax=ax, column='CAR', legend=True, scheme='quantiles', k=10, cmap='viridis') buildings.plot(ax=ax, color='white', alpha=0.5) ax.set_axis_off() plt.show()
Because we know building heights for our buildings gdf, we can also calculate FAR. This part of New York has height data, only stored as strings, so we have to convert them to floats (or int) and fill NaN values with zero.
buildings
floats
int
NaN
FAR requires floor areas for building gdf instead of covered area.
[7]:
buildings['height'] = buildings['height'].fillna(0).astype(float) buildings['floor_area'] = momepy.FloorArea(buildings, 'height').series
[8]:
tessellation['FAR'] = momepy.AreaRatio(tessellation, buildings, 'area', 'floor_area', 'uID').series
[9]:
f, ax = plt.subplots(figsize=(10, 10)) tessellation.plot(ax=ax, column='FAR', legend=True, scheme='quantiles', k=10, cmap='viridis') buildings.plot(ax=ax, color='white', alpha=0.5) ax.set_axis_off() plt.show()
Location-based density is described in examples using spatial weights.