Flood Exposure Notebook

Explore S2S Metadata and API

This notebook walks through an example that explores the Space2Stats Metadata, and then uses the API to fetch flood data for various provinces within a country.

# !pip install space2stats-client matplotlib contextily 

import pandas as pd
import geopandas as gpd
from shapely import from_geojson
import matplotlib.pyplot as plt # pip install matplotlib contextily
import contextily as ctx
from space2stats_client import Space2StatsClient

# Initialize the client
client = Space2StatsClient(verify_ssl=False)  # Set verify_ssl=False if you encounter SSL issues

# Supress InsecureRequestWarning
import urllib3
urllib3.disable_warnings(urllib3.exceptions.InsecureRequestWarning)

Query Metadata

Each dataset in Space2Stats is stored as a STAC item. Metadata for each item can be explored through the following browser.

The get_topics function retrieves a list of dictionaries with key details for each dataset. The keys in each dictionary are the item ids.

topics = client.get_topics()
pd.options.display.max_colwidth = None
topics

We can extract additional metadata like fields and descriptions using the item id.

properties = client.get_properties("flood_exposure_15cm_1in100")
properties

Alternatively, we can also explore the fields avaialble via the API fields endpoint:

fields = client.get_fields()
fields

Extract H3 Data

Let’s work with the subset of fields from the flood exposure item: ['pop', 'pop_flood', 'pop_flood_pct']

flood_vars = ['pop', 'pop_flood', 'pop_flood_pct']

We will define our AOIs by fetching admin boundaries from the GeoBoundaries project.

# Try fetching the boundaries
ISO3 = "UGA" # Uganda
ADM = "ADM2" # Level 2 administrative boundaries
adm_boundaries = client.fetch_admin_boundaries(ISO3, ADM)
adm_boundaries.plot()

client.get_summary?

Run API Calls

df = client.get_summary(
    gdf=adm_boundaries, 
    spatial_join_method="centroid", 
    fields=flood_vars,
    geometry="polygon"
)

pd.reset_option('display.max_colwidth')
df.head()

Check that there are no duplicate hexagon ids

df['hex_id'].duplicated().sum()

Convert geometry column from geojson into shapely polygons

df["geometry"] = df["geometry"].apply(lambda geom: from_geojson(geom))
gdf = gpd.GeoDataFrame(df, geometry="geometry", crs="EPSG:4326")
gdf.head()

Map hexagon data

fig, ax = plt.subplots(1, 1, figsize=(10, 10))
gdf.plot(ax=ax, column="pop_flood_pct", 
         legend=True, cmap="Reds", alpha=0.75, 
         scheme="equal_interval", k=5, 
         legend_kwds=dict(title='% Pop. Exposed', fmt="{:.0%}"),
         linewidth=0)
ctx.add_basemap(ax, source=ctx.providers.Esri.WorldPhysical, crs=3857)
plt.axis("off")
plt.show()

Extract Admin Summaries

adm_boundaries_zs = client.get_aggregate(
    gdf=adm_boundaries, 
    spatial_join_method="centroid", 
    fields=['pop', 'pop_flood'], 
    aggregation_type="sum"
)

adm_boundaries_zs.head()

Recalculate share of population exposed with aggregate data

adm_boundaries_zs.loc[:, "pop_flood_pct"] = adm_boundaries_zs["pop_flood"] / adm_boundaries_zs["pop"]

fig, ax = plt.subplots(1, 1, figsize=(10, 10))
adm_boundaries_zs.plot(
    ax=ax, column="pop_flood_pct", legend=True, 
    cmap="Reds", scheme="natural_breaks", 
    k=5, legend_kwds=dict(title='% Pop. Exposed', fmt="{:.0%}"),
    linewidth=0.2, edgecolor='black')
ctx.add_basemap(ax, source=ctx.providers.Esri.WorldPhysical, crs='EPSG:4326')
plt.axis("off")
plt.show()

List top 20 provinces by population exposed

table = adm_boundaries_zs.sort_values('pop_flood_pct', ascending=False).head(20)[['shapeName', 'pop_flood', 'pop_flood_pct']].rename(
    columns={
        'shapeName': 'Province'
        })
table.loc[:, "Population Exposed"] = table.loc[:, "pop_flood"].apply(lambda x: f"{x:,.0f}")
table.loc[:, "Population Exposed (%)"] = table.loc[:, "pop_flood_pct"].apply(lambda x: f"{x:.2%}")
table.reset_index(drop=True, inplace=True)
display(table[['Province', 'Population Exposed', 'Population Exposed (%)']])