Exposure in Python
For this tutorial you need to install matplotlib, folium, mapclassify and geopandas. Also the software package databaselibfrom GFZ is needed.
# install necessary libraries - remove # to install
#!pip config set global.extra-index-url https://git.gfz-potsdam.de/api/v4/projects/2940/packages/pypi/simple
#!pip install databaselib matplotlib folium mapclassify geopandas
Connect to a database¶
import geopandas
from shapely import from_wkt
from databaselib import SpatiaLiteDatabase
db = SpatiaLiteDatabase('data/2023-011_Schorlemmer-et-al_AND.Andorra.db')
db.connect(init_spatial_metadata=False)
Select data¶
If we want to know any information about an entity (a tile or a building), for example the information at a certain location (or Quadkey tile), we first need to query the ID of the feature. A simple example is given here:
quadkey = '120222210222313200'
db.cursor.execute(f"SELECT id, ST_AsText(geom) FROM entity WHERE quadkey = '{quadkey}'")
entity_id,geom = db.cursor.fetchone()
print(entity_id, geom)
521 MULTIPOLYGON(((1.484528 42.572298, 1.484528 42.57331, 1.483154 42.57331, 1.483154 42.572298, 1.484528 42.572298)))
Info on an entity¶
The entity is just a geographical description of a feature, but does not contain any semantic information. This is stored in the asset table. As we are making a probabilistic assessment, multiple assets could exist for one entity.
If we would like to know the total population and structural value in the Quadkey that has been defined previously, we can simply query all assets with this entity ID.
db.cursor.execute(f"SELECT SUM(night), SUM(structural) FROM asset WHERE entity_id = {entity_id}")
population, structural_value = db.cursor.fetchone()
print(f"Expected population: {population}")
print(f"Expected total structural value: ā¬{structural_value:,.2f}")
Expected population: 37.6324613306 Expected total structural value: ā¬1,332,212.74
Inner join the assets and the entities. Instead of querying by entity ID, we can also query by Quadkey directly, by doing an inner join on the entity table, giving the same result.
db.cursor.execute(f"""
SELECT SUM(night), SUM(structural)
FROM entity
INNER JOIN asset
ON entity.id = asset.entity_id
WHERE quadkey = '{quadkey}'
""")
population, structural_value = db.cursor.fetchone()
print(f"Expected population: {population}")
print(f"Expected total structural value: ā¬{structural_value:,.2f}")
Expected population: 37.6324613306 Expected total structural value: ā¬1,332,212.74
Asset types in a tile¶
On top of population and monetary values, we can query all structural types inside a tile, by joining the asset table with the taxonomy table, that stores all building using the GEM Taxonomy. This is done as follows:
db.cursor.execute(f"""
SELECT taxonomy_string, number
FROM asset
INNER JOIN taxonomy ON asset.taxonomy_id = taxonomy.id
WHERE entity_id = {entity_id}
ORDER BY number DESC
""")
total_number_of_buildings = 0
print("Probabilistically, in Quadkey `{quadkey}`, there are the following number of buildings expected for each of the structural types:")
for taxonomy_string, number in db.cursor:
print(f"{number:.4f} for type `{taxonomy_string}`.")
total_number_of_buildings += number
print(f"In total, there are {total_number_of_buildings} buildings expected in the tile.")
Probabilistically, in Quadkey `{quadkey}`, there are the following number of buildings expected for each of the structural types:
5.0827 for type `#/MR/LWAL+CDL/H:2//RES////////`.
1.6930 for type `#/CR/LFINF+CDL/HBET:3-5//RES////////`.
1.6930 for type `#/CR/LFINF+CDN/HBET:3-5//RES////////`.
1.1276 for type `#/MR/LWAL+CDL/HBET:3-5//COM////////`.
1.1027 for type `#/S/LFM+CDL/H:1//RES////////`.
0.7775 for type `#/MR/LWAL+CDL/H:1//COM////////`.
0.5380 for type `#/CR/LFM+CDL/HBET:3-5//COM////////`.
0.4411 for type `#/S/LWAL+CDL/H:1//RES////////`.
0.3307 for type `#/MR/LWAL+CDL/H:2//COM////////`.
0.2691 for type `#/S/LFM+CDN/H:1//IND////////`.
0.2487 for type `#/CR/LFM+CDL/H:1//COM////////`.
0.2205 for type `#/W/LFM+CDL/H:1//RES////////`.
0.2153 for type `#/CR+PC/LFM+CDN/H:1//IND////////`.
0.1296 for type `#/S/LFM+CDL/HBET:3-5//COM////////`.
0.1128 for type `#/CR/LFM+CDL/H:2//COM////////`.
0.0875 for type `#/S+SL/LFM+CDN/H:1//IND////////`.
0.0499 for type `#/S/LWAL+CDL/HBET:3-5//COM////////`.
0.0444 for type `#/S/LFM+CDL/H:1//COM////////`.
0.0434 for type `#/W/LFM+CDL/H:1//COM////////`.
0.0404 for type `#/CR/LFM+CDN/H:1//IND////////`.
0.0336 for type `#/S/LFBR+CDN/H:1//IND////////`.
0.0315 for type `#/W/LFM+CDL/H:2//COM////////`.
0.0269 for type `#/MUR/LWAL+CDN/H:1//IND////////`.
0.0174 for type `#/S/LWAL+CDL/H:1//COM////////`.
0.0048 for type `#/S/LFM+CDL/H:2//COM////////`.
0.0019 for type `#/S/LWAL+CDL/H:2//COM////////`.
In total, there are 14.3639267692 buildings expected in the tile.
Spatial operations¶
Usually, one does not care too much about a single tile or building, but a geographical area. We can query the same information, but for a larger area.
bounding_box = f"1.4,42.4,1.6,42.6"
db.cursor.execute(f"""
SELECT taxonomy_string, SUM(number) as number, SUM(night) AS night, SUM(structural) AS structural
FROM entity
INNER JOIN asset ON entity.id = asset.entity_id
INNER JOIN taxonomy ON asset.taxonomy_id = taxonomy.id
WHERE ST_Intersects(entity.geom, BuildMbr({bounding_box}, 4326))
GROUP BY taxonomy_string
ORDER BY number DESC
""")
total_number_of_buildings = 0
i = 0
print("Probabilistically, in Quadkey `{quadkey}`, there are the following number of buildings expected for each of the structural types:")
for taxonomy_string, number, night, structural in db.cursor:
if i < 5:
print(f"{number:.4f} for type `{taxonomy_string}`,\nwith a total night population of {night:.2f}\n")
total_number_of_buildings += number
i += 1
if i >= 5:
print("more...")
print()
print(f"In total, there are {total_number_of_buildings} buildings expected in the area.")
Probabilistically, in Quadkey `{quadkey}`, there are the following number of buildings expected for each of the structural types:
3743.8335 for type `#/MR/LWAL+CDL/H:2//RES////////`,
with a total night population of 20581.31
1264.7322 for type `#/CR/LFINF+CDN/HBET:3-5//RES////////`,
with a total night population of 21653.49
1264.7322 for type `#/CR/LFINF+CDL/HBET:3-5//RES////////`,
with a total night population of 21653.49
805.6699 for type `#/S/LFM+CDL/H:1//RES////////`,
with a total night population of 3637.37
410.7774 for type `#/MR/LWAL+CDL/H:1//COM////////`,
with a total night population of 71.58
more...
In total, there are 9416.438219204278 buildings expected in the area.
bounding_box = f"1.4,42.4,1.6,42.6"
sql_query = f"""
SELECT id, ST_AsText(geom), SUM(asset.night), SUM(asset.structural), SUM(asset.number)
FROM entity
INNER JOIN asset ON entity.id = asset.entity_id
GROUP BY id
"""
db.cursor.execute(sql_query)
df_dict = {"id": [], "geometry": [], "population": [], "structural_value": [], "n_buildings": []}
for idx, geom, population,structural_value,n_buildings in db.cursor:
df_dict["id"].append(idx)
df_dict["geometry"].append(from_wkt(geom))
df_dict["population"].append(population)
df_dict["structural_value"].append(structural_value)
df_dict["n_buildings"].append(n_buildings)
gdf = geopandas.GeoDataFrame(df_dict, crs='EPSG:4326')
gdf.explore("n_buildings")
db.close()