Search
Examples of usage

Spatial weights are used across momepy. This notebook will illustrate its use on three examples.

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.

import osmnx as ox

gdf = ox.footprints.footprints_from_place(place='Kahla, Germany')
gdf_projected = ox.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, 4684.51it/s]
Changing geometry: 100%|██████████| 31/31 [00:00<00:00, 89.86it/s]
Loop 2 out of 2.
Identifying changes: 100%|██████████| 2520/2520 [00:00<00:00, 4689.22it/s]
Changing geometry: 100%|██████████| 2/2 [00:00<00:00, 60.83it/s]
Inward offset...
Discretization...
  1%|▏         | 36/2521 [00:00<00:06, 358.07it/s]
Generating input point array...
100%|██████████| 2521/2521 [00:03<00:00, 715.83it/s]
Generating Voronoi diagram...
Generating GeoDataFrame...
Vertices to Polygons: 100%|██████████| 267595/267595 [00:07<00:00, 35900.20it/s]
Dissolving Voronoi polygons...
Preparing limit for edge resolving...
Building R-tree...
  8%|▊         | 31/371 [00:00<00:01, 309.27it/s]
Identifying edge cells...
100%|██████████| 371/371 [00:01<00:00, 342.48it/s]
 22%|██▏       | 53/238 [00:00<00:00, 527.54it/s]
Cutting...
100%|██████████| 238/238 [00:00<00:00, 502.71it/s]

First order contiguity

Distance to neighbours

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.

sw1 = momepy.sw_high(k=1, gdf=tessellation, ids='uID')
buildings['neighbour_dist'] = momepy.NeighborDistance(buildings, sw1, 'uID').series
100%|██████████| 2518/2518 [00:04<00:00, 605.43it/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.

buildings = buildings.dropna(subset=['neighbour_dist'])
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()

Higher order / distance

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.

sw3 = momepy.sw_high(k=3, gdf=tessellation, ids='uID')

Average character

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.

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%|██████████| 2518/2518 [00:03<00:00, 730.57it/s]
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.

tessellation['mean_area_iq'] = momepy.AverageCharacter(
    tessellation, areas, sw3, 'uID', rng=(25, 75)).mean
100%|██████████| 2518/2518 [00:03<00:00, 663.96it/s]
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:

tessellation['med_area'] = momepy.AverageCharacter(
    tessellation, areas, sw3, 'uID', mode='median').median
100%|██████████| 2518/2518 [00:03<00:00, 774.48it/s]
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()

Weighted character

The weighted average is another example using the same spatial weights. For illustration, we can try area-weighted circular compactness:

circular_compactness = momepy.CircularCompactness(buildings)
buildings['weighted_circom'] = momepy.WeightedCharacter(
    buildings, circular_compactness.series, sw3, 'uID', momepy.Area(buildings).series).series
100%|██████████| 2518/2518 [00:03<00:00, 732.57it/s]
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()

Density

We will again use our Manhattan case study to illustrate Density.

point = (40.731603, -73.977857)
dist = 1000
gdf = ox.footprints.footprints_from_point(point=point, distance=dist)
gdf_projected = ox.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%|██████████| 3198/3198 [00:01<00:00, 3187.50it/s]
Changing geometry: 100%|██████████| 20/20 [00:00<00:00, 70.47it/s]
Loop 2 out of 2.
Identifying changes: 100%|██████████| 3165/3165 [00:00<00:00, 4689.01it/s]
Changing geometry: 100%|██████████| 2/2 [00:00<00:00, 54.70it/s]
Inward offset...
Discretization...
  1%|          | 33/3163 [00:00<00:09, 328.11it/s]
Generating input point array...
100%|██████████| 3163/3163 [00:06<00:00, 475.58it/s]
Generating Voronoi diagram...
Generating GeoDataFrame...
Vertices to Polygons: 100%|██████████| 496719/496719 [00:12<00:00, 39978.03it/s]
Dissolving Voronoi polygons...
Preparing limit for edge resolving...
Building R-tree...
 22%|██▏       | 36/163 [00:00<00:00, 358.46it/s]
Identifying edge cells...
100%|██████████| 163/163 [00:00<00:00, 357.33it/s]
100%|██████████| 95/95 [00:00<00:00, 540.44it/s]
/Users/martin/Strathcloud/Personal Folders/momepy/momepy/momepy/elements.py:434: UserWarning: Tessellation contains MultiPolygon elements. Initial objects should be edited. unique_id of affected elements: [830, 1976, 2867, 2872, 3129]
  "unique_id of affected elements: {}".format(list(uids))
Cutting...
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:

buildings['height'] = buildings['height'].fillna(0).astype(float)
buildings['floor_area'] = momepy.FloorArea(buildings, 'height').series

Now we merge floor areas to tessellation based on shared unique ID and generate spatial weights.

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.

gross = momepy.Density(
    tessellation, values='floor_area', spatial_weights=sw, unique_id='uID')
tessellation['gross_density'] = gross.series
100%|██████████| 3163/3163 [00:03<00:00, 801.99it/s]
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.

buildings['area'] = momepy.Area(buildings).series
tessellation = tessellation.merge(buildings[['uID', 'area']])
coverage = momepy.Density(
    tessellation, values='area', spatial_weights=sw, unique_id='uID')
tessellation['gross_coverage'] = coverage.series
100%|██████████| 3163/3163 [00:04<00:00, 748.59it/s]
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()