From f76a435c5e00d55e793a833715f6d615dbf9dee8 Mon Sep 17 00:00:00 2001 From: Tom Russell Date: Mon, 11 Nov 2024 14:44:31 +0000 Subject: [PATCH] Exploratory notebook for landslides NbS --- .gitignore | 1 + notebooks/explore-landslides.ipynb | 175 +++++++++++++++++++++++++++++ 2 files changed, 176 insertions(+) diff --git a/.gitignore b/.gitignore index 9f61de52..fe3e0923 100644 --- a/.gitignore +++ b/.gitignore @@ -10,6 +10,7 @@ build # Data data/* results/* +notebooks/*.geojson notebooks/*.gpkg notebooks/*.png notebooks/*.qgz diff --git a/notebooks/explore-landslides.ipynb b/notebooks/explore-landslides.ipynb index 6b7d9caa..334b22e1 100644 --- a/notebooks/explore-landslides.ipynb +++ b/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": {