Note
Spatial weights are used across momepy. This notebook will illustrate its use on three examples.
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 gdf = ox.geometries.geometries_from_place('Kahla, Germany', tags={'building':True}) gdf_projected = ox.projection.project_gdf(gdf) 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%|██████████| 2932/2932 [00:00<00:00, 7214.32it/s] Changing geometry: 100%|██████████| 630/630 [00:04<00:00, 148.53it/s]
Loop 2 out of 2.
Identifying changes: 100%|██████████| 2115/2115 [00:00<00:00, 10114.93it/s] Changing geometry: 100%|██████████| 172/172 [00:00<00:00, 185.18it/s]
Inward offset... Discretization...
6%|▌ | 114/2008 [00:00<00:01, 1133.81it/s]
Generating input point array...
100%|██████████| 2008/2008 [00:01<00:00, 1907.59it/s]
Generating Voronoi diagram... Generating GeoDataFrame...
Vertices to Polygons: 100%|██████████| 249298/249298 [00:04<00:00, 53635.00it/s]
Dissolving Voronoi polygons... Preparing limit for edge resolving... Building R-tree... Identifying edge cells... Cutting...
0it [00:00, ?it/s]
To calculate the mean distance to neighbouring buildings, we need queen contiguity weights of the first order capturing the relationship between immediate neighbours. Relationship between buildings is here represented by relationships between their tessellation cells.
[3]:
sw1 = momepy.sw_high(k=1, gdf=tessellation, ids='uID')
[4]:
buildings['neighbour_dist'] = momepy.NeighborDistance(buildings, sw1, 'uID').series
100%|██████████| 2005/2005 [00:01<00:00, 1761.93it/s]
Note: If there is no neighbour for a building denoted in spatial_weights, the value is set to np.nan. In GeoPandas older than 0.7.0, rows with NaN have to be removed before plotting with natural breaks scheme.
spatial_weights
np.nan
[5]:
buildings = buildings.dropna(subset=['neighbour_dist'])
[6]:
f, ax = plt.subplots(figsize=(10, 10)) buildings.plot(ax=ax, column='neighbour_dist', scheme='naturalbreaks', k=15, legend=True, cmap='Spectral') ax.set_axis_off() plt.show()
However, typical usage of spatial weights is to capture the vicinity of each feature. As illustrated in the previous notebook, there are multiple options on how to capture it. In this example, we will use queen contiguity of the higher order (3) based on morphological tessellation.
[7]:
sw3 = momepy.sw_high(k=3, gdf=tessellation, ids='uID')
Mean value of selected character within a vicinity of each cell (or building, plot) is a simple example. AverageCharacter can measure mean, median or mode and defaults to all. Each of them can be accessed using .mean, .median or .mode. .series will return mean.
AverageCharacter
all
.mean
.median
.mode
.series
[8]:
areas = momepy.Area(tessellation).series mean_area = momepy.AverageCharacter( tessellation, values=areas, spatial_weights=sw3, unique_id='uID') tessellation['mean_area'] = mean_area.mean
100%|██████████| 2005/2005 [00:01<00:00, 1901.77it/s]
[9]:
f, ax = plt.subplots(figsize=(10, 10)) tessellation.plot(ax=ax, column='mean_area', legend=True, scheme='quantiles', k=15, cmap='Blues_r') buildings.plot(ax=ax, color="white", alpha=0.4) ax.set_axis_off() plt.show()
In some cases, we might want to eliminate the effect of outliers. To do so, we can specify the range on which should AverageCharacter calculate mean. Below we will measure only interquartile mean.
[10]:
tessellation['mean_area_iq'] = momepy.AverageCharacter( tessellation, areas, sw3, 'uID', rng=(25, 75)).mean
100%|██████████| 2005/2005 [00:01<00:00, 1603.13it/s]
[11]:
f, ax = plt.subplots(figsize=(10, 10)) tessellation.plot(ax=ax, column='mean_area_iq', legend=True, scheme='quantiles', k=15, cmap='Greens_r') buildings.plot(ax=ax, color="white", alpha=0.4) ax.set_axis_off() plt.show()
Another option would be to calculate median only:
[12]:
tessellation['med_area'] = momepy.AverageCharacter( tessellation, areas, sw3, 'uID', mode='median').median
100%|██████████| 2005/2005 [00:00<00:00, 3063.14it/s]
[13]:
f, ax = plt.subplots(figsize=(10, 10)) tessellation.plot(ax=ax, column='med_area', legend=True, scheme='quantiles', k=15, cmap='Purples_r') buildings.plot(ax=ax, color="white", alpha=0.4) ax.set_axis_off() plt.show()
The weighted average is another example using the same spatial weights. For illustration, we can try area-weighted circular compactness:
[14]:
circular_compactness = momepy.CircularCompactness(buildings) buildings['weighted_circom'] = momepy.WeightedCharacter( buildings, circular_compactness.series, sw3, 'uID', momepy.Area(buildings).series).series
100%|██████████| 2003/2003 [00:00<00:00, 2473.83it/s]
[15]:
f, ax = plt.subplots(figsize=(10, 10)) buildings.plot(ax=ax, column='weighted_circom', legend=True, scheme='quantiles', k=15, cmap='viridis') ax.set_axis_off() plt.show()
We will again use our Manhattan case study to illustrate Density.
Density
[16]:
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
Identifying changes: 100%|██████████| 3200/3200 [00:01<00:00, 2176.30it/s] Changing geometry: 100%|██████████| 2157/2157 [00:18<00:00, 115.38it/s]
Identifying changes: 100%|██████████| 1459/1459 [00:00<00:00, 2216.25it/s] Changing geometry: 100%|██████████| 863/863 [00:05<00:00, 164.14it/s]
8%|▊ | 70/833 [00:00<00:01, 692.05it/s]
100%|██████████| 833/833 [00:01<00:00, 583.48it/s]
Vertices to Polygons: 100%|██████████| 356220/356220 [00:06<00:00, 54018.76it/s]
Dissolving Voronoi polygons...
100%|██████████| 3/3 [00:00<00:00, 681.04it/s]
Preparing limit for edge resolving... Building R-tree... Identifying edge cells... Cutting...
[17]:
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()
To get gross density, we need to know floor areas:
[18]:
def clean_heights(x): try: return float(x) except ValueError: return 0 buildings['height'] = buildings['height'].fillna(0).apply(clean_heights) buildings['floor_area'] = momepy.FloorArea(buildings, 'height').series
Now we merge floor areas to tessellation based on shared unique ID and generate spatial weights.
[19]:
tessellation = tessellation.merge(buildings[['uID', 'floor_area']]) sw = momepy.sw_high(k=3, gdf=tessellation, ids='uID')
Density is then following the same principle as illustrated above.
[20]:
gross = momepy.Density( tessellation, values='floor_area', spatial_weights=sw, unique_id='uID') tessellation['gross_density'] = gross.series
100%|██████████| 833/833 [00:00<00:00, 1148.61it/s]
[21]:
f, ax = plt.subplots(figsize=(10, 10)) tessellation.plot(ax=ax, column='gross_density', legend=True, scheme='naturalbreaks', k=15) buildings.plot(ax=ax, color='white', alpha=.5) ax.set_axis_off() plt.show()
In a similar way can be done gross coverage.
[22]:
buildings['area'] = momepy.Area(buildings).series tessellation = tessellation.merge(buildings[['uID', 'area']])
[23]:
coverage = momepy.Density( tessellation, values='area', spatial_weights=sw, unique_id='uID') tessellation['gross_coverage'] = coverage.series
100%|██████████| 833/833 [00:00<00:00, 1022.98it/s]
[24]:
f, ax = plt.subplots(figsize=(10, 10)) tessellation.plot(ax=ax, column='gross_coverage', legend=True) buildings.plot(ax=ax, color='white', alpha=.5) ax.set_axis_off() plt.show()