pseudo-3D visualisation with Spatial Analysis

This notebook will:

1. produce an interactive pseudo-3D Building Model visualization - which a user can navigate, query, share that;

i) colour buildings by type (to easily visualize building stock)
ii) includes additional features (parks, bus rapid transit route, etc.)

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

i) use estimated average occupancy / household size to calculate a population estimate –with a previous census metric population growth rate and projected (future) population are also possible and
ii) quantify Building Volume per Capita.

3. further applications of Spatial Data Science

- calculate percentage homes and population with direct access to on-site renewable energy infrastructure –rooftop photovoltaic panels (PV) and solar water heaters (SWH).
- calculate the Annual Average Solar (photovoltaic) Potential, per home.

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

#- load the magic

%matplotlib inline
import time
from datetime import timedelta
import os
from pathlib import Path

import requests
import fiona
import overpass
import numpy as np
import json
import geojson
import pandas as pd

import shapely
from shapely.geometry import Polygon, shape, mapping
from shapely.strtree import STRtree

import city3D

import matplotlib.pyplot as plt
from matplotlib.patches import Polygon as MplPolygon
from matplotlib.collections import PatchCollection

#- works fine
Tstart = time.time()
import warnings
warnings.filterwarnings('ignore')

1. Interactive Visualization

Harvest OpenStreetMap - Query the Overpass API from within Jupyter and convert to .geojson.

Set an area-of-interest:

This is done: large area -> focus area or State (Province) -> Village (neighborhood / campus)

large = 'Western Cape'
focus = 'Mamre'                      # |'University Estate' |   Walmer Estate
osm_type = 'relation'

query = """
     [out:json][timeout:360];
    // --when areas have duplicate names given the world has a limited amount of uniquely named places
    area[name='{0}'] ->.b;
    // -- target area ~ can be way or relation
    wr(area.b)[name='{1}'];
    map_to_area -> .a;
        // I want all buildings
        (way['building'](area.a);
        // and relation type=multipolygon ~ to removed courtyards from buildings
        relation["building"]["type"="multipolygon"](area.a);
    );
    out count;
    out geom 2500;
    """.format(large, focus)

#- execute function from city3D, harvest buildings and return GeoDataFrameLite | home-baked gdf
gdf = city3D.overpass_to_gdf(query)
gdf.head(2)

	building	building:levels	addr:city	addr:postcode	addr:street	addr:suburb	amenity	denomination	diocese	heritage	...	street_vendor	brand	opening_date	self_service	capacity	construction	type	geometry	osm_id	osm_type
0	yes	1	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	POLYGON ((18.4706033 -33.5057876, 18.4706995 -...	328118446	way
1	church	2	Cape Town	7347	Kerk Street	Mamre	place_of_worship	moravian	Cape Town North	building	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	POLYGON ((18.4709153 -33.5059178, 18.4709287 -...	328118447	way

2 rows × 63 columns

NOTICE:

village will return a maximum of 2 500 buildings in any focus area.

#- some print statements
print('')
print(focus, 'has', len(gdf), 'buildings')
if int(len(gdf)) < 2500:
    print('\n\033[1mAll the buildings\033[0m  in', focus, 'have been harvested')
else:
    print('\n', int(len(gdf))-2500, "buildings have not been harvested.")

#-- try focus = 'Salt River' to see how a small urban suburb will perform or go over to the Suburb folder

Mamre has 2372 buildings

All the buildings  in Mamre have been harvested

Please do not burden the OpenStreetMap server with excessive calls for data.

If you need to investigate a larger area (> 2 500 buildings); choose suburb please.

Calculate building height:

We assume a building level is 2.8 meters high and add another 1.3 meters (to account for the roof) and create a new attribute height.

^{_{image adapted from the 3D geoinformation group at TUDelft}}

The Python code to execute the .calc_Bldheight function is in the city3D.py script

# -- execute function. calculate building height
gdf = city3D.bldHeights(gdf)

#- get the location for the pseudo-3D viz. combine all geometries
geom = shapely.unary_union(gdf['geometry'])
# centroid
xy = (geom.centroid.x, geom.centroid.y)

# bounding box
minx, miny, maxx, maxy = geom.bounds
bbox = [minx, miny, maxx, maxy]

~ In order to make the most of the semantic data we need to extract the osm_tags from the dictionary: and add it as tooltips to the visualization.

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

#- look
build_df = gdf.copy()
build_df.head(2)

	osm_id	address	building	building:levels	building:use	residential	amenity	operator	building_height	min_height	plus_code	footprint	geometry
0	328118446	None	yes	1	NaN	NaN	NaN	NaN	4.1	0.0	4FRWFFVC+P79	[[(18.471, -33.506), (18.471, -33.506), (18.47...	POLYGON ((18.4706033 -33.5057876, 18.4706995 -...
1	328118447	Moravian Church South Africa Kerk Street Mamre...	church	2	NaN	NaN	place_of_worship	NaN	6.9	0.0	4FRWFFVC+H9Q	[[(18.471, -33.506), (18.471, -33.506), (18.47...	POLYGON ((18.4709153 -33.5059178, 18.4709287 -...

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

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

len(build_df)

2344

Caveat

Besides the spatial aspect (the geometry) that defines the shape of a feature OpenStreetMap is a rich source of data. We have already introduced the standard and well-know building=* and building:levels=*. We now introduce another^*.

In many communities building’s, due to their age and the nature of change, are no longer used for their intended purpose. Consider warehouses along a waterfront that have transformed from storage to apartments or offices as an example.

To account for refurbishment –be as representative as possible– and conform to the OpenStreetMap Guide we typically tag these:
building=* ~ the original purpose + building:use=* ~ the current use.

^{*_{cell 26 and cell 44 introduce tags specific to buildings.
The bldHeights(gdf) function harvests a few more tags; such as: amenity=*, social_facility=* and the components of addr=*. Can you identify where the function does that? Hint: look in city3D.py}}

#- some data wrangling to account for when building:use is different from the original purpose 
#- (building=warehouse now loft apartments or =church now office, etc.)
df2 = build_df.copy()

# The entire operation is a single line using .loc
df2.loc[
    # The condition to find rows where 'building:use' is 'residential'
    # This check ensures the column exists, preventing a KeyError
    (df2['building:use'] == 'residential') & ~df2['building:use'].isna() 
    if 'building:use' in df2.columns else [False] * len(df2), 
    
    # The column to be updated
    'building'] = (
    # The value to assign to the 'building' column
    df2['building:use'] 
    if 'building:use' in df2.columns else None
)

#-- colour the building stock based on building:type

## we define specific colors
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 == '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 == 'civic' \
    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]

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

df2.head(2)

	osm_id	address	building	building:levels	building:use	residential	amenity	operator	building_height	min_height	plus_code	footprint	geometry	fill_color
0	328118446	None	yes	1	NaN	NaN	NaN	NaN	4.1	0.0	4FRWFFVC+P79	[[(18.471, -33.506), (18.471, -33.506), (18.47...	POLYGON ((18.4706033 -33.5057876, 18.4706995 -...	[255, 255, 204]
1	328118447	Moravian Church South Africa Kerk Street Mamre...	church	2	NaN	NaN	place_of_worship	NaN	6.9	0.0	4FRWFFVC+H9Q	[[(18.471, -33.506), (18.471, -33.506), (18.47...	POLYGON ((18.4709153 -33.5059178, 18.4709287 -...	[225, 225, 51]

Additional Features:

To show the potential and power of 3D City Models we can add additional features to the visualization; namely: bus rapid transit, parks, agricultural land and waterways (streams). We get this from OpenStreetMap as well..

query = """[out:json][timeout:360];
            // --main area
            area[name='{0}']->.b;
            // -- target area ~ can be way or relation
            wr(area.b)[name='{1}'];
            map_to_area -> .a;
                (
                // query
                //way["sport"](area.a);
                way[leisure~'track|pitch|park'](area.a);
                // relations (multipolygons)
                relation["leisure"~"track|pitch|park"](area.a);
                );
           // print results
           out geom;
           """.format(large, focus)

green_spaces = city3D.overpass_to_gdf(query)#, geojson=True)

query = """[out:json][timeout:360];
            // --main area
            area[name='{0}']->.b;
            // -- target area ~ can be way or relation
            wr(area.b)[name='{1}'];
            map_to_area -> .a;
                (
                // query
                way['waterway'='stream'](area.a);
                way['water'](area.a);
                );
            // print results
            out geom;
            """.format(large, focus)

water_spaces = city3D.overpass_to_gdf(query)#, geojson=True)

#query = """[out:json][timeout:180];
#        // --main area
#        area[name='{0}']->.b;
#        // -- target area ~ can be way or relation
#        wr(area.b)[name='{1}'];
#        map_to_area -> .a;
#            (
#            // query
#            way['landuse'='farmland'](area.a);
#            );
#        // print results
#        out geom;
#        """.format(large, focus)
#
#p_spaces = city3D.overpass_to_gdf(query, geojson=True)

# the bus route ~~ note we only choose routes with a 'colour' tag
query = """
[out:json][timeout:360];
area[name='{0}'];
// -- target area ~ can be way or relation
    // gather results
    (
    // query part for: “"bus route"”
    relation["type"="route"]["route"="bus"]['operator'="MyCiTi"]['colour'](area);
    );
// print results
out geom;
""".format(large)

Rgdf = city3D.overpass_to_gdf(query)

#green_spaces

# have a look at a random bus route
Rgdf.head(2)

	colour	from	name	network	operator	public_transport:version	ref	route	to	type	...	fee	interval	opening_hours	direction	note	website	description	geometry	osm_id	osm_type
0	#AACDD2	Civic Centre	Bus A01: Civic Centre – Airport [Suspended]	Cape Town IRT	MyCiTi	2	A01	bus	Airport	route	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	LINESTRING (18.4287357 -33.9197677, 18.4288394...	947075	relation
1	#AACDD2	Airport	Bus A01: Airport – Civic Centre [Suspended]	Cape Town IRT	MyCiTi	2	A01	bus	Civic Centre	route	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	LINESTRING (18.5958488 -33.9695727, 18.59574 -...	947076	relation

2 rows × 23 columns

# extract path and assign colour ~~ so the visualization matches the official documentation
Rgdf = Rgdf[Rgdf['colour'].notna()]

def hex_to_rgb(h):
    h = h.lstrip("#")
    #h = h.replace('#', '')
    return tuple(int(h[i : i + 2], 16) for i in (0, 2, 4))

Rgdf["colour"] = Rgdf["colour"].apply(hex_to_rgb)

Rgdf.head(2)

	colour	from	name	network	operator	public_transport:version	ref	route	to	type	...	fee	interval	opening_hours	direction	note	website	description	geometry	osm_id	osm_type
0	(170, 205, 210)	Civic Centre	Bus A01: Civic Centre – Airport [Suspended]	Cape Town IRT	MyCiTi	2	A01	bus	Airport	route	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	LINESTRING (18.4287357 -33.9197677, 18.4288394...	947075	relation
1	(170, 205, 210)	Airport	Bus A01: Airport – Civic Centre [Suspended]	Cape Town IRT	MyCiTi	2	A01	bus	Civic Centre	route	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	LINESTRING (18.5958488 -33.9695727, 18.59574 -...	947076	relation

2 rows × 23 columns

#- create pseudo-3D viz

file = './result/interactiveOnly.html' # will name and save html here

html = city3D.create_maplibre_3Dviz(
    result_dir = file,
    buildings_gdf = df2,
    water_gdf = water_spaces,
    green_gdf = green_spaces,
    brt_gdf = Rgdf,
    center = xy
)

#- uncomment to show
#city3D.show_interactive_html(html)

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.

---

GO FURTHER

2. Spatial Data Science (demography and housing)

Now that we have a visualization of building stock (buildings colorized by `use`); lets do some basic spatial analysis:

• We’ll calculate a population estimate, within our area of interest, and then

• calculate the Building Volume Per Capita (BVPC).

While calculating a population estimate 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

#-- lets have a look at the data we have
df2.head(2)

	osm_id	address	building	building:levels	building:use	residential	amenity	operator	building_height	min_height	plus_code	footprint	geometry	fill_color
0	328118446	None	yes	1	NaN	NaN	NaN	NaN	4.1	0.0	4FRWFFVC+P79	[[(18.471, -33.506), (18.471, -33.506), (18.47...	POLYGON ((18.4706033 -33.5057876, 18.4706995 -...	[255, 255, 204]
1	328118447	Moravian Church South Africa Kerk Street Mamre...	church	2	NaN	NaN	place_of_worship	NaN	6.9	0.0	4FRWFFVC+H9Q	[[(18.471, -33.506), (18.471, -33.506), (18.47...	POLYGON ((18.4709153 -33.5059178, 18.4709287 -...	[225, 225, 51]

#--we only want building=house or =apartment or =residential, etc.
#gdf = df2[df2["building"].isin(['house', 'semidetached_house', 'terrace', 'terraced', 'apartments', 'residential', 'dormitory', 'cabin', 'garage'])].copy()
#gdf.head(2)

2. a) Calculate a Population Estimate:

(with population growth rate and population projection possible too)

#- some data wrangling to replace 'bld:residential' to 'bld:student' if 'residential:student' -- like cell [15]:
gdf2 = df2.copy()
gdf2.loc[
    # The condition to find rows where 'residential' is 'student'
    # This check ensures the column exists, preventing a KeyError
    (gdf2['residential'] == 'student') & ~gdf2['residential'].isna() 
    if 'residential' in gdf2.columns else [False] * len(gdf2), 
    
    # The column to be updated
    'building'] = (
    # The value to assign to the 'building' column
    gdf2['residential'] 
    if 'residential' in gdf2.columns else None
)

#- some more data wrangling - Convert numeric columns
with pd.option_context("future.no_silent_downcasting", True):
    gdf2 = gdf2.assign(**{
        col: pd.to_numeric(
            gdf2[col].fillna(0).infer_objects(copy=False), errors='coerce'
        )
        for col in ['building:flats', 'building:units', 'beds', 'rooms', 'building:levels']
        if col in gdf2.columns
    })

print(len(gdf2))

2344

gdf2.head(2)

	osm_id	address	building	building:levels	building:use	residential	amenity	operator	building_height	min_height	plus_code	footprint	geometry	fill_color
0	328118446	None	yes	1	NaN	NaN	NaN	NaN	4.1	0.0	4FRWFFVC+P79	[[(18.471, -33.506), (18.471, -33.506), (18.47...	POLYGON ((18.4706033 -33.5057876, 18.4706995 -...	[255, 255, 204]
1	328118447	Moravian Church South Africa Kerk Street Mamre...	church	2	NaN	NaN	place_of_worship	NaN	6.9	0.0	4FRWFFVC+H9Q	[[(18.471, -33.506), (18.471, -33.506), (18.47...	POLYGON ((18.4709153 -33.5059178, 18.4709287 -...	[225, 225, 51]

gdf2['building'].value_counts()

building
house           1642
cabin            382
yes              157
garage           116
retail            11
civic              7
roof               7
school             5
church             3
public             2
construction       2
detached           2
office             1
clinic             1
greenhouse         1
kindergarten       1
clubhouse          1
guest_house        1
service            1
shed               1
Name: count, dtype: int64

This area is urban with single level housing units. To calculate a population estimate 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

#- average number of residents per formal house
f_house = 6
#- average number of residents per informal structure / social housing
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:apartments 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 (please follow the Guide) and local knowledge

c = gdf2.columns

def pop(row):
    #- formal house
    if row['building'] == 'house' or row['building'] == 'semidetached_house':
        return f_house
    if row['building'] == 'terrace' or row['building'] == 'terraced':
        if 'building:units' in c and row['building:units'] != 0:
            return row['building:units'] * f_house
        else:
            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 c and row['rooms'] != 0:
                return row['rooms']
            if 'building:flats' in c and row['building:flats'] != 0:
                return row['building:flats'] * inf_structure
        else:
            return inf_structure
    #-- social facility [shelter / carehome]
    if row['building'] == 'residential' and row['social_facility'] is not np.nan:
        if 'building:units' in c and row['building:units'] != 0:
            return row['building:units'] * inf_structure
        else: 
            return row['beds']
                
    #- formal apartment
    if row['building'] == 'apartments':
        return row['building:flats'] * 3
        
    #- private student residence 
    if row['building'] == 'student':
        if row['building:levels'] > 1:
            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 c and row['rooms'] != 0:
                return row['rooms']
            if 'beds' in c and row['beds'] != 0:
                return row['beds']
        else:
            return 3

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

est_pop = int(gdf2['pop'].sum())
print('The calculated population estimate is:', est_pop)

The estimated population is: 11380

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{11 120^{*}}{9048}}]}{12} * 100 = 1.47\%\]

* _{***Notice!*** The calculated population estimate (11176) is **NOT** the number in the formula (11 120). This community is frequently updated on OpenStreetMap and variations are common.}

Your Participation!

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

#- previous population
start_population = 9048

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

#-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.64 %

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

\[p = P_o * (1 + r)^{t} = p = 11120 * (1 + 0.0147)^{10} = 12870\]

Your Participation!

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

#- period in years from now
years = 10

#- account for non-residential areas without failure
#- helper function
def safe_population_estimate(est_pop, r, years):
    try:
        p = est_pop * (1 + (r / 100))**years
        return int(p)
    except Exception as e:
        print(f"Population estimate failed: {e}")
        return None  # keeps notebook running

#- execute function
p = safe_population_estimate(est_pop, r, years)

#- shows error and moves on
if p is not None:
    print(f"calculated population estimate {years} years from now: {p}")

estimated population 10 years from now: 13387

2. b) Building Volume Per Capita (BVPC)

BVPC: total building volume divided by population of a community

gdf2.crs

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

We need BVPC in a cubic meter but our data is in Decimal Degrees. We need to convert coordinates from a Geographic to a local Projected system.

To keep this extremely generic (use anywhere) we go with WGS84 / Universal Transverse Mercator (UTM).

#- internal geopandas function
gdf2.estimate_utm_crs()  

<Projected CRS: EPSG:32734>
Name: WGS 84 / UTM zone 34S
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: Between 18°E and 24°E, southern hemisphere between 80°S and equator, onshore and offshore. Angola. Botswana. Democratic Republic of the Congo (Zaire). Namibia. South Africa. Zambia.
- bounds: (18.0, -80.0, 24.0, 0.0)
Coordinate Operation:
- name: UTM zone 34S
- method: Transverse Mercator
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

Your Participation!

Fill in the proper epsg in the cell below

#- first line above: <Projected CRS: EPSG:xxxxx>
epsg = 32734

We first need to check the quality —topology— of the dataset

# prepare to plot (more buildings = more time) 
start = time.time()

gdf2 = gdf2.to_crs(epsg)
gdf2_copy = gdf2.copy()

geoms = gdf2_copy["geometry"].tolist()
tree = STRtree(geoms)

#- Query the tree for overlaps
# This returns two arrays: 'i' (index in geoms we are checking) and 'j' (index in the tree it overlapped with)
# Vectorized query: fastest way to find overlaps
i, j = tree.query(geoms, predicate="overlaps")

#- filter self-matches and get unique indices of all involved buildings
mask = i != j
overlap_idx = set(i[mask])

new_df1 = gdf2_copy.iloc[list(overlap_idx)].reset_index(drop=True)

end = time.time()
print('runtime:', str(timedelta(seconds=(end - start))))

runtime: 0:00:00.238659

Plot

Browse the saved './data/topologyFig' at your leisure

#- plot
def plot_geometries(df, ax=None, facecolor='none', edgecolor='purple', alpha=0.5):
    if ax is None:
        fig, ax = plt.subplots(figsize=(10,10))

    patches = []

    for geom in df['geometry']:
        if geom is None:
            continue

        if isinstance(geom, Polygon):
            # Exterior ring
            patches.append(MplPolygon(list(geom.exterior.coords), closed=True))
            # Interiors (holes)
            for interior in geom.interiors:
                patches.append(MplPolygon(list(interior.coords), closed=True))
        elif isinstance(geom, MultiPolygon):
            for poly in geom.geoms:
                patches.append(MplPolygon(list(poly.exterior.coords), closed=True))
                for interior in poly.interiors:
                    patches.append(MplPolygon(list(interior.coords), closed=True))

    pc = PatchCollection(patches, facecolor=facecolor, edgecolor=edgecolor, alpha=alpha)
    ax.add_collection(pc)
    ax.autoscale()
    ax.set_aspect('equal')
    return ax

# Example usage:
fig, ax = plt.subplots(figsize=(11, 11))
plot_geometries(gdf2_copy, ax=ax, facecolor='none', edgecolor='purple', alpha=0.2)
if len(new_df1) > 0:
    plot_geometries(new_df1, ax=ax, facecolor='none', edgecolor='red', alpha=0.5)
#-- save
plt.savefig('./data/topologyFig', dpi=300)
#plt.show()

	Challenges will be highlight in ‘Red’
	or none

Why are we doing this?

If you continue without fixing the challenges, the BVPC value will not be true.
The goal is to be as representative as possible.

If necessary; edit OpenStreetMap and fix the challenge please.

And remember to give the OpenStreetMap server at least a day before attempting the process again.

Alchemy is a process. Please be patient.

BVPC = total building volume divided by population of a community

#- area and volume
#gdf2['area'] = gdf2['geometry'].area
gdf2['area'] = gdf2['geometry'].apply(lambda geom: geom.area if geom else 0)
gdf2['volume'] = gdf2['area'] * gdf2['building_height']

#- remove the volume of the ground floor (unoccupied) when building:levels > 7 [this is an arbitrary number based on local knowledge]
#- typically this space is reserved for some other function: retail, etc. 
gdf2['volume'] = [
    (row['volume'] - row['area'] * 2.8) if (
        ('social_facility' not in gdf2.columns or pd.isna(row.get('social_facility')))
        and row['building:levels'] > 7
        and row['building'] in ['residential', 'apartments', 'student']
    ) else row['volume']
    for _, row in gdf2.iterrows()
]
 
gdf2['bvpc'] = np.where(
    gdf2['pop'] > 0,
    gdf2['volume'] / gdf2['pop'],
    np.nan
)
gdf2.tail(2)

	osm_id	address	building	building:levels	building:use	residential	amenity	operator	building_height	min_height	plus_code	footprint	geometry	fill_color	pop	area	volume	bvpc
2369	12289266	22 Clarkeson Street Mamre 7347 Cape Town	house	1	NaN	NaN	NaN	NaN	4.1	0.0	4FRWFFMF+W7J	[[(18.473, -33.515), (18.473, -33.515), (18.47...	POLYGON ((265302.8715881496 6288753.930217357,...	[255, 255, 204]	6.0	344.679759	1413.187012	235.531169
2370	12357148	2 Tol Street Mamre 7347 Cape Town	house	1	NaN	NaN	NaN	NaN	4.1	0.0	4FRWFFPJ+P7W	[[(18.481, -33.513), (18.481, -33.513), (18.48...	POLYGON ((265989.6485932828 6288995.734908563,...	[255, 255, 204]	6.0	327.808101	1344.013215	224.002202

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

count    2024.000000
mean       74.937556
std        49.834958
min         7.419412
25%        32.781743
50%        65.428197
75%       101.510961
max       399.892663
Name: bvpc, dtype: float64

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

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

Building Volume Per Capita (BVPC): 89.588

This BVPC value is for all the buildings; we only want buildings people live in (homes).

And we can seperate building:house from building:cabin and building:residential to understand 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 their home

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

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

#- student
stu = gdf2[gdf2["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:  9852.0  with Building Volume Per Capita (BVPC): 83.837

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

INFORMAL: Population:  1528.0  with Building Volume Per Capita (BVPC) 36.682

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: geo3D seperates building:cabin (shack) from building:residential to more precisely represent informal structures without typical roof trussess but account for social housing that does

Have a look at LoD2geo3D; to understand the performance of LoD2 models within the geo3D framework.

3. Further examples of Spatial Data Science (renewable energy):

Let’s attempt to understand the % of homes and population served with renewable energy.

SDG indicators are typically calculated at region and national scales.
Here, because we are working with highly detailed, local data, we can explore what a Tier 3 local indicator might look like at a neighbourhood level.

In this section 3. we evaluate SDG 7: Ensure access to affordable, reliable, sustainable and modern energy for all at a community level and calculate the proportion of residential units and population that have direct access to on-site renewable energy infrastructure –rooftop photovoltaic panels (PV) and solar water heaters (SWH).

a. Percentage of households served by rooftop renewable energy
b. Percentage of the population served by rooftop renewable energy
c. And then we go even further to calculate the Annual Solar Potential in MWh (theoretical maximum electricity) that homes can harvest from the sun over the course of one year.

#- harvest rooftop solar

query = """
     [out:json][timeout:360];
    // --when areas have duplicate names given the world has a limited amount of uniquely named places
    area[name='{0}'] ->.b;
    // -- target area ~ can be way or relation
    wr(area.b)[name='{1}'];
    map_to_area -> .a;
    (
        way["power"="generator"]["generator:source"="solar"](area.a);                // Catches simple generators
       //  way["power"="solar_photovoltaic_panel"](area.a);                          // Catches the alternate tag
    );
    out geom;
    """.format(large, focus)

#- execute function from city3D, harvest generator solar and return GeoDataFrameLite | home-baked gdf
sol = city3D.overpass_to_gdf(query)
sol = sol.to_crs(epsg)

if len(sol) < 0:
    print("\033[0m No rooftop solar are mapped in", focus)

#- look
sol.head(2)

	generator:method	generator:output:hot_water	generator:source	generator:type	location	power	start_date	area	generator:output:electricity	geometry	osm_id	osm_type
0	thermal	yes	solar	solar_thermal_collector	roof	generator	NaN	NaN	NaN	POLYGON ((266074.7976815852 6288576.002370844,...	1095737675	way
1	thermal	yes	solar	solar_thermal_collector	roof	generator	NaN	NaN	NaN	POLYGON ((266194.47713799647 6288764.539267579...	1095739760	way

#- we only want rooftop or roof
#sol = sol[sol['location'].isin(['rooftop', 'roof'])]
#sol = sol.reset_index(drop=True)
# how many of each type PV and SHW on all blds

#- the number of renewable in the AREA
sol['generator:method'].value_counts()

generator:method
thermal         57
photovoltaic     1
Name: count, dtype: int64

# join (link) rooftop renewable energy to the appropriate bld
def buildings_with_solar(gdf_buildings, gdf_solar):
    # Prepare output arrays
    solar_ids_per_building = [[] for _ in range(len(gdf_buildings["geometry"]))]
    solar_types_per_building = [[] for _ in range(len(gdf_buildings))]
    
    for i, b_geom in enumerate(gdf_buildings["geometry"]):
        for j, s_geom in enumerate(gdf_solar["geometry"]):
            #if b_geom.intersects(s_geom):
            if b_geom.contains(s_geom):
                solar_ids_per_building[i].append(gdf_solar["osm_id"].iloc[j])
                solar_types_per_building[i].append(gdf_solar["generator:method"].iloc[j])

    #- keep only unique values
    #unique_methods_per_building = [list(set(lst)) for lst in solar_types_per_building]
    
    gdf_buildings["solar_ids"] = solar_ids_per_building
    gdf_buildings["generator:method"] = solar_types_per_building #unique_methods_per_building
    gdf_buildings["has_solar"] = [len(lst) > 0 for lst in solar_ids_per_building]
    gdf_buildings["solar_ids"] = solar_ids_per_building
    
    return gdf_buildings

blds = buildings_with_solar(gdf2, sol)
blds.head(2)

	osm_id	address	building	building:levels	building:use	residential	amenity	operator	building_height	min_height	...	footprint	geometry	fill_color	pop	area	volume	bvpc	solar_ids	generator:method	has_solar
0	328118446	None	yes	1	NaN	NaN	NaN	NaN	4.1	0.0	...	[[(18.471, -33.506), (18.471, -33.506), (18.47...	POLYGON ((265041.6351569728 6289775.0592048075...	[255, 255, 204]	NaN	268.409550	1100.479153	NaN	[]	[]	False
1	328118447	Moravian Church South Africa Kerk Street Mamre...	church	2	NaN	NaN	place_of_worship	NaN	6.9	0.0	...	[[(18.471, -33.506), (18.471, -33.506), (18.47...	POLYGON ((265070.9766813367 6289761.325392997,...	[225, 225, 51]	NaN	371.168875	2561.065235	NaN	[]	[]	False

2 rows × 21 columns

#- we only want buildings people live in (homes)
blds = blds[blds["building"].isin(['house', 'semidetached_house', 'terrace', 'terraced', 'apartments', 'residential', 'dormitory', 'cabin', 'garage'])].copy()

3. a) Household rooftop solar

Percentage of households served by rooftop renewable energy

\[ \text{\% homes with renewable energy} = \frac{\text{Number of dwellings with mapped solar PV or SWH}}{\text{Total number of dwellings}} \times 100 \]

#solHms = round(len(sol) / len(gdf2) * 100, 2)

#- harvest columns
with_solar = sum(blds["has_solar"])
pop = est_pop #gdf["pop"]
total_homes = len(blds)

solHms = round((with_solar / total_homes) * 100, 2)

print('\033[1mPercentage homes\033[0m, in', focus,', with rooftop photovoltaic panels (PV) and solar water heaters (SWH):', solHms)

Percentage homes, in Mamre , with rooftop photovoltaic panels (PV) and solar water heaters (SWH): 2.62

NB: this number includes the OpenStreetMap building=garage building type. Go to Cell ±44 (above) to exclude this building type from the estimate.

3. b) Rooftop solar population

Percentage of population served by rooftop renewable energy

\[ \text{\% population with renewable energy} = \frac{\text{Number of residents with mapped solar PV or SWH}}{\text{Calculated population estimate}} \times 100 \]

#solPop = round(len(sol) / len(gdf2) * 100, 2)

pop_total = blds["pop"].sum()
pop_solar = blds["pop"][blds["has_solar"]].sum()

solPop = round((pop_solar / pop_total) * 100, 2)

print('\033[1mPercentage population\033[0m , in', focus,', with rooftop photovoltaic panels (PV) and solar water heaters (SWH):', solPop)

Percentage population , in Mamre , with rooftop photovoltaic panels (PV) and solar water heaters (SWH): 2.95

#- number of solar renewable per residential unit (home). 
#blds["generator:method"] = blds["generator:method"].apply(
#    lambda lst: np.nan if len(lst) == 0 else lst
#)
# number of solar renewable on HOMES
blds['generator:method'].explode().value_counts()

generator:method
thermal         55
photovoltaic     1
Name: count, dtype: int64

3. c) Solar potential (MWh)

In this section, we attempt to understand how much ‘fuel’ a rooftop can get from the sun.

We are calling the NASA POWER API (Prediction Of Worldwide Energy Resources). This is a global dataset that uses NASA satellite observations and weather models to tell us exactly how much solar radiation hits a specific coordinate on Earth.

What we are requesting:

Parameter: ALLSKY_SFC_SW_DWN (Global Horizontal Irradiance): a 30-year historical average of solar radiation. This ensures our communities solar potential is based on long-term climate trends rather than a single year of weather.

Source: NASA POWER Climatology API

Goal: To calculate the Annual Total GHI (kWh/m2/year). This value tells us the cumulative ‘solar pressure’ hitting our rooftops over an entire year, which we then use to calculate how many Megawatt-hours (MWh) of clean electricity our neighborhood can generate.

def get_ghi_data(lat, lon):#, year="2020"):
    #url = "https://power.larc.nasa.gov/api/temporal/daily/point"
    #url = "https://power.larc.nasa.gov/api/temporal/monthly/point"
    url = "https://power.larc.nasa.gov/api/temporal/climatology/point"
    params = {
        "parameters": "ALLSKY_SFC_SW_DWN", # This is NASA's GHI code
        "community": "RE",                 # Renewable Energy community
        "longitude": lon,
        "latitude": lat,
        "format": "JSON"
    }
    
    response = requests.get(url, params=params)
    data = response.json()
    
    # Extract the GHI values into a list: 13 values. one per month and the last one; ANN the average for the year.
    ghi_values = data['properties']['parameter']['ALLSKY_SFC_SW_DWN']
    long_term_monthly = ghi_values['ANN']
    
    # CONVERSION: Multiply by 365 to get the Yearly Total Sum
    annual_total_sum = long_term_monthly * 365
    
    return annual_total_sum

lat, lon = xy[1], xy[0] 
annual_avg = get_ghi_data(lat, lon)#, year)
print(f"Annual Average GHI: {round(annual_avg, 2)} kWh/m²/year")

Annual Average GHI: 2029.76 kWh/m²/year

Annual Solar Potential (MWh)

We use a simplified formula to provide a clear baseline.

\[ \text{Potential (MWh)} = \frac{(\text{Surface Area} \times \text{utilization factor}) \times \text{GHI}_{\text{annual}} \times 0.2}{1000} \]

_{Theoretical Framework: The annual energy output of a photovoltaic system (E) is determined by the product of the total solar resource (GHI), the active area of the array (A), and the system’s overall efficiency (η), adjusted by a Performance Ratio (PR) to account for real-world losses. — based on NREL (2022) & IEC 61724-1. We then adapt this formula and represents a combined value of 25% nominal panel efficiency and a 0.80 Performance Ratio with a single 0.20 system efficiency value and account for usable area, a heuristic for gabled roofs.}

Your Participation!

Fill in the utilization_factor below

As a ‘rule-of-thumb’ a community with traditional gabled houses: utilization_factor = 0.4 (less than half), while a high-density suburb with flat-roofed apartments: utilization_factor = 0.6

# on average, how much of the roof faces the sun? adjust based on roof types
utilization_factor = 0.4  

blds['solar_mwh'] = (((blds['area'] * utilization_factor) * (annual_avg) * 0.20) / 1000)
blds.head(2)

	osm_id	address	building	building:levels	building:use	residential	amenity	operator	building_height	min_height	...	geometry	fill_color	pop	area	volume	bvpc	solar_ids	generator:method	has_solar	solar_mwh
7	656840974	39 Dove Lane Mamre 7347 Cape Town	house	1	NaN	NaN	NaN	NaN	4.1	0.0	...	POLYGON ((264864.98193212313 6288741.008398377...	[255, 255, 204]	6.0	38.658416	158.499506	26.416584	[]	[]	False	6.277400
8	656840975	37 Dove Lane Mamre 7347 Cape Town	house	1	NaN	NaN	NaN	NaN	4.1	0.0	...	POLYGON ((264865.3504529182 6288749.140981195,...	[255, 255, 204]	6.0	39.833421	163.317028	27.219505	[1096772168]	[thermal]	True	6.468199

2 rows × 22 columns

# Calculate the average annual solar potential
average_solar_potential = blds['solar_mwh'].mean()

print(" \033[1mThe average solar potential, per home\033[0m , for", focus, "is:", round(average_solar_potential,2), "MWh/year")

 The average solar potential, per home , for Mamre is: 17.01 MWh/year

NB: this number includes the OpenStreetMap building=garage building type. Go to Cell ±42 to exclude this building type from the estimate.

What does this MWh/year value mean?

To put the value in context, 15 MWh/year:

• is enough to provide 100% of the electricity for 4 to 5 average UK homes (which use ~3.4 MWh each) or 1.5 average US homes (~10.7 MWh each)

• is enough power to drive an Electric Vehicle for 75,000 kilometers –that’s almost two full trips around the Earth.

• saves roughly 10 metric tons of Carbon Dioxide from entering the atmosphere.

Sanity Check!

• In Cell ±44 we excluded non-residential building types =office, commercial, retail, warehouse, industrial, etc. from the analysis.

• Cape Town typically yields ~1.6–1.7 MWh per installed kWp per year; the higher per-household values reported here reflect rooftop potential derived from available area –that considers a utilization_factor, not a 1 kWp system.

  A 1 kWp solar PV system requires approximately 5–8 m² of panel area (e.g. panels of roughly 1 m × 1.7 m, depending on technology).

  We are NOT asking: How much energy (MWh/year) would a single 1 kWp PV system generate on a roof?
  We are asking: How much energy (MWh/year) could these roofs harvest, given their available area and a realistic utilization factor?

• The BVPC Warning applies here too. These are LoD1 3D City Models, which represent buildings as simple extrusions.

  LoD2 models —that capture roof form (e.g. gable, hipped, mansard, domes)— would provide more representative estimates of both BVPC and Average Annual Solar Potential. In such cases, the utilization_factor becomes less critical, as usable roof geometry is explicitly modelled.

Have a look at LoD2geo3D; to understand the performance of LoD2 models within the geo3D framework.

4. Possible Secondary and Tertiary level conversations starters:

communicate and exchange ideas and understanding

Topic	Secondary Level Questions	Tertiary Level Questions
Geography	- Talk about the main difference between a globe and a map, and why we use map projections to represent the Earth on a flat surface - Explain why different map projections are used for different purposes. For example, why might a Mercator projection be useful for navigation, but not for comparing the sizes of countries?	- Discuss why it is necessary to convert geographic coordinates (latitude and longitude) to a projected coordinate system in the context of the geospatial sciences. What are some potential issues if this conversion is not done? - How does geodesy contribute to the geospatial sciences?
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.

---

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

focus = ’University Estate’ or ’Salt River’ or ’Observatory’ (with residents per formal house = 4 | 5 for Salt River and residents per informal structure = 3) and ’Cape Peninsula University of Technology (Bellville Campus)’
osm_type = ’relation’ with CPUT (Bellville Campus) as ’way’

Salt River and Observatory 2013 population was 6 577 and 9 207 (as per City of Cape Town Open Data). University Estate is a tiny urban neighbourhood and population data is typically not available at that scale.

Tend = time.time()
print('runtime:', str(timedelta(seconds=(Tend - Tstart))))