Spatial Data Science with CityJSON¶

The purpose of this Notebook is to *work with* the product of osm_LoD1_3DCityModel; a previously created CityJSON city model.

This notebook will:

1. allow the user to execute an application of Spatial Data Science

a) population estimation --with a previous census metric population growth rate and projected (future) population are also possible and
b) a measure of Building Volume per Capita

2. produce an interactive visualization -via pydeck- which a user can navigate, query and share that;

a) colour buildings by type (to easily visualize building stock)

3. propose several Geography and Sustainable Development Education conversation starters for Secondary and Tertiary level students

Please Note:

*The village* processing option is meant for areas with no more than for 2 500 buildings.

In [1]:

#- load the magic

%matplotlib inline
import os
from pathlib import Path

import numpy as np
import pandas as pd
import geopandas as gpd
import shapely
from shapely.geometry import Polygon, shape, mapping
import json
import geojson

from cjio import cityjson

import matplotlib.pyplot as plt
import pydeck as pdk

In [2]:

import warnings
warnings.filterwarnings('ignore')

The area under investigation is Mamre, Cape Town. South Africa.

In [3]:

#- change to harvest the appropriate CityJSON

#jparams = json.load(open('osm3DUEstate_param.json'))
#jparams = json.load(open('osm3DwEstate_param.json'))
#jparams = json.load(open('osm3DCPUT_param.json'))
jparams = json.load(open('osm3DMamre_param.json'))
#jparams = json.load(open('osm3Dobs_param.json'))
#jparams = json.load(open('osm3DsRiver_param.json'))

In [4]:

cm = cityjson.load(path=jparams['cjsn_solid']) #-- citjsnClean_uEstate10m.json in the result folder

In [5]:

print(cm)

CityJSON version = 1.1
EPSG = 32734
bbox = [263698.0877859225, 6287743.433441361, 145.1199951171875, 6287743.433441361, 6287743.433441361, 252.57000732421875]
=== CityObjects ===
|-- TINRelief (1)
|-- Building (2187)
===================
materials = False
textures = False

In [6]:

df = cm.to_dataframe()
#- remove the first feature: the terrain
df = df[1:]    

#- harvest the crs
theinfo = cm.get_info()
crs = theinfo[1]

gdf = gpd.GeoDataFrame(df, geometry=[shape(d) for d in df.pop("footprint")], crs=crs[7:])#jparams['crs'])
#gdf.head(2)

1. Spatial Data Science¶

We start with basic spatial analysis

We'll estimate the population, within our area of interest, and then
calculate the Building Volume Per Capita (BVPC).

While estimating population is well documented; recent investigations to understand overcrowding have led to newer measurements.

The most noteable of these is Building Volume Per Capita (BVPC) (Ghosh, T; et al. 2020). BVPC is the cubic meters of building per person. BVPC tells us how much space one person has per residential living unit (a house / apartment / etc.). It is *a proxy measure of economic inequality and a direct measure of housing inequality*.

BVPC builds on the work of (Reddy, A and Leslie, T.F., 2013) and attempts to integrate with several Sustainable Development Goals (most noteably: SDG 11: Developing sustainable cities and communities) and captures the average *'living space'* each person has in their home.

These analysis expect the user to have some basic knowledge about the environment under inquiry / investigation

In [7]:

gdf.head(2)

Out[7]:

	osm_id	building	building:levels	plus_code	ground_height	building_height	roof_height	address	amenity	residential	bottom_roof_height	building:use	geometry
328118446	328118446.0	yes	1	CMXPX3XW+XQR	179.43	4.1	183.53	NaN	NaN	NaN	NaN	NaN	POLYGON ((265041.635 6289775.059, 265050.531 6...
328118447	328118447.0	church	2	CPXWXXX3+XFX	179.04	6.9	185.94	Kerk Street Mamre 7347 Cape Town	place_of_worship	NaN	NaN	NaN	POLYGON ((265070.977 6289761.325, 265072.315 6...

In [8]:

#gdf.plot()
# have a look at the building type and amenities available
gdf['building'].unique()

Out[8]:

array(['yes', 'church', 'house', 'public', 'civic', 'office', 'retail',
       'clinic', 'school', 'cabin', 'garage', 'greenhouse', 'roof',
       'kindergarten', 'clubhouse', 'service', 'detached', 'shed'],
      dtype=object)

1. a) Estimate Population:

(with population growth rate and population projection possible too)

In [9]:

#--we only want building=house or =apartment or =residential
gdf_pop = gdf[gdf["building"].isin(['house', 'semidetached_house', 'terrace', 'apartments', 'residential', 'dormitory', 'cabin'])].copy()

In [10]:

#- some data wrangling to replace 'bld:residential' to 'bld:student' if 'residential:student'
gdf2 = gdf_pop.copy()

if 'residential' in gdf2.columns:
    df_res = gdf2[gdf2['residential'] == 'student']
    #df_res = df2[df2['building:use'] != None]
    df_res = df_res[~df_res['residential'].isna()]
    gdf_pop.loc[df_res.index, 'building'] = df_res['residential'] 

#- some more data wrangling
with pd.option_context("future.no_silent_downcasting", True):
    if 'building:flats' in gdf_pop.columns: 
        gdf_pop['building:flats'] = pd.to_numeric(gdf_pop['building:flats'].fillna(0).infer_objects(copy=False))
    if 'building:units' in gdf_pop.columns:    
        gdf_pop['building:units'] = pd.to_numeric(gdf_pop['building:units'].fillna(0).infer_objects(copy=False))
    if 'beds' in gdf_pop.columns:   
        gdf_pop['beds'] = pd.to_numeric(gdf_pop['beds'].fillna(0).infer_objects(copy=False))
    if 'rooms' in gdf_pop.columns:   
        gdf_pop['rooms'] = pd.to_numeric(gdf_pop['rooms'].fillna(0).infer_objects(copy=False))

gdf_pop["building:levels"] = pd.to_numeric(gdf_pop["building:levels"])

gdf_pop.head(2)

Out[10]:

	osm_id	building	building:levels	plus_code	ground_height	building_height	roof_height	address	amenity	residential	bottom_roof_height	building:use	geometry
656840974	656840974.0	house	1	C6XCX2X5+XGV	170.33	4.1	174.43	39 Dove Lane Mamre 7347 Cape Town	NaN	NaN	NaN	NaN	POLYGON ((264864.982 6288741.008, 264864.617 6...
656840975	656840975.0	house	1	C6XCXCXJ+XQX	169.70	4.1	173.80	37 Dove Lane Mamre 7347 Cape Town	NaN	NaN	NaN	NaN	POLYGON ((264865.350 6288749.141, 264865.026 6...

In [11]:

gdf_pop['building'].value_counts()

Out[11]:

building
house    1628
cabin     242
Name: count, dtype: int64

This area is urban with single level housing units. To estimate population is thus pretty straight forward.

We start with local knowledge.

On average there are roughly 6 people per building:house in this area.

An *informal* structure (shack) is tagged building:cabin and houses 4 people.

Your Participation!

We will execute the calculation programmatically. Fill in the relevant variables in the cell below

In [12]:

#- average number of residents per formal house
f_house = 6
#- average number of residents per informal structure
inf_structure = 4

Furthermore:
- social housing is tagged building:residential with the number of occupants iether the number of informal structure occupants or building:flats * inf_structure
- A social_facility (carehome, shelter, etc.) harvests the beds 'key:value' pair.
- building:apartment harvests the building:flats 'key:value' pair (the number of units) to calculate *3 people per apartment.
- *Student accomodation*:

University owed: is tagged building:dormitory with residential:university and harvests the beds 'key:value' pair.

Private for-profit: is tagged building:residential or :dormitory with residential:student and then harvests the building:flats or :rooms 'key:value' pair (the number of units) to calculate *1 people per apartment; if level: > 1 else *3 people in a house share.

The tagging scheme and numbers is based on how your community is mapped and local knowledge

In [13]:

c = gdf_pop.columns

def pop(row):
    #- formal house
    if row['building'] == 'house' or row['building'] == 'semidetached_house':
        return f_house
    if row['building'] == 'terrace' and 'building:units' in c:
        return row['building:units'] * f_house

    #- informal structure (shack)
    if row['building'] == 'cabin':
        return inf_structure

    #- in this case social housing
    if row['building'] == 'residential' and 'social_facility' in c and row['social_facility'] is np.nan:
        if row['building:levels'] > 1:
            if 'rooms' in row and row['rooms'] != 0:
                return row['rooms']
            if 'building:flats' in row and row['building:flats'] != 0:
                return row['building:flats'] * inf_structure
        else:
            return inf_structure

    #-- social facility [shelter / carehome]
    if row['building'] == 'residential' and 'social_facility' in c and row['social_facility'] is not np.nan:
        if row['building:levels'] > 1:
            if 'building:units' in row and row['building:units'] != 0:
                return row['building:units'] 
            if 'beds' in row and row['beds'] != 0:
                return row['beds']
        else:
            return inf_structure

    #- formal apartment
    if row['building'] == 'apartments':
        if 'rooms' in row and row['rooms'] != 0:
            return row['rooms']
        else:
            return row['building:flats'] * 3
        
    #- private student residence 
    if row['building'] == 'student':
        if row['building:levels'] > 1:
            if 'rooms' in row and row['rooms'] != 0:
                return row['rooms']
            else:
                return row['building:flats']
        else:
            return 3
            
    # university owned student residence
    if row['building'] == 'dormitory' and row['residential'] == 'university':
        if row['building:levels'] > 1:
            if 'rooms' in row and row['rooms'] != 0:
                return row['rooms']
            if 'beds' in row and row['beds'] != 0:
                return row['beds']
        else:
            return 3

gdf_pop['pop'] = gdf_pop.apply(lambda x: pop(x), axis=1)

est_pop = gdf_pop['pop'].sum()
print('The estimated population is:', est_pop)

The estimated population is: 10736

The official STATSSA 2011 census figure, for this community, is 9048.

We can calculate the annual population growth rate using the formula for Annual population growth:

$$r = \frac{\ln{[\frac{End Population}{Start Population}}]}{n} * 100 = \frac{\ln{[\frac{10 736}{9048}}]}{12} * 100 = 1.43\%$$

Your Participation!

It is possible to execute the calculation programmatically. Fill in the relevant variables in the cell below

In [14]:

#- previous population
start_population = 9048

#- period in years from the previous census
years = 12

In [15]:

#-execute
r = (np.log(est_pop/start_population)/years) * 100
print('population growth rate of approximately:', round(r, 2), '%')

population growth rate of approximately: 1.43 %

To conclude; we can project into the future with a very basic formula to estimate the population x-years from now:

$$p = P_o * (1 + r)^{t} = p = 10736 * (1 + 0.0143)^{10} = 12 368$$

Your Participation!

It is possible to execute the calculation programmatically. Fill in the variables in the cell below

In [16]:

#- period in years from now
years = 10

In [17]:

p = est_pop * (1 + (r/100))**years

print('estimated population', years ,'years from now:', int(p))

estimated population 10 years from now: 12368

1. b) Building Volume Per Capita (BVPC): BVPC = total population of a community divided by sum of building volume

In [18]:

#gdf_pop.head(3)

In [19]:

gdf_pop['area'] = gdf_pop['geometry'].area#\.map(lambda p: p.area)
gdf_pop['volume'] = gdf_pop['area'] * gdf_pop['building_height']

#- remove the volume of the ground floor (unoccupied) when building:levels > 7 [this is an arbitrary number based on local knowledge]
#- typically the space is reserved for some other function: retail, etc. 
gdf_pop['volume'] = [
    (row['volume'] - row['area'] * 2.8) if (
        'social_facility' in row and row['social_facility'] is np.nan and row['building:levels'] > 7 and
        row['building'] in ['residential', 'apartments', 'student']
    ) else row['volume']
    for _, row in gdf_pop.iterrows()
]

gdf_pop['bvpc'] =  gdf_pop['volume'] / gdf_pop['pop']

gdf_pop.tail(2)

Out[19]:

	osm_id	building	building:levels	plus_code	ground_height	building_height	roof_height	address	amenity	residential	bottom_roof_height	building:use	geometry	pop	area	volume	bvpc
12289266	12289266.0	house	1	CCX5X6XM+X5X	185.75	4.1	189.85	22 Clarkeson Street Mamre 7347 Cape Town	NaN	NaN	NaN	NaN	POLYGON ((265305.077 6288750.422, 265302.872 6...	6	344.679759	1413.187012	235.531169
12357148	12357148.0	house	1	C8XXXCX5+XQR	193.26	4.1	197.36	2 Tol Street Mamre 7347 Cape Town	NaN	NaN	NaN	NaN	POLYGON ((266000.758 6288992.653, 265989.649 6...	6	329.618113	1351.434264	225.239044

In [20]:

print(gdf_pop['bvpc'].describe())

count    1870.000000
mean       77.901895
std        50.554879
min        11.821442
25%        37.163535
50%        68.396996
75%       104.641182
max       400.313316
Name: bvpc, dtype: float64

In [21]:

bvpc = round(gdf_pop['volume'].sum() / est_pop, 3)

print('Building Volume Per Capita (BVPC):', bvpc)

Building Volume Per Capita (BVPC): 79.495

This BVPC value is general.

We can seperate building:house from building:cabin and building:residential to undertand the differences between *formal and informal* housing in this area.

We want to understand the living space (the cubic-meter BVPC value) each person has in thier home

In [22]:

formal = gdf_pop[gdf_pop["building"].isin(['house', 'semidetached_house', 'terrace', 'apartment'])].copy()
f_pop = formal['pop'].sum()
#f_area = formal['area'].mean()

informal = gdf_pop[gdf_pop["building"].isin(['residential', 'cabin'])].copy()
inf_pop = informal['pop'].sum()
#inf_area = formal['area'].mean()

#- student
stu = gdf_pop[gdf_pop["building"].isin(['student', 'dormitory'])].copy()
stu_pop = stu['pop'].sum()

bvpc_formal = round(formal['volume'].sum() / formal['pop'].sum()if formal['pop'].sum() != 0 else 0, 3)
bvpc_informal = round(informal['volume'].sum() / informal['pop'].sum() if informal['pop'].sum() != 0 else 0, 3)
bvpc_stu = round(stu['volume'].sum() / stu['pop'].sum() if stu['pop'].sum() != 0 else 0, 3)

print('FORMAL: Population: ', f_pop, ' with Building Volume Per Capita (BVPC):', bvpc_formal)
print('')
print('STUDENT RESIDENCE: Population: ', stu_pop, ' with Building Volume Per Capita (BVPC):', bvpc_stu)
print('')
print('INFORMAL: Population: ', inf_pop, ' with Building Volume Per Capita (BVPC)', bvpc_informal)

FORMAL: Population:  9768  with Building Volume Per Capita (BVPC): 83.156

STUDENT RESIDENCE: Population:  0  with Building Volume Per Capita (BVPC): 0

INFORMAL: Population:  968  with Building Volume Per Capita (BVPC) 42.559

Warning:

These are LoD1 3D City Models and works well in these types of areas.
LoD2 would offer a more representative BVpC (Ghosh, T; et al. 2020) value; when the complexity of the built environment increases.

Think about a house with living space in the roof structure, so called 'attic living', or an apartment / residential building with different levels, loft apartments and/or units in the turrets of a building.

*consider*: this area seperates building:cabin from building:residential to more precisely represent informal structures without typical roof trussess but account for social housing that does

To understand the performance in an Urban setting change cell [2] above:

jparams = osm3DUEstate ('University Estate') or osm3DsRiver ('Salt River') or osm3Dobs ('Observatory') (with residents per formal house = 4 and residents per informal structure = 3) and osm3DCPUT ('Cape Peninsula University of Technology (Bellville Campus))'

2. Interactive Visualization¶

You might want to create and share an html visualization.

In this example we identify building stock by color but you are limited only through your imagination and the data you have access too

In [23]:

#- pydeck needs geographic coords
gdf = gdf.to_crs(4326)

In [24]:

# -- get the location for pydeck
with warnings.catch_warnings():
    warnings.simplefilter('ignore')
    [xy] = gdf.dissolve().centroid

bbox = [gdf.total_bounds[0], gdf.total_bounds[1], 
        gdf.total_bounds[2], gdf.total_bounds[3]]

In [25]:

# have a look at the building type and amenities available
gdf['building'].unique()

Out[25]:

array(['yes', 'church', 'house', 'public', 'civic', 'office', 'retail',
       'clinic', 'school', 'cabin', 'garage', 'greenhouse', 'roof',
       'kindergarten', 'clubhouse', 'service', 'detached', 'shed'],
      dtype=object)

Building Stock: To differentiate a school, housing, retail, healthcare and community focused facilities (library, municipal office, community centre) we color the buildings - we harvest the osm tags [amenity and building type] directly.

In [26]:

# colour buildings based on use / amenity
def color(bld):
    #- formal house
    if bld == 'house' or bld == 'semidetached_house' or bld == 'terrace': #- add maisonette, duplex, etc. 
        return [255, 255, 204]                        #-grey
    if bld == 'apartments':
        return [252, 194, 3]                          #-orange 
    #- informal structure / social housing / student
    if bld == 'residential' or bld == 'dormitory' or bld == 'student' or bld == 'cabin':
        return [119, 3, 252]                          #-purple
        
    if bld == 'garage' or bld == 'parking':
        return [3, 132, 252]                          #-blue        
    if bld == 'retail' or bld == 'supermarket':
        return [253, 141, 60]
    if bld == 'office' or bld == 'commercial':
        return [185, 206, 37]
    if bld == 'school' or bld == 'kindergarten' or bld == 'university' or bld == 'college':
        return [128, 0, 38]
    if bld == 'clinic' or bld == 'doctors' or bld == 'hospital':
        return [89, 182, 178]
    if bld == 'community_centre' or bld == 'service' or bld == 'post_office' or bld == 'hall' \
    or bld ==  'townhall' or bld == 'police' or bld == 'library' or bld == 'fire_station' :
        return [181, 182, 89]
    if bld == 'warehouse' or bld == 'industrial':
        return [193, 255, 193]
    if bld == 'hotel':
        return [139, 117, 0]
    if bld == 'church' or bld == 'mosque' or bld == 'synagogue':
        return [225, 225, 51]
    else:
        return [255, 255, 204]

gdf["fill_color"] = gdf['building'].apply(lambda x: color(x))

In [27]:

#- look
gdf.head(2)

Out[27]:

	osm_id	building	building:levels	plus_code	ground_height	building_height	roof_height	address	amenity	residential	bottom_roof_height	building:use	geometry	fill_color
328118446	328118446.0	yes	1	CMXPX3XW+XQR	179.43	4.1	183.53	NaN	NaN	NaN	NaN	NaN	POLYGON ((18.47060 -33.50579, 18.47070 -33.505...	[255, 255, 204]
328118447	328118447.0	church	2	CPXWXXX3+XFX	179.04	6.9	185.94	Kerk Street Mamre 7347 Cape Town	place_of_worship	NaN	NaN	NaN	POLYGON ((18.47092 -33.50592, 18.47093 -33.505...	[225, 225, 51]

In [ ]:

## ~ (x, y) - bl, tl, tr, br  ~~ or ~~ sw, nw, ne, se
#area = [[[18.4377, -33.9307], [18.4377, -33.9283], [18.4418, -33.9283], [18.4418, -33.9307]]]
area = [[[bbox[0], bbox[1]], [bbox[0], bbox[3]], 
         [bbox[2], bbox[3]], [bbox[2], bbox[1]]]]

## ~ (y, x)
view_state = pdk.ViewState(latitude=xy.y, longitude=xy.x, zoom=16.5, max_zoom=19, pitch=72, 
                                   bearing=80)

land = pdk.Layer(
    "PolygonLayer",
    area,
    stroked=False,
    # processes the data as a flat longitude-latitude pair
    get_polygon="-",
    get_fill_color=[0, 0, 0, 1],
    #material = True,
    #shadowEnabled = True
)
building_layer = pdk.Layer(
    "PolygonLayer",
    gdf,
    #id="geojson",
    opacity=0.3,
    stroked=False,
    get_polygon="geometry.coordinates",
    filled=True,
    extruded=True,
    wireframe=False,
    get_elevation="building_height",
    #get_fill_color="[255, 255, 255]", #255, 255, 255
    get_fill_color="fill_color",
    get_line_color="fill_color",#[255, 255, 255],
    #material = True, 
    #shadowEnabled = True, 
    auto_highlight=True,
    pickable=True,
)

tooltip = {"html": "<b>Levels:</b> {building:levels} <br/> <b>Address:</b> {address}\
<br/> <b>Plus Code:</b> {plus_code} <br/> <b>Building Type:</b> {building}"}

#change the tooltip to show bus routes and comment out the previous
#tooltip = {"html": "<b>Route:</b> {name} <br/>"}

r = pdk.Deck(layers=[land,  building_layer],#, greenspaces_layer, p_layer, water_layer, r_layer], #
             #views=[{"@@type": "MapView", "controller": True}],
             initial_view_state=view_state,
             map_style = 'dark_no_labels', #pdk.map_styles.LIGHT,
             tooltip=tooltip)
#save
r.to_html("./result/interactiveOnly.html", offline=True)

on a laptop without a mouse:

trackpad left-click drag-left and -right;
Ctrl left-click drag-up, -down, -left and -right to rotate and so-on and
+ next to Backspace zoom-in and - next to + zoom-out.

Now you do your community. ~ If your area needs OpenStreetMap data and you want to contribute please follow the Guide.

3. Possible Secondary and Tertiary level conversations starters:

Topic	Secondary Level Questions	Tertiary Level Questions
Basic Understanding and Observations	- What types of buildings are most common in the area (houses, apartments, retail, etc.)? - Can you identify any patterns in the distribution of different types of buildings (e.g., are retail stores concentrated in certain areas)?	- How does the building stock composition (e.g., ratio of houses) correlate with the population? demographics (e.g., age distribution, household size) for the area will strengthen the analysis! - Analyze the relationship between building density and population. What urban planning theories can explain this relationship?
Spatial Relationships and Impacts	- How does the location of residential areas compare to the location of retail and commercial areas? - What impact might the density and distribution of buildings have on local traffic and transportation? - How might the population distribution affect the demand for local services such as schools, hospitals, and parks?	- Evaluate the accessibility of essential services (e.g., healthcare, education) in relation to the population and building types. - Assess the potential social and economic impacts of a proposed new residential or commercial development in the area.
Socioeconomic and Environmental Considerations	- Are there any correlations between the types of housing available and the household size? additional demographics (e.g., income level) for the area will strengthen the analysis! - How might the current building stock and population influence the local economy? demographics (e.g., age distribution, household size) for the area will strengthen the analysis! - What are some potential environmental impacts of the current building distribution, such as green space availability or pollution levels?	- How does the current building stock support or hinder sustainable development goals (e.g., energy efficiency, reduced carbon footprint)? - What strategies could be implemented to increase the resilience of the community to environmental or economic changes?
Future Planning and Development	- Based on the current building stock and population metrics, what areas might benefit from additional housing or commercial development? - How could urban planners use this information to improve the quality of life in the area? - What changes would you recommend to better balance residential, commercial, and recreational spaces?	- How might different zoning regulations impact the distribution of residential, commercial, and industrial buildings in the future? - Propose urban design solutions that could improve the sustainability and livability of the area, considering both current metrics and future projections.
Quantitative and Qualitative Research		- Design a research study to investigate the impact of building type diversity on community wellbeing. What methodologies would you use? - Analyze historical data to understand trends in building development and population growth. How have these trends shaped the current urban landscape? - Conduct a SWOT analysis (Strengths, Weaknesses, Opportunities, Threats) of the area based on the building stock and population metrics.

In [ ]: