Note

This page was generated from user_guide/weights/examples.ipynb.
Interactive online version: Binder badge

Examples of usage

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

[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.

[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]

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.

[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.

[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()
../../_images/user_guide_weights_examples_9_0.png

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.

[7]:
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.

[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()
../../_images/user_guide_weights_examples_14_0.png

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()
../../_images/user_guide_weights_examples_17_0.png

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()
../../_images/user_guide_weights_examples_20_0.png

Weighted character

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()
../../_images/user_guide_weights_examples_23_0.png

Density

We will again use our Manhattan case study to illustrate 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
Loop 1 out of 2.
Identifying changes: 100%|██████████| 3200/3200 [00:01<00:00, 2176.30it/s]
Changing geometry: 100%|██████████| 2157/2157 [00:18<00:00, 115.38it/s]
Loop 2 out of 2.
Identifying changes: 100%|██████████| 1459/1459 [00:00<00:00, 2216.25it/s]
Changing geometry: 100%|██████████| 863/863 [00:05<00:00, 164.14it/s]
Inward offset...
Discretization...
  8%|▊         | 70/833 [00:00<00:01, 692.05it/s]
Generating input point array...
100%|██████████| 833/833 [00:01<00:00, 583.48it/s]
Generating Voronoi diagram...
Generating GeoDataFrame...
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()

../../_images/user_guide_weights_examples_26_0.png

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()
../../_images/user_guide_weights_examples_33_0.png

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()
../../_images/user_guide_weights_examples_37_0.png