<jupyter_start><jupyter_text>3D Building Models from Volunteered Public Data**We:**- **generate Level of Detail 1 (LoD1) 3D Building Models from [OpenStreetMap](https://en.wikipedia.org/wiki/OpenStreetMap) data,**- **export to [CityJSON](https://www.cityjson.org/) and** - **visualize with [pydeck](https://pydeck.gl/index.html).** This notebook differs from the previous in area-of-interest while the code accounts for "holes" (courtyards) in buildings.<jupyter_code>#load the magic

%matplotlib inline
import os
from pathlib import Path
import requests

import overpass
import osm2geojson

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

import ogr

import pydeck as pdk
path = Path('./data')<jupyter_output><empty_output><jupyter_text>Harvest OpenStreetMapWe query the [Overpass API](https://wiki.openstreetmap.org/wiki/Overpass_API) from within [Jupyter](https://jupyter.org/) and covert to [.geojson](https://geojson.org/)<jupyter_code>query = """
[out:json][timeout:30];
area[name='Cape Peninsula University of Technology (Bellville Campus)']->.a;

(
  (
    // I want all buildings
    way[building](area.a);

    // plus every building:part
    way["building:part"](area.a);
    // and relation type=multipolygon ~ to removed courtyards from buildings
    relation["building"]["type"="multipolygon"](area.a);
  );
-
  // excluding buildings with relation type=building role=outline
  // to remove outlines that surround building:part
  (
    // for every way in the input set select the relations of which it is an "outline" member
    rel(bw:"outline")["type"="building"];
    // back to the ways with role "outline"
    way(r:"outline");
  );
);
out body;
>;
out skel qt;
"""

url = "http://overpass-api.de/api/interpreter"
r = requests.get(url, params={'data': query})
#rr = r.read()
gj = osm2geojson.json2geojson(r.json())
# have a look at a random feature
#gj['features'][1]<jupyter_output><empty_output><jupyter_text>Calculate `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`. ~ original script at [hdb3d-code](https://github.com/ualsg/hdb3d-code/blob/master/hdb2d.py).<jupyter_code>storeyheight = 2.8

#-- iterate through the list of buildings and create GeoJSON features rich in attributes
footprints = {
    "type": "FeatureCollection",
    "features": []
    }

for i in gj['features']:
    f = {
    "type" : "Feature"
    }
    # at a minimum we only want building:levels tagged
    if 'building:levels' in i['properties']['tags']:
        f["properties"] = {}
        
        for p in i["properties"]:
        #print(i["properties"]['tags']['building:levels'])
        #break
        #-- store all OSM attributes and prefix them with osm_
            f["properties"]["osm_%s" % p] = i["properties"][p]
            osm_shape = shape(i["geometry"])
            #-- a few buildings are not polygons, rather linestrings. This converts them to polygons
            #-- rare, but if not done it breaks the code later
            if osm_shape.type == 'LineString':
                osm_shape = Polygon(osm_shape)
                
                #-- and multipolygons must be accounted for
            elif osm_shape.type == 'MultiPolygon':
                osm_shape = Polygon(osm_shape[0])
                #-- convert the shapely object to geojson
            f["geometry"] = mapping(osm_shape)

    #-- finally calculate the height and store it as an attribute (extrusion of geometry 
    ## -- will be done in the next script)
            f["properties"]['height'] = float(i["properties"]['tags']['building:levels']) * storeyheight + 1.3    
            footprints['features'].append(f)
    
#-- store the data as GeoJSON
with open(path/'cput_footprints.geojson', 'w') as outfile:
    json.dump(footprints, outfile)<jupyter_output><empty_output><jupyter_text>Reproject to a local coordinate systemThe [gdal](https://gdal.org/) [ogr2ogr](https://gdal.org/programs/ogr2ogr.html) program is a swiss-army knife for converting geospatial data and can be executed from jupyter directly.<jupyter_code>!ogr2ogr -f "GeoJSON" \
C:/{entire_path}/cput_footprints_r.geojson \
C:/{entire_path}/cput_footprints.geojson \
-s_srs EPSG:4326 -t_srs EPSG:32733<jupyter_output><empty_output><jupyter_text>ExtrudeWe modify the original [hdb3d.py](https://github.com/ualsg/hdb3d-code/blob/master/hdb3d.py) and execute here.<jupyter_code>!python cput_extruder.py<jupyter_output><empty_output><jupyter_text>New semantically rich CityJSON files should be created as per `line 29` of `cput_extruder.py` We can upgrade to CityJSON v1.0 and validate against the schema with...<jupyter_code>!cjio \
C:/{entire_path}/cput3dbuildings.json \
upgrade_version \
save C:/{entire_path}/cput3dbuildingsV1.json<jupyter_output><empty_output><jupyter_text>... and<jupyter_code>!cjio \
C:/{entire_path}/cput3dbuildingsV1.json \
validate<jupyter_output><empty_output><jupyter_text>CityJSON can be viewed with [azul](https://apps.apple.com/nl/app/azul/id1173239678?mt=12), the [CityJSON QGIS plugin](https://github.com/cityjson/cityjson-qgis-plugin), or the web-viewer [ninja](https://ninja.cityjson.org/).**We can have a look right here with [pydeck](https://pydeck.gl/index.html); the [deck.gl](https://deck.gl//) port for Python.** ~ which gives us the incredible ability to add tooltips.**pydeck needs an area and a center:** *we harvest from OpenStreetMap (this can be a suburb)*<jupyter_code>query = """
[out:json][timeout:25];
//gather results
(
  // query part for: “university”
  way['name'='Cape Peninsula University of Technology (Bellville Campus)'];
);
//print results
out body;
>;
out skel qt;
"""

url = "http://overpass-api.de/api/interpreter"
r = requests.get(url, params={'data': query})
#rr = r.read()
area = osm2geojson.json2geojson(r.json())
gdf = gpd.GeoDataFrame.from_features(area['features'])
bounds = gdf.geometry.bounds
x = gdf.centroid.x
y = gdf.centroid.y
bbox = [(bounds.minx, bounds.miny), (bounds.minx, bounds.maxy), 
        (bounds.maxx, bounds.maxy), (bounds.maxx, bounds.miny)]<jupyter_output><empty_output><jupyter_text>*I want additional visual effects to show the potential and power of 3D City Models; namely: parks, agricultural land and waterways (streams). We get this from OpenStreetMap as well.*<jupyter_code>query = """
[out:json][timeout:25];
area[name='Cape Peninsula University of Technology (Bellville Campus)'];
(
  // query
  way["leisure"~'pitch'](area);
);
// print results
out body;
>;
out skel qt;
"""

url = "http://overpass-api.de/api/interpreter"
p = requests.get(url, params={'data': query})
#rr = r.read()
green_spaces = osm2geojson.json2geojson(p.json())

query = """
[out:json][timeout:25];
area[name='Cape Peninsula University of Technology (Bellville Campus)'];
(
  // query
  way["natural"~'water'](area);
);
// print results
out body;
>;
out skel qt;
"""

url = "http://overpass-api.de/api/interpreter"
w = requests.get(url, params={'data': query})
#rr = r.read()
water_spaces = osm2geojson.json2geojson(w.json())

query = """
[out:json][timeout:25];
area[name='Cape Peninsula University of Technology (Bellville Campus)'];
(
  // query
  way['parking'](area);
);
// print results
out body;
>;
out skel qt;
"""

url = "http://overpass-api.de/api/interpreter"
p = requests.get(url, params={'data': query})
#rr = r.read()
p_spaces = osm2geojson.json2geojson(p.json())<jupyter_output><empty_output><jupyter_text>I experienced challenges with geojson and was forced to convert to pandas. ~ In order to make the most of the semantic data we need to extract the `osm_tags` from the dictionary: this can be any data you have.<jupyter_code>data = './data/cput_footprints.geojson'
json = pd.read_json(data)
build_df = pd.DataFrame()
# Parse the geometry out to Pandas
build_df["coordinates"] = json["features"].apply(lambda row: row["geometry"]["coordinates"])
build_df["height"] = round(json["features"].apply(lambda row: row["properties"]["height"]), 1)
#build_df["tags"] = json["features"].apply(lambda row: row["properties"]["osm_tags"])

#we want to display data so extract values from the dictionary 
build_df["tags"] = json["features"].apply(lambda row: row["properties"]["osm_tags"])
build_df['level'] = build_df['tags'].apply(lambda x: x.get('building:levels'))
build_df['name'] = build_df['tags'].apply(lambda x: x.get('name'))
## ~ (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][0][0], bbox[0][1][0]], [bbox[1][0][0], bbox[1][1][0]], 
         [bbox[2][0][0], bbox[2][1][0]], [bbox[3][0][0], bbox[3][1][0]]]]

## ~ (y, x)
view_state = pdk.ViewState(latitude=y[0], longitude=x[0], 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",
    build_df,
    #id="geojson",
    opacity=0.3,
    stroked=False,
    get_polygon="coordinates",
    filled=True,
    extruded=True,
    wireframe=False,
    get_elevation="height",
    get_fill_color="[255, 255, 255]", #255, 255, 255
    get_line_color=[255, 255, 255],
    #material = True, 
    #shadowEnabled = True, 
    auto_highlight=True,
    pickable=True,
)

greenspaces_layer =  pdk.Layer(
    "GeoJsonLayer",
    green_spaces,
    opacity=0.5,
    stroked=False,
    filled=True,
    wireframe=True,
    get_fill_color="[14, 140, 58]",
    get_line_color='[14, 140, 58]',
)

p_layer =  pdk.Layer(
    "GeoJsonLayer",
    p_spaces,
    opacity=0.2,
    stroked=False,
    filled=True,
    wireframe=True,
    get_fill_color="[128,128,128]",
    get_line_color='[128,128,128]',
)

water_layer =  pdk.Layer(
    "GeoJsonLayer",
    water_spaces,
    opacity=0.8,
    stroked=False,
    filled=True,
    wireframe=True,
    get_fill_color="[35, 35, 142]",
    get_line_color='[35, 35, 142]',
)

tooltip = {"html": "<b>Levels:</b> {level} <br/> <b>Name:</b> {name}"}

r = pdk.Deck(layers=[land, building_layer, greenspaces_layer, p_layer, water_layer], 
             #views=[{"@@type": "MapView", "controller": True}],
             initial_view_state=view_state,
             map_style = 'dark_no_labels', 
             tooltip=tooltip)

#save
r.to_html("./img/dark_cputBellville.html")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1558172397550636, 'fraction_numerical': 0.02298634008478568, 'mean_word_length': 4.153398058252427, 'pattern_counts': {'":': 5, '<': 46, '<?xml version=': 0, '>': 52, 'https://': 15, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '390798', 'license': 'MIT', 'n_tokens_mistral': 3875, 'n_tokens_neox': 3667, 'n_words': 871, 'path': 'osm_lod1_3dbuildingmodel_cput.ipynb', 'repo_name': 'AdrianKriger/osm_LoD1_3Dbuildings', 'size': '10615'}