Annex 01 - Data

import datetime as dt
import os

import contextily as cx
import geopandas as gpd
import numpy as np
import osmnx as ox
import pandas as pd
import pooch
import pystac_client
import rasterio as rio
import seaborn as sns
import swisstopopy as stp
import zenodo_client
from meteora import utils
from meteora.clients import AWELClient, MeteoSwissClient, NetatmoClient
from odc import stac as odc_stac
from rasterio import mask, plot

# region data
region = "Zurich"
crs = "epsg:2056"  # projected CRS in meters
buffer_dist = 2000  # in meters

# Netatmo auth
client_id = os.getenv("NETATMO_CLIENT_ID", "")
client_secret = os.getenv("NETATMO_CLIENT_SECRET", "")

# for heatwave detection
year = 2022
start_month = 6
end_month = 8
heatwave_t_threshold = 25
heatwave_n_consecutive_days = 3
station_agg_func = "mean"
inter_station_agg_func = "max"

# target temperature data
target_hour = 21  # 9 PM

# species richness dataset
species_richness_filepath = (
    "https://www.envidat.ch/dataset/0f81429a-204e-4720-929e-5b251e9ab1c4/resource/"
    "f3abc8d4-c92e-4a57-aea4-39ba8c93158d/download/dataset_5_species_richness_all.txt"
)

# land use/land cover (LULC) survey
# sls_column = "LU18_4"
lulc_url = (
    "https://access.yareta.unige.ch/access/metadata/"
    "6ab4b715-904f-4cb9-961c-6a25b4c1116b/data/8744b6b3-777e-498b-b4ae-d082fd0ebb30/"
    "download"
)

# TesseraCH embeddings (Swiss Data Cube STAC)
tessera_stac_url = "https://explorer.swissdatacube.org/stac"
tessera_stac_crs = "epsg:4326"
tessera_collection = "tesserach"
tessera_collection_crs = "epsg:2056"
tessera_year = 2024
tessera_indexes = range(1, 17)
tessera_search_limit = 200
tessera_resolution = 10

# output files
dst_study_area_filepath = "data/study-area.gpkg"
dst_stations_filepath = "data/stations.gpkg"
dst_netatmo_cws_filepath = "data/netatmo-cws.gpkg"
dst_ts_df_filepath = "data/heatwaves-2022.csv"
dst_bird_richness_filepath = "data/bird-richness.gpkg"
dst_lulc_filepath = "data/lulc.tif"
dst_buildings_gdf_filepath = "data/buildings.gpkg"
dst_dem_filepath = "data/dem.tif"
dst_tree_canopy_filepath = "data/tree-canopy.tif"
dst_tessera_filepath = "data/tesserach-zurich.tif"

# buffer around the region so that we can compute focal features
region_buffer_gser = ox.geocode_to_gdf(region).to_crs(crs).buffer(buffer_dist)

# save study area to a file
region_buffer_gser.to_file(dst_study_area_filepath)

Meteoroloical data

Meteorological stations

# get meteorological stations for the (non-buffered) region
client_dict = {
    "AWEL": AWELClient(region),
    "MeteoSwiss": MeteoSwissClient(region),
}
# rename elevation column to have a common label
ELEV_COL_RENAME_DICT = {
    "masl": "elevation",  # for AWEL
    "station_height_masl": "elevation",  # for MeteoSwiss
}
stations_gdf = pd.concat(
    [
        client.stations_gdf.to_crs(crs)
        .rename(columns=ELEV_COL_RENAME_DICT)[["elevation", "geometry"]]
        .assign(source=source)
        for source, client in client_dict.items()
    ]
)

# plot them on a map
ax = stations_gdf.plot("source", legend=True)
cx.add_basemap(ax, crs=stations_gdf.crs)

# save stations to file
stations_gdf.to_file(dst_stations_filepath)

Citizen weather stations

netatmo_cws_gdf = (
    NetatmoClient(region, client_id, client_secret)
    .stations_gdf.rename(columns={"altitude": "elevation"})
    .to_crs(stations_gdf.crs)
)

# plot them on a map
ax = pd.concat([netatmo_cws_gdf.assign(source="Netatmo"), stations_gdf]).plot(
    "source", legend=True
)
cx.add_basemap(ax, crs=stations_gdf.crs)

netatmo_cws_gdf.to_file(dst_netatmo_cws_filepath)

Temperature data

We will get the air temperature data for all the heatwaves of the target year.

# get the JJA time series of temperature measurements, convert it to wide format, drop
# stations with no data and resample to hourly mean
ts_df = (
    utils.long_to_wide(
        pd.concat(
            [
                client.get_ts_df(
                    "temperature",
                    dt.date(year, start_month, 1),
                    # get latest moment of the latest day of the month
                    # ACHTUNG: this won't work if `end_month` is 12 (see commented code
                    # below)
                    dt.datetime.combine(
                        dt.date(year, end_month + 1, 1) - dt.timedelta(days=1),
                        dt.time.max,
                    ),
                )
                for client in client_dict.values()
            ]
        )
    )
    .dropna(how="all", axis="columns")
    .resample("h")
    .mean()
)
ts_df.head()

# get only heatwave data
heatwave_ts_df = utils.get_heatwave_ts_df(
    ts_df,
    heatwave_t_threshold=heatwave_t_threshold,
    heatwave_n_consecutive_days=heatwave_n_consecutive_days,
    station_agg_func=station_agg_func,
    inter_station_agg_func=inter_station_agg_func,
)
heatwave_ts_df.head()

# plot hourly time series
plot_ts_df = heatwave_ts_df.stack().rename("temperature").reset_index()
sns.lineplot(
    data=plot_ts_df.assign(**{"hour": plot_ts_df["time"].dt.hour}),
    x="hour",
    y="temperature",
    hue="heatwave",
)

# save temperature time series to file (only target hour)
heatwave_ts_df.droplevel("heatwave")[
    heatwave_ts_df.index.get_level_values("time").hour == target_hour
].to_csv(dst_ts_df_filepath)

Bird richness data

From (Casanelles-Abella et al., 2021).

bird_gdf = pd.read_csv(species_richness_filepath, sep="\t")
bird_gdf = bird_gdf[bird_gdf["GROUP"] == "Birds"]
bird_gdf = gpd.GeoDataFrame(
    bird_gdf, geometry=gpd.points_from_xy(bird_gdf["X"], bird_gdf["Y"]), crs="epsg:2056"
)
ax = bird_gdf.plot("n.species", legend=True, legend_kwds={"label": "N. species"})
cx.add_basemap(ax, crs=bird_gdf.crs)

# save to a file
bird_gdf.to_file(dst_bird_richness_filepath)

Land use data

Version of the Swiss land use statistics downscaled to 25 m resolution by (Giuliani et al., 2022).

# get Swiss LULC data and clip it to the target (buffered) region only
lulc_filepath = pooch.retrieve(lulc_url, known_hash=None)
with rio.open(lulc_filepath) as src:
    dst_img, dst_transform = mask.mask(src, [region_buffer_gser.iloc[0]], crop=True)
    dst_meta = src.meta.copy()

# replace "nodata" with 0, the "nodata" value in the Swiss land use statistics
dst_nodata = 0  # src.nodata
dst_meta.update(
    {
        "driver": "GTiff",
        "nodata": dst_nodata,
        "height": dst_img.shape[1],
        "width": dst_img.shape[2],
        "transform": dst_transform,
    }
)

# reclassify to 4 main domains only - ACHTUNG: ugly hardcoded procedure
dst_img = np.digitize(
    dst_img, [1, 3, 15, 24, 31, 37, 41, 42, 45, 50, 57, 58, 61, 62, 64, 69, 72]
)

dst_img = np.where(dst_img < 17, dst_img, dst_nodata)
with rio.open(dst_lulc_filepath, "w", **dst_meta) as dst:
    dst.write(dst_img)

# plot the region data from the target LULC survey
with rio.open(dst_lulc_filepath) as src:
    plot.show(src)

Buildings with heights

# get building footprints with heights (from Swisstopo) for the (buffered) region
buildings_gdf = stp.get_bldg_gdf(region_buffer_gser)
# compute area and add it as a column (note that building_gdf is already projected)
buildings_gdf = buildings_gdf.assign(**{"area": buildings_gdf.area})
# add a "volume" column
buildings_gdf["volume"] = buildings_gdf["area"] * buildings_gdf["height"]
# drop the "element" index level (index by osmid only)
buildings_gdf = buildings_gdf.droplevel("element")

# plot buildings with heights
ax = buildings_gdf.plot(
    column="height",
    cmap="viridis",
    legend=True,
    legend_kwds={"label": "Building height (m)"},
)
cx.add_basemap(ax, crs=buildings_gdf.crs)

# save to a file
buildings_gdf.to_file(dst_buildings_gdf_filepath)

Digital Elevation model

# get DEM raster (from Swisstopo) for the (buffered) region and save to a file
stp.get_dem_raster(region_buffer_gser, dst_dem_filepath, alti3d_res=2)

# # plot DEM raster
# fig, ax = plt.subplots()
# with rio.open(dst_dem_filepath) as src:
#     retted = plot.show(src, ax=ax)
# fig.colorbar(retted.get_images()[0], ax=ax)

Tree canopy cover

# get tree canopy cover raster (from Swisstopo) for the (buffered) region and save to a
# file
stp.get_tree_canopy_raster(
    region_buffer_gser, dst_tree_canopy_filepath, pooch_retrieve_kwargs={"path": "las"}
)

# plot tree canopy map
with rio.open(dst_tree_canopy_filepath) as src:
    retted = plot.show(src)

TesseraCH embeddings

We can also prepare a local multi-band TesseraCH raster for the buffered study area. TesseraCH is the Swiss Data Cube embedding product derived from Sentinel-2 imagery and distributed as one GeoTIFF per embedding dimension in EPSG:2056 at 10 m resolution.

Here we use pystac-client to search the official Swiss Data Cube STAC API and odc-stac to load and mosaic the matching TesseraCH items directly into a single xarray object, which we then clip to the buffered study area and save as a local multi-band GeoTIFF.

# search the STAC API for TesseraCH items intersecting the buffered study area, then
# load the requested dimensions with odc-stac, clip the mosaic, and save it as a local
# multi-band raster
region_buffer_stac_gser = region_buffer_gser.to_crs(tessera_stac_crs)
tessera_band_names = [f"Dim{index}" for index in tessera_indexes]

stac_client = pystac_client.Client.open(tessera_stac_url)
search = stac_client.search(
    collections=[tessera_collection],
    intersects=region_buffer_stac_gser.iloc[0],
    datetime=f"{tessera_year}-01-01T00:00:00Z/{tessera_year}-12-31T23:59:59Z",
    limit=tessera_search_limit,
)
items = search.item_collection()
print(f"Found {len(items)} TesseraCH STAC item(s).")

tessera_ds = odc_stac.load(
    items,
    bands=tessera_band_names,
    bbox=region_buffer_stac_gser.total_bounds.tolist(),
    crs=crs,
    resolution=tessera_resolution,
    groupby="solar_day",
)

if "time" in tessera_ds.dims:
    # select first (time) snapshot
    tessera_ds = tessera_ds.isel(time=0, drop=True)

tessera_da = tessera_ds.to_array(dim="band").assign_coords(band=tessera_band_names)
tessera_da

# clip to the region extent (collection CRS is already EPSG:2056) and save to file
tessera_clipped_da = tessera_da.rio.clip(region_buffer_gser.geometry)
tessera_clipped_da.rio.to_raster(dst_tessera_filepath)

# plot the first embedding dimension
with rio.open(dst_tessera_filepath) as src:
    plot.show(src.read(1))

Upload to zenodo

# first upload
# define the metadata that will be used on initial upload
# data = zenodo_client.Metadata(
#     title="focalpy user guide data",
#     upload_type="dataset",
#     description="Data used for the focalpy user guide",
#     creators=[
#         zenodo_client.Creator(
#             name="Bosch, Martí",
#             affiliation="University of Bern",
#             orcid="0000-0001-8735-9144",
#         ),
#     ],
# )
# res = zenodo_client.ensure_zenodo(
#     # this is a unique key you pick that will be used to store the numeric deposition
#     # ID on your local system's cache
#     key="focalpy",
#     data=data,
#     paths=[
#         dst_study_area_filepath,
#         dst_stations_filepath,
#         dst_ts_df_filepath,
#         dst_lulc_filepath,
#         dst_buildings_gdf_filepath,
#         dst_dem_filepath,
#         dst_tree_canopy_filepath,
#         dst_tessera_filepath,
#     ],
#     # sandbox=True,  # remove this when you're ready to upload to real Zenodo
# )
# update existing repository
DEP_ID = "17725958"
res = zenodo_client.update_zenodo(
    DEP_ID,
    # add the list of files to update
    [dst_tessera_filepath],
)

# show the links to the generated zenodo repository
links = res.json()["links"]
print(f"URL: {links['html']}\nDOI: {links['doi']}")