Exploratory notebook for landslides NbS

.gitignore

@@ -10,6 +10,7 @@ build
 # Data
 data/*
 results/*
+notebooks/*.geojson
 notebooks/*.gpkg
 notebooks/*.png
 notebooks/*.qgz

notebooks/explore-landslides.ipynb

@@ -13,6 +13,7 @@
     "import numpy\n",
     "import pandas\n",
     "import rasterio\n",
+    "import orjson as json\n",
     "from snail.intersection import get_cell_indices\n",
     "from tqdm.auto import tqdm\n",
     "\n",
@@ -165,6 +166,180 @@
     "    ) as dataset:\n",
     "    dataset.write(data, 1)"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Demo options"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "list(Path(\".\").glob(\"*.gpkg\"))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "basins = geopandas.read_file('selection__hybas_lev12_v1c.gpkg')\n",
+    "minx, miny, maxx, maxy = basins.centroid.total_bounds"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "all_basins = geopandas.read_parquet(\"../results/input/hydrobasins/hybas_lev12_v1c.geoparquet\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "basins = all_basins.cx[minx:maxx, miny:maxy]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "splits = geopandas.read_file('intersection__ls_nbs_current__split_ead.gpkg').drop(columns=\"DN\").drop_duplicates()\n",
+    "splits.shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "hazard_cols = [c for c in splits.columns if \"EAD\" in c]\n",
+    "hazard_cols"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "options = geopandas.read_file('extract__ls_nbs_current.gpkg').query(\"DN != 0\").reset_index(drop=True)  # in case we polygonised some holes, drop DN==0\n",
+    "print(options.shape)\n",
+    "options[\"feature_id\"] = range(len(options))\n",
+    "options[\"option_landuse\"] = options[\"DN\"].map({1: \"crops\", 2: \"other\"})\n",
+    "options = (\n",
+    "    options\n",
+    "    .sjoin(basins[[\"HYBAS_ID\", \"geometry\"]], predicate=\"intersects\")\n",
+    "    .drop(columns=[\"index_right\", \"DN\"])\n",
+    "    .drop_duplicates(subset=\"feature_id\")\n",
+    "    .sort_values(by=\"feature_id\")\n",
+    ")\n",
+    "buf_geom = options.geometry.to_crs(\"ESRI:54009\").buffer(500).to_crs(\"EPSG:4326\")\n",
+    "buf = options.copy()\n",
+    "buf.geometry = buf_geom"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "options"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "options.iloc[0]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "joined = buf.sjoin(splits, how=\"left\")\n",
+    "options_damages = joined[[\"feature_id\"] + hazard_cols].groupby(\"feature_id\").sum()\n",
+    "\n",
+    "def json_nanlist(series):\n",
+    "    \"\"\"Aggregation function, returns a JSON-encoded list of unique values\n",
+    "    from the series (excluding NaNs) decoded to str (assuming json.dumps\n",
+    "    is returning bytes, as in orjson)\n",
+    "    \"\"\"\n",
+    "    return json.dumps(list(series.dropna().unique())).decode()\n",
+    "\n",
+    "options_ids = joined[[\"feature_id\", \"id\"]].groupby(\"feature_id\").agg({\"id\": json_nanlist}).rename(columns={\"id\": \"feature_ids\"})\n",
+    "\n",
+    "options_ead = (\n",
+    "    options\n",
+    "    .set_index(\"feature_id\")\n",
+    "    [[\"option_landuse\", \"HYBAS_ID\", \"geometry\"]]\n",
+    "    .join(options_ids)\n",
+    "    .join(options_damages)\n",
+    ").reset_index()\n",
+    "\n",
+    "options_ead[\"area_m2\"] = options_ead.geometry.to_crs(\"ESRI:54009\").area\n",
+    "options_ead[\"area_ha\"] = options_ead[\"area_m2\"] * 1e-4\n",
+    "# On average, self-planted projects cost the landowner 40-80p per tree plus VAT,\n",
+    "# or £1-£1.80 per tree plus VAT where planted by a Trust-arranged contractor. A\n",
+    "# 1ha site at the recommended 1,000-1,600 trees per hectare would therefore cost\n",
+    "# the landowner around £900-2,000.\n",
+    "# https://www.woodlandtrust.org.uk/plant-trees/trees-for-landowners-and-farmers/morewoods/#:~:text=On%20average%2C%20self%2Dplanted%20projects,landowner%20around%20%C2%A3900%2D2%2C000.\n",
+    "# As of 2024-11: 900 GBP is ~1160 USD\n",
+    "#               2000         2580\n",
+    "options_ead[\"cost_usd_amin\"] = options_ead[\"area_ha\"] * 1160  # demo cost_per_ha\n",
+    "options_ead[\"cost_usd_mean\"] = options_ead[\"area_ha\"] * ((1160 + 2580) / 2)  # demo cost_per_ha\n",
+    "options_ead[\"cost_usd_amax\"] = options_ead[\"area_ha\"] * 2580  # demo cost_per_ha\n",
+    "# On average, one hectare of native broadleaf woodland will store 300 - 350\n",
+    "# tonnes of carbon over a 100-year period.\n",
+    "# https://www.woodlandtrust.org.uk/plant-trees/woodland-carbon-farmers-and-landowners/#:~:text=How%20much%20carbon%20do%20trees,over%20a%20100%2Dyear%20period.\n",
+    "options_ead[\"carbon_capture_t_amin\"] = options_ead[\"area_ha\"] * 300  # demo carbon benefit\n",
+    "options_ead[\"carbon_capture_t_mean\"] = options_ead[\"area_ha\"] * 325  # demo carbon benefit\n",
+    "options_ead[\"carbon_capture_t_amax\"] = options_ead[\"area_ha\"] * 350  # demo carbon benefit\n",
+    "\n",
+    "\n",
+    "options_ead.rename(columns={\n",
+    "    \"hazard-ls_eq_tiled__road_damage_fraction_EAD\": \"hazard-ls_eq__avoided_ead_amin\",\n",
+    "    \"hazard-ls_eq_tiled__road_lower_EAD\": \"hazard-ls_eq__avoided_ead_mean\",\n",
+    "    \"hazard-ls_eq_tiled__road_upper_EAD\": \"hazard-ls_eq__avoided_ead_amax\",\n",
+    "    \"hazard-LS_RF_Median_1980-2018__road_damage_fraction_EAD\": \"hazard-ls_rf__avoided_ead_amin\",\n",
+    "    \"hazard-LS_RF_Median_1980-2018__road_lower_EAD\": \"hazard-ls_rf__avoided_ead_mean\",\n",
+    "    \"hazard-LS_RF_Median_1980-2018__road_upper_EAD\": \"hazard-ls_rf__avoided_ead_amax\",\n",
+    "    \"hazard-_landslide_sum__road_damage_fraction_EAD\": \"hazard-ls_sum__avoided_ead_amin\",\n",
+    "    \"hazard-_landslide_sum__road_lower_EAD\": \"hazard-ls_sum__avoided_ead_mean\",\n",
+    "    \"hazard-_landslide_sum__road_upper_EAD\": \"hazard-ls_sum__avoided_ead_amax\",\n",
+    "}, inplace=True)\n",
+    "options_ead.iloc[11]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "options_ead.to_file(\"joined__ls_nbs_current__split_ead.geojson\")\n",
+    "options_ead.to_file(\"joined__ls_nbs_current__split_ead.gpkg\")"
+   ]
   }
  ],
  "metadata": {