{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# Annex 01 - Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import datetime as dt\n",
    "import os\n",
    "\n",
    "import contextily as cx\n",
    "import geopandas as gpd\n",
    "import numpy as np\n",
    "import osmnx as ox\n",
    "import pandas as pd\n",
    "import pooch\n",
    "import pystac_client\n",
    "import rasterio as rio\n",
    "import seaborn as sns\n",
    "import swisstopopy as stp\n",
    "import zenodo_client\n",
    "from meteora import utils\n",
    "from meteora.clients import AWELClient, MeteoSwissClient, NetatmoClient\n",
    "from odc import stac as odc_stac\n",
    "from rasterio import mask, plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# region data\n",
    "region = \"Zurich\"\n",
    "crs = \"epsg:2056\"  # projected CRS in meters\n",
    "buffer_dist = 2000  # in meters\n",
    "\n",
    "# Netatmo auth\n",
    "client_id = os.getenv(\"NETATMO_CLIENT_ID\", \"\")\n",
    "client_secret = os.getenv(\"NETATMO_CLIENT_SECRET\", \"\")\n",
    "\n",
    "# for heatwave detection\n",
    "year = 2022\n",
    "start_month = 6\n",
    "end_month = 8\n",
    "heatwave_t_threshold = 25\n",
    "heatwave_n_consecutive_days = 3\n",
    "station_agg_func = \"mean\"\n",
    "inter_station_agg_func = \"max\"\n",
    "\n",
    "# target temperature data\n",
    "target_hour = 21  # 9 PM\n",
    "\n",
    "# species richness dataset\n",
    "species_richness_filepath = (\n",
    "    \"https://www.envidat.ch/dataset/0f81429a-204e-4720-929e-5b251e9ab1c4/resource/\"\n",
    "    \"f3abc8d4-c92e-4a57-aea4-39ba8c93158d/download/dataset_5_species_richness_all.txt\"\n",
    ")\n",
    "\n",
    "# land use/land cover (LULC) survey\n",
    "# sls_column = \"LU18_4\"\n",
    "lulc_url = (\n",
    "    \"https://access.yareta.unige.ch/access/metadata/\"\n",
    "    \"6ab4b715-904f-4cb9-961c-6a25b4c1116b/data/8744b6b3-777e-498b-b4ae-d082fd0ebb30/\"\n",
    "    \"download\"\n",
    ")\n",
    "\n",
    "# TesseraCH embeddings (Swiss Data Cube STAC)\n",
    "tessera_stac_url = \"https://explorer.swissdatacube.org/stac\"\n",
    "tessera_stac_crs = \"epsg:4326\"\n",
    "tessera_collection = \"tesserach\"\n",
    "tessera_collection_crs = \"epsg:2056\"\n",
    "tessera_year = 2024\n",
    "tessera_indexes = range(1, 17)\n",
    "tessera_search_limit = 200\n",
    "tessera_resolution = 10\n",
    "\n",
    "# output files\n",
    "dst_study_area_filepath = \"data/study-area.gpkg\"\n",
    "dst_stations_filepath = \"data/stations.gpkg\"\n",
    "dst_netatmo_cws_filepath = \"data/netatmo-cws.gpkg\"\n",
    "dst_ts_df_filepath = \"data/heatwaves-2022.csv\"\n",
    "dst_bird_richness_filepath = \"data/bird-richness.gpkg\"\n",
    "dst_lulc_filepath = \"data/lulc.tif\"\n",
    "dst_buildings_gdf_filepath = \"data/buildings.gpkg\"\n",
    "dst_dem_filepath = \"data/dem.tif\"\n",
    "dst_tree_canopy_filepath = \"data/tree-canopy.tif\"\n",
    "dst_tessera_filepath = \"data/tesserach-zurich.tif\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# buffer around the region so that we can compute focal features\n",
    "region_buffer_gser = ox.geocode_to_gdf(region).to_crs(crs).buffer(buffer_dist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save study area to a file\n",
    "region_buffer_gser.to_file(dst_study_area_filepath)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5",
   "metadata": {},
   "source": [
    "## Meteoroloical data\n",
    "\n",
    "### Meteorological stations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get meteorological stations for the (non-buffered) region\n",
    "client_dict = {\n",
    "    \"AWEL\": AWELClient(region),\n",
    "    \"MeteoSwiss\": MeteoSwissClient(region),\n",
    "}\n",
    "# rename elevation column to have a common label\n",
    "ELEV_COL_RENAME_DICT = {\n",
    "    \"masl\": \"elevation\",  # for AWEL\n",
    "    \"station_height_masl\": \"elevation\",  # for MeteoSwiss\n",
    "}\n",
    "stations_gdf = pd.concat(\n",
    "    [\n",
    "        client.stations_gdf.to_crs(crs)\n",
    "        .rename(columns=ELEV_COL_RENAME_DICT)[[\"elevation\", \"geometry\"]]\n",
    "        .assign(source=source)\n",
    "        for source, client in client_dict.items()\n",
    "    ]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot them on a map\n",
    "ax = stations_gdf.plot(\"source\", legend=True)\n",
    "cx.add_basemap(ax, crs=stations_gdf.crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save stations to file\n",
    "stations_gdf.to_file(dst_stations_filepath)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9",
   "metadata": {},
   "source": [
    "### Citizen weather stations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [],
   "source": [
    "netatmo_cws_gdf = (\n",
    "    NetatmoClient(region, client_id, client_secret)\n",
    "    .stations_gdf.rename(columns={\"altitude\": \"elevation\"})\n",
    "    .to_crs(stations_gdf.crs)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot them on a map\n",
    "ax = pd.concat([netatmo_cws_gdf.assign(source=\"Netatmo\"), stations_gdf]).plot(\n",
    "    \"source\", legend=True\n",
    ")\n",
    "cx.add_basemap(ax, crs=stations_gdf.crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "netatmo_cws_gdf.to_file(dst_netatmo_cws_filepath)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13",
   "metadata": {},
   "source": [
    "### Temperature data\n",
    "\n",
    "We will get the air temperature data [for all the heatwaves](https://www.meteoswiss.admin.ch/weather/weather-and-climate-from-a-to-z/heat-warnings.html) of the target year."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the JJA time series of temperature measurements, convert it to wide format, drop\n",
    "# stations with no data and resample to hourly mean\n",
    "ts_df = (\n",
    "    utils.long_to_wide(\n",
    "        pd.concat(\n",
    "            [\n",
    "                client.get_ts_df(\n",
    "                    \"temperature\",\n",
    "                    dt.date(year, start_month, 1),\n",
    "                    # get latest moment of the latest day of the month\n",
    "                    # ACHTUNG: this won't work if `end_month` is 12 (see commented code\n",
    "                    # below)\n",
    "                    dt.datetime.combine(\n",
    "                        dt.date(year, end_month + 1, 1) - dt.timedelta(days=1),\n",
    "                        dt.time.max,\n",
    "                    ),\n",
    "                )\n",
    "                for client in client_dict.values()\n",
    "            ]\n",
    "        )\n",
    "    )\n",
    "    .dropna(how=\"all\", axis=\"columns\")\n",
    "    .resample(\"h\")\n",
    "    .mean()\n",
    ")\n",
    "ts_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get only heatwave data\n",
    "heatwave_ts_df = utils.get_heatwave_ts_df(\n",
    "    ts_df,\n",
    "    heatwave_t_threshold=heatwave_t_threshold,\n",
    "    heatwave_n_consecutive_days=heatwave_n_consecutive_days,\n",
    "    station_agg_func=station_agg_func,\n",
    "    inter_station_agg_func=inter_station_agg_func,\n",
    ")\n",
    "heatwave_ts_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot hourly time series\n",
    "plot_ts_df = heatwave_ts_df.stack().rename(\"temperature\").reset_index()\n",
    "sns.lineplot(\n",
    "    data=plot_ts_df.assign(**{\"hour\": plot_ts_df[\"time\"].dt.hour}),\n",
    "    x=\"hour\",\n",
    "    y=\"temperature\",\n",
    "    hue=\"heatwave\",\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save temperature time series to file (only target hour)\n",
    "heatwave_ts_df.droplevel(\"heatwave\")[\n",
    "    heatwave_ts_df.index.get_level_values(\"time\").hour == target_hour\n",
    "].to_csv(dst_ts_df_filepath)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18",
   "metadata": {},
   "source": [
    "## Bird richness data\n",
    "\n",
    "From <cite data-cite=\"casanelles2021applying\"/>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19",
   "metadata": {},
   "outputs": [],
   "source": [
    "bird_gdf = pd.read_csv(species_richness_filepath, sep=\"\\t\")\n",
    "bird_gdf = bird_gdf[bird_gdf[\"GROUP\"] == \"Birds\"]\n",
    "bird_gdf = gpd.GeoDataFrame(\n",
    "    bird_gdf, geometry=gpd.points_from_xy(bird_gdf[\"X\"], bird_gdf[\"Y\"]), crs=\"epsg:2056\"\n",
    ")\n",
    "ax = bird_gdf.plot(\"n.species\", legend=True, legend_kwds={\"label\": \"N. species\"})\n",
    "cx.add_basemap(ax, crs=bird_gdf.crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save to a file\n",
    "bird_gdf.to_file(dst_bird_richness_filepath)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21",
   "metadata": {},
   "source": [
    "## Land use data\n",
    "\n",
    "Version of the [Swiss land use statistics](https://www.bfs.admin.ch/bfs/en/home/dienstleistungen/geostat/geodaten-bundesstatistik/boden-nutzung-bedeckung-eignung/arealstatistik-schweiz.html) downscaled to 25 m resolution by <cite data-cite=\"giuliani2022downscaling\"/>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get Swiss LULC data and clip it to the target (buffered) region only\n",
    "lulc_filepath = pooch.retrieve(lulc_url, known_hash=None)\n",
    "with rio.open(lulc_filepath) as src:\n",
    "    dst_img, dst_transform = mask.mask(src, [region_buffer_gser.iloc[0]], crop=True)\n",
    "    dst_meta = src.meta.copy()\n",
    "\n",
    "# replace \"nodata\" with 0, the \"nodata\" value in the Swiss land use statistics\n",
    "dst_nodata = 0  # src.nodata\n",
    "dst_meta.update(\n",
    "    {\n",
    "        \"driver\": \"GTiff\",\n",
    "        \"nodata\": dst_nodata,\n",
    "        \"height\": dst_img.shape[1],\n",
    "        \"width\": dst_img.shape[2],\n",
    "        \"transform\": dst_transform,\n",
    "    }\n",
    ")\n",
    "\n",
    "# reclassify to 4 main domains only - ACHTUNG: ugly hardcoded procedure\n",
    "dst_img = np.digitize(\n",
    "    dst_img, [1, 3, 15, 24, 31, 37, 41, 42, 45, 50, 57, 58, 61, 62, 64, 69, 72]\n",
    ")\n",
    "\n",
    "dst_img = np.where(dst_img < 17, dst_img, dst_nodata)\n",
    "with rio.open(dst_lulc_filepath, \"w\", **dst_meta) as dst:\n",
    "    dst.write(dst_img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the region data from the target LULC survey\n",
    "with rio.open(dst_lulc_filepath) as src:\n",
    "    plot.show(src)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24",
   "metadata": {},
   "source": [
    "## Buildings with heights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get building footprints with heights (from Swisstopo) for the (buffered) region\n",
    "buildings_gdf = stp.get_bldg_gdf(region_buffer_gser)\n",
    "# compute area and add it as a column (note that building_gdf is already projected)\n",
    "buildings_gdf = buildings_gdf.assign(**{\"area\": buildings_gdf.area})\n",
    "# add a \"volume\" column\n",
    "buildings_gdf[\"volume\"] = buildings_gdf[\"area\"] * buildings_gdf[\"height\"]\n",
    "# drop the \"element\" index level (index by osmid only)\n",
    "buildings_gdf = buildings_gdf.droplevel(\"element\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot buildings with heights\n",
    "ax = buildings_gdf.plot(\n",
    "    column=\"height\",\n",
    "    cmap=\"viridis\",\n",
    "    legend=True,\n",
    "    legend_kwds={\"label\": \"Building height (m)\"},\n",
    ")\n",
    "cx.add_basemap(ax, crs=buildings_gdf.crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save to a file\n",
    "buildings_gdf.to_file(dst_buildings_gdf_filepath)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28",
   "metadata": {},
   "source": [
    "## Digital Elevation model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get DEM raster (from Swisstopo) for the (buffered) region and save to a file\n",
    "stp.get_dem_raster(region_buffer_gser, dst_dem_filepath, alti3d_res=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# # plot DEM raster\n",
    "# fig, ax = plt.subplots()\n",
    "# with rio.open(dst_dem_filepath) as src:\n",
    "#     retted = plot.show(src, ax=ax)\n",
    "# fig.colorbar(retted.get_images()[0], ax=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31",
   "metadata": {},
   "source": [
    "## Tree canopy cover"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get tree canopy cover raster (from Swisstopo) for the (buffered) region and save to a\n",
    "# file\n",
    "stp.get_tree_canopy_raster(\n",
    "    region_buffer_gser, dst_tree_canopy_filepath, pooch_retrieve_kwargs={\"path\": \"las\"}\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot tree canopy map\n",
    "with rio.open(dst_tree_canopy_filepath) as src:\n",
    "    retted = plot.show(src)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34",
   "metadata": {},
   "source": [
    "## TesseraCH embeddings\n",
    "\n",
    "We can also prepare a local multi-band TesseraCH raster for the buffered study area.\n",
    "TesseraCH is the Swiss Data Cube embedding product derived from Sentinel-2 imagery and\n",
    "distributed as one GeoTIFF per embedding dimension in `EPSG:2056` at 10 m resolution.\n",
    "\n",
    "Here we use `pystac-client` to search the official Swiss Data Cube STAC API and\n",
    "`odc-stac` to load and mosaic the matching TesseraCH items directly into a single\n",
    "xarray object, which we then clip to the buffered study area and save as a local\n",
    "multi-band GeoTIFF.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35",
   "metadata": {},
   "outputs": [],
   "source": [
    "# search the STAC API for TesseraCH items intersecting the buffered study area, then\n",
    "# load the requested dimensions with odc-stac, clip the mosaic, and save it as a local\n",
    "# multi-band raster\n",
    "region_buffer_stac_gser = region_buffer_gser.to_crs(tessera_stac_crs)\n",
    "tessera_band_names = [f\"Dim{index}\" for index in tessera_indexes]\n",
    "\n",
    "stac_client = pystac_client.Client.open(tessera_stac_url)\n",
    "search = stac_client.search(\n",
    "    collections=[tessera_collection],\n",
    "    intersects=region_buffer_stac_gser.iloc[0],\n",
    "    datetime=f\"{tessera_year}-01-01T00:00:00Z/{tessera_year}-12-31T23:59:59Z\",\n",
    "    limit=tessera_search_limit,\n",
    ")\n",
    "items = search.item_collection()\n",
    "print(f\"Found {len(items)} TesseraCH STAC item(s).\")\n",
    "\n",
    "tessera_ds = odc_stac.load(\n",
    "    items,\n",
    "    bands=tessera_band_names,\n",
    "    bbox=region_buffer_stac_gser.total_bounds.tolist(),\n",
    "    crs=crs,\n",
    "    resolution=tessera_resolution,\n",
    "    groupby=\"solar_day\",\n",
    ")\n",
    "\n",
    "if \"time\" in tessera_ds.dims:\n",
    "    # select first (time) snapshot\n",
    "    tessera_ds = tessera_ds.isel(time=0, drop=True)\n",
    "\n",
    "tessera_da = tessera_ds.to_array(dim=\"band\").assign_coords(band=tessera_band_names)\n",
    "tessera_da"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36",
   "metadata": {},
   "outputs": [],
   "source": [
    "# clip to the region extent (collection CRS is already EPSG:2056) and save to file\n",
    "tessera_clipped_da = tessera_da.rio.clip(region_buffer_gser.geometry)\n",
    "tessera_clipped_da.rio.to_raster(dst_tessera_filepath)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the first embedding dimension\n",
    "with rio.open(dst_tessera_filepath) as src:\n",
    "    plot.show(src.read(1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38",
   "metadata": {},
   "source": [
    "## Upload to zenodo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39",
   "metadata": {},
   "outputs": [],
   "source": [
    "# first upload\n",
    "# define the metadata that will be used on initial upload\n",
    "# data = zenodo_client.Metadata(\n",
    "#     title=\"focalpy user guide data\",\n",
    "#     upload_type=\"dataset\",\n",
    "#     description=\"Data used for the focalpy user guide\",\n",
    "#     creators=[\n",
    "#         zenodo_client.Creator(\n",
    "#             name=\"Bosch, Martí\",\n",
    "#             affiliation=\"University of Bern\",\n",
    "#             orcid=\"0000-0001-8735-9144\",\n",
    "#         ),\n",
    "#     ],\n",
    "# )\n",
    "# res = zenodo_client.ensure_zenodo(\n",
    "#     # this is a unique key you pick that will be used to store the numeric deposition\n",
    "#     # ID on your local system's cache\n",
    "#     key=\"focalpy\",\n",
    "#     data=data,\n",
    "#     paths=[\n",
    "#         dst_study_area_filepath,\n",
    "#         dst_stations_filepath,\n",
    "#         dst_ts_df_filepath,\n",
    "#         dst_lulc_filepath,\n",
    "#         dst_buildings_gdf_filepath,\n",
    "#         dst_dem_filepath,\n",
    "#         dst_tree_canopy_filepath,\n",
    "#         dst_tessera_filepath,\n",
    "#     ],\n",
    "#     # sandbox=True,  # remove this when you're ready to upload to real Zenodo\n",
    "# )\n",
    "# update existing repository\n",
    "DEP_ID = \"17725958\"\n",
    "res = zenodo_client.update_zenodo(\n",
    "    DEP_ID,\n",
    "    # add the list of files to update\n",
    "    [dst_tessera_filepath],\n",
    ")\n",
    "\n",
    "# show the links to the generated zenodo repository\n",
    "links = res.json()[\"links\"]\n",
    "print(f\"URL: {links['html']}\\nDOI: {links['doi']}\")"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "-all"
  },
  "kernelspec": {
   "display_name": "Python (Pixi)",
   "language": "python",
   "name": "pixi-kernel-python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.9"
  },
  "pixi-kernel": {
   "environment": "user-guide"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}