osm_LoD1_3DCityModel (No Internet)

This is a stand-alone Notebook for use in low resource settings. NO internet connection is necessary.

While an internet connection is NOT necessary you will NEED to have sourced an osm.pbf.

The purpose of this notebook is to walk a user through osm_LoD1_3DCityModel.

This notebook will:

1. allow the user to create a Level-of-Detail 1 (LoD1) 3D City Model.

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

The suburb processing option is meant for areas with more than for 2 500 buildings.

#- load the magic
import time
from datetime import timedelta
import tempfile

import os
from itertools import chain
import math
import requests
import overpass
import copy
import json
import numpy as np
import pandas as pd
import topojson as tp

import shapely
from shapely.geometry import Point, Polygon, MultiPolygon, polygon 
from shapely.ops import snap, transform
from shapely.strtree import STRtree

import city3D

import pyproj
from osgeo import gdal, ogr, osr

import triangle as tr

from openlocationcode import openlocationcode as olc

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

Tstart = time.time()

import warnings
warnings.filterwarnings('ignore')

A parameter.json defines the path and files.

jparams = json.load(open('wStock_param.json'))       
#jparams = json.load(open('sRiver_param25m.json'))       
#jparams = json.load(open('saao_param.json'))  
#jparams = json.load(open('mamre_param.json'))  

area of interest	elevation model	CityJSON and metadata

Harvest OpenStreetMap - interogate an osm.pbf (“Protocolbuffer Binary Format”) from within Jupyter and convert to .geojson.

PLEASE SUPPLY YOUR OWN osm.pbf.

Either crop an area directly from OpenStreetMap with the official tool, select a predefined area from any number of providers, such as Geofabrik, or…

… download your own. Provincial extracts for South Africa are available here: http://download.openstreetmap.fr/extracts/africa/south_africa/

# Input OSM PBF file
input_pbf = "./data/CapeTown.osm.pbf"

#input_pbf = "./data/south-africa-latest.osm.pbf"

Lets first harvest the boundary of the area; we want to interogate

start = time.time()

#- execute function from city3D and return GeoDataFrameLite | home-baked gdf
aoi = city3D.extract_boundaries_by_name(input_pbf, jparams)

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

#- suppose 'aoi' is your GeoDataFrameLite or list of geometries
aoi.head(2)

runtime: 0:00:02.450105

	boundary	geometry	name	osm_id	other_tags	place	type
0	place	MULTIPOLYGON (((18.4407807 -33.9286635, 18.440...	Woodstock	2034285	"wikidata"=>"Q3644460"	suburb	boundary

# gt the bounding box (BBOX) of the boundary

geoms = aoi['geometry'].tolist()
#- combine all geometries into a single union
combined_geom = shapely.unary_union(geoms)  # returns Polygon or MultiPolygon
#- compute bounding box
minx, miny, maxx, maxy = combined_geom.bounds
#extent = [minx - 250, miny - 250,maxx + 250, maxy + 250]

Only harvest what we need from the osm.pbf.

start = time.time()

gdal.UseExceptions()
gdal.SetConfigOption("OGR_GEOMETRY_ACCEPT_UNCLOSED_RING", "NO") 
#gdal.SetConfigOption("USE_CUSTOM_INDEXING", "NO")

# GDAL Virtual File System (VSI) to avoid writing to disk
geojson_vsimem = "/vsimem/temp.geojson"

#- GDAL VectorTranslate to extract only buildings & fix geometries
gdal.VectorTranslate(
    geojson_vsimem,                                           # Output as in-memory GeoJSON
    input_pbf,                                                # Source OSM PBF file
    format="GeoJSON",                                         # Output format
    layers=["multipolygons"],                                 # Extract only multipolygons
    options=["-where", "building IS NOT NULL", "-makevalid", 
                          "-spat", str(minx), str(miny), str(maxx), str(maxy)]  # Filter buildings & fix geometries
)

#- execute and return GeoDataFrameLite | home-baked gdf
gdf = city3D.read_vsimem_geojson(geojson_vsimem)

#- cleanup VSI Memory
gdal.Unlink(geojson_vsimem)

# show gdf
#gdf.head()

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

ERROR 1: Non closed ring detected.
ERROR 1: Non closed ring detected.

runtime: 0:00:01.764839

gdf.head(2)
#len(gdf)

	amenity	building	craft	geometry	historic	leisure	man_made	name	office	osm_id	osm_way_id	other_tags	shop	sport	tourism	type
0	None	school	None	MULTIPOLYGON (((18.4390374 -33.9341504, 18.439...	None	None	None	None	None	11029147	None	"building:levels"=>"2","ref:ZA:emis"=>"1033103...	None	None	None	multipolygon
1	marketplace	retail	None	MULTIPOLYGON (((18.4580079 -33.9273091, 18.458...	None	None	None	The Neighbourgoods Market	None	12227309	None	"addr:city"=>"Cape Town","addr:suburb"=>"Woods...	None	None	None	multipolygon

# Convert valid strings, ignore None/NaN
def safe_convert(tag_string):
    if isinstance(tag_string, str):
        try:
            # Replace "=>" with ":" and fix newlines
            formatted_string = "{" + tag_string.replace("=>", ":").replace("\n", " ") + "}"
            return json.loads(formatted_string)  # Parse safely
        except json.JSONDecodeError:
            return {}  # Return empty dict on failure
    return {}  # Return empty dict if NaN or None

# Apply conversion function
gdf["tags"] = gdf["other_tags"].apply(safe_convert)

# Extract values safely - Normalize the 'tags' column to create a new DataFrame
tags_df = pd.json_normalize(gdf['tags'])
# Join the new columns back to the original GeoDataFrame
gdf = pd.concat([gdf, tags_df], axis=1)
# (Optional) Drop the original 'tags' column
gdf = gdf.drop(columns=['other_tags'])

# Ensure a single 'osm_id' column
if 'osm_id' in gdf.columns:
    if 'osm_way_id' in gdf.columns:
        gdf['osm_id'] = [o if pd.notna(o) else w 
                         for o, w in zip(gdf['osm_id'], gdf['osm_way_id'])]
        gdf = gdf.drop(columns=['osm_way_id'])
elif 'osm_way_id' in gdf.columns:
    gdf = gdf.rename(columns={'osm_way_id': 'osm_id'})

#gdf = gdf[gdf.geometry.apply(lambda x: x.within(aoi.unary_union))]
gdf = gdf[gdf.geometry.apply(lambda x: x.within(shapely.unary_union(aoi.geometry)))]
gdf.crs = "EPSG:4326"

gdf.head(2)

	amenity	building	craft	geometry	historic	leisure	man_made	name	office	osm_id	...	studio	guest_house	abandoned:building	motorcycle:rental	second_hand	drink:coffee	bus	network	opening_date	unisex
1	marketplace	retail	None	MULTIPOLYGON (((18.4580079 -33.9273091, 18.458...	None	None	None	The Neighbourgoods Market	None	12227309	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
2	None	apartments	None	MULTIPOLYGON (((18.450818 -33.9279241, 18.4508...	None	None	None	Church Square	None	12249345	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

2 rows × 149 columns

ts = gdf[gdf['building'].notna()]
#len(ts)
print('\n', len(ts), "buildings have been harvested from", input_pbf)

 3797 buildings have been harvested from ./data/CapeTown.osm.pbf

ts.head(2)

	amenity	building	craft	geometry	historic	leisure	man_made	name	office	osm_id	...	studio	guest_house	abandoned:building	motorcycle:rental	second_hand	drink:coffee	bus	network	opening_date	unisex
1	marketplace	retail	None	MULTIPOLYGON (((18.4580079 -33.9273091, 18.458...	None	None	None	The Neighbourgoods Market	None	12227309	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
2	None	apartments	None	MULTIPOLYGON (((18.450818 -33.9279241, 18.4508...	None	None	None	Church Square	None	12249345	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

2 rows × 149 columns

# basic cleaning to harvest building=* (no building:part=*) and building=levels tags only

#- we only want buildings with =levels data
ts['building:levels'] = pd.to_numeric(ts['building:levels'], errors='coerce')
ts = ts[ts['building:levels'] > 0]

#- without building:part
ts = ts[ts.get("building:part").isnull()] if "building:part" in ts else ts
print('\n\033[1m', jparams['FocusArea'], 'has \033[0m', len(ts), 'buildings')

 Woodstock has  3715 buildings

# have a look
ts.tail(2)

	amenity	building	craft	geometry	historic	leisure	man_made	name	office	osm_id	...	studio	guest_house	abandoned:building	motorcycle:rental	second_hand	drink:coffee	bus	network	opening_date	unisex
4388	None	semidetached_house	None	MULTIPOLYGON (((18.4497301 -33.9341967, 18.449...	None	None	None	None	None	1340498943	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
4389	None	warehouse	None	MULTIPOLYGON (((18.4479047 -33.9257214, 18.447...	None	None	None	None	None	1423081172	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

2 rows × 149 columns

#- coordinate reference system
ts.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 the Projected Coordinate Reference System.

#- estimate utm: internal geopandas function
ts.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 espg in the cell below

#- fill <Projected CRS: EPSG:32734> from above here epsg = EPSG:32734
epsg = 'EPSG:32734'

#project blds
ts = ts.to_crs(epsg)
#project aoi
aoi = aoi.to_crs(epsg)

1. Create LoD1 3D City Model

Now we process.

aoibuffer = aoi.copy()

def buffer01(row):
    with np.errstate(invalid='ignore'):
        return row.geometry.buffer(150, cap_style=3, join_style=2)

aoibuffer['geometry'] = aoibuffer.apply(buffer01, axis=1)

#- suppose 'aoi' is your GeoDataFrameLite or list of geometries
geoms = aoibuffer['geometry'].tolist()
#- combine all geometries into a single union
combined_geom = shapely.unary_union(geoms)  # returns Polygon or MultiPolygon
#- compute bounding box
minx, miny, maxx, maxy = combined_geom.bounds

extent = [minx - 250, miny - 250,
          maxx + 250, maxy + 250]

Now the DEM
one is available at raster

gdal.SetConfigOption("GTIFF_SRS_SOURCE", "GEOKEYS")
gdal.UseExceptions() 

# set the path and nodata
OutTile = gdal.Warp(jparams['projClip_raster'], 
                    jparams['in_raster'],
                    dstSRS=epsg,
                    srcNodata = jparams['nodata'],
                    #-  dstNodata = 0,
                    #-- outputBounds=[minX, minY, maxX, maxY]
                    outputBounds = [extent[0], extent[1], extent[2], extent[3]])
OutTile = None 

#- convert raster to XYZ in-memory

#- virtual in-memory path
xyz_mem_path = "/vsimem/temp_xyz.xyz"  
gdal.Translate(xyz_mem_path, jparams['projClip_raster'], format="XYZ")  

#- read XYZ from GDAL's in-memory file
xyz_vsimem = gdal.VSIFOpenL(xyz_mem_path, "rb")
xyz_bytes = gdal.VSIFReadL(1, gdal.VSIStatL(xyz_mem_path).size, xyz_vsimem)
gdal.VSIFCloseL(xyz_vsimem)
#- cleanup in-memory file
gdal.Unlink(xyz_mem_path) 

0

Buildings

#- simplify geometry
ts = city3D.GeoDataFrameLite(ts)
geojson_dict = json.loads(ts.to_json())

for feat in geojson_dict["features"]:
    if feat.get("type") is None:
        feat["type"] = "multipolygon"
    if feat.get("geometry") is None:
        feat["geometry"] = {"type":"MultiPolygon","coordinates":[]}

topo = tp.Topology(geojson_dict, prequantize=False, winding_order='CCW_CW')
simplified_geojson = topo.toposimplify(0.25).to_geojson()

#- back into home-baked gdf
ts = city3D.GeoDataFrameLite.from_json(simplified_geojson)
ts.crs = epsg

#- highlight crossing features (buildings). more buildings = more time
start = time.time()

ts = ts.to_crs(epsg)
ts_copy = ts.copy()

geoms = ts_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 = ts_copy.iloc[list(overlap_idx)].reset_index(drop=True)

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

runtime: 0:00:00.513966

Plot

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

#%matplotlib

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(ts_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()

	Typical challenges will be highlight in Red.
	or none
If you continue without fixing the challenges, the LoD1 City Model will not be error free.	Our goal is a topologically correct surface.

Why are we doing this?

If you continue without fixing the challenges, the LoD1 3D City Model will NOT conform to the ISO 19107 spatial schema for 3D primatives.
It is possible to create a very high quality product from minimal resources.
Please create a high quality product.

If necessary; edit OpenStreetMap and fix the challenge please.

And remember.
Many Planet.osm mirrors release a fresh .osm.pbf EVERYDAY!
Give the OpenStreetMap server at least a day before attempting the process again.

Alchemy is a process. Please be patient.

# set the path to the projected, cliped elevation
src_filename = jparams['projClip_raster']

src_ds = gdal.Open(src_filename) 
gt_forward = src_ds.GetGeoTransform()
rb = src_ds.GetRasterBand(1)

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

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

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

# -- execute function. write geoJSON
dis = city3D.bldHeights(ts)

start = time.time()

dis_c = dis.copy()
dis_c.drop(dis.index[dis['building'] == 'bridge'], inplace = True)
dis_c.drop(dis.index[dis['building'] == 'roof'], inplace = True)

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

runtime: 0:00:00.007202

dis.head(2)

	osm_id	address	building	building:levels	building:use	building:flats	building:units	beds	rooms	residential	amenity	social_facility	operator	building_height	min_height	plus_code	footprint	geometry
0	12227309	The Neighbourgoods Market Woodstock Cape Town	retail	1.0	NaN	NaN	NaN	NaN	NaN	NaN	marketplace	NaN	NaN	4.1	0.0	4FRW3FF5+27C	[[(265023.481, 6242993.344), (265021.289, 6242...	POLYGON ((265023.481386 6242993.344091, 265021...
1	12249345	Church Square 34 Dickson Street 7915 Woodstoc...	apartments	6.0	NaN	105	NaN	NaN	NaN	NaN	None	NaN	NaN	18.1	0.0	4FRW3FC2+R5G	[[(264360.391, 6242908.633), (264361.062, 6242...	POLYGON ((264360.390586 6242908.633245, 264361...

prepare the elevation for the TIN

#- 
#dis_c = dis.copy()

#- prepare xyz (more buildings = more time)
start = time.time()

# Convert bytes to DataFrame
xyz_str = xyz_bytes.decode("utf-8")  # Decode to string

dtype_spec = {
    "x": np.float32,  # Reduce precision from float64 to float32 (saves memory)
    "y": np.float32,
    "z": np.float32
}

#df = pd.read_csv(jparams['xyz'], delimiter = ' ', header=None, names=["x", "y", "z"])                                                                                
df = pd.read_csv(pd.io.common.StringIO(xyz_str), delimiter=" ",  header=None, 
                 names=["x", "y", "z"], dtype=dtype_spec)                     # in memory fastest

#- Create the shapely 'geometry' column directly (Vectorized) and GeoDataFrameLite | home-baked gdf
df['geometry'] = df.apply(lambda row: Point(row['x'], row['y']), axis=1) 
gdf = city3D.GeoDataFrameLite(df) 
gdf.crs = epsg

# --- cleanup ---
gdf = gdf[gdf['z'] != jparams['nodata']] 
gdf.reset_index(drop=True, inplace=True)
gdf = gdf.round(2)
#print(len(gdf))

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

runtime: 0:00:00.152461

#dis.tail(2)

Prepare for Triangle:

The Python code to execute the city3D.functions are in the city3D.py script

#- harvest the building vertices, combine with the elevation, create regions and segments for Triangle
coords, regions, segments = city3D.prepareTri(gdf, dis_c, aoibuffer)

Triangle

A = dict(vertices=np.array(coords), segments=np.array(segments), #holes=np.array(holes), 
         regions=np.array(regions))

# 'p'  = Triangulate the PSLG: Delauney triangulation with segments (building outlines) as constraints.
# 'Y' = Do NOT add Steiner points
# 'A'  = Attribute triangles with region IDs
# 'z'  = Zero-based indexing (prevents index errors)
Tr = tr.triangulate(A, 'pYAz')

#- the vertices 
final_verts_2d = Tr['vertices']

#- 
z_cache = {(row.x, row.y): row.z for row in gdf.itertuples()}

## -- we triangulate in 2D and project into 3D space. the vertices of the building outlines need a 'z'-value
final_verts_3d = []
for x, y in final_verts_2d:
    x_r, y_r = x, y
    
    # 2. Check if we already have the Z value in our GDF points
    if (x_r, y_r) in z_cache:
        z = z_cache[(x_r, y_r)]
    else:
        # 3. Only query the raster if the point is a new vector/Steiner vertex
        z = float(city3D.rasterQuery2(x, y, gt_forward, rb))

    final_verts_3d.append([x, y, z])

final_verts_3d = np.array(final_verts_3d)

#s- eparate triangles by their Region ID for CityJSON
tris = Tr['triangles']
tri_attr = Tr['triangle_attributes'].flatten()

CityJSON

#- 
minz = gdf['z'].min()
maxz = gdf['z'].max()

create CityJSON

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

# -- execute function. create CityJSON
crs = epsg[5:]

city3D.output_cityjson(extent, minz, maxz, tris, tri_attr, final_verts_3d, dis, jparams, gt_forward, rb, crs)

src_ds = None

Go over to Ninja the online CityJSON viewer and explore!

You are welcome to further investigate the quality of a 3D Model.

The val3dity web app will test CityJSON geometric primitives.

If you parse the result of this notebook through val3dity it will return a report with an invalid TINRelief and error.

This particular area contains Buildings with courtyards. The courtyards (polygons) are islands of terrain disconnected from the larger TINRelief (shell); thus the error *.

^{* _{Don’t take my word for it. Test and see for yourself! saao_param.json (South African Royal Observatory, Cape Town) will produce a 100% topologically correct Open Geospatial Consortium (OGC) standard LoD1 3D model that conforms to the ISO 19107 spatial schema for 3D primatives [connecting and planar surfaces, correct orientation of the surfaces and watertight volumes]}}

NOW

To understand the value and usefulness of a 3D City Model; parse the result of this Notebook through CityJSONspatialDataScience.ipynb to workthrough an example of:

• calculate a population estimate,

• quantify Building Volume per Capita, and

• calculate the Annual Average Solar (photovoltaic) Potential, per home.

2. Possible Secondary and Tertiary level conversations starters:

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?

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

runtime: 0:00:22.644852