Notebook for landslide results exploration

.gitignore

@@ -10,6 +10,11 @@ build
 # Data
 data/*
 results/*
+notebooks/*.gpkg
+notebooks/*.png
+notebooks/*.qgz
+notebooks/*.tif
+notebooks/*.xml
 *.bak
 
 # snakemake

notebooks/explore-landslides.ipynb

@@ -0,0 +1,191 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from functools import partial\n",
+    "from pathlib import Path\n",
+    "\n",
+    "import geopandas\n",
+    "import numpy\n",
+    "import pandas\n",
+    "import rasterio\n",
+    "from snail.intersection import get_cell_indices\n",
+    "from tqdm.auto import tqdm\n",
+    "\n",
+    "tqdm.pandas()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "base_dir = Path(\"../results\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "damage_cost = geopandas.read_parquet(base_dir/\"africa-latest_filter-road-tertiary/hazard-landslide-arup/damage_EAD_and_cost_per_trigger.geoparquet\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "damage_cost.to_file(\"africa-latest_damage_EAD_and_cost_per_trigger.gpkg\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "splits_dir = base_dir / 'direct_damages/africa-latest_filter-road-tertiary/hazard-landslide-arup/split_EAD_and_cost_per_trigger'\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dfs = [geopandas.read_parquet(fname) for fname in splits_dir.glob(\"*.geoparquet\")]\n",
+    "\n",
+    "splits = pandas.concat(dfs)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "splits.to_file(\"split_EAD_and_cost_per_trigger.gpkg\", driver=\"GPKG\", engine=\"pyogrio\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "splits.columns"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "with rasterio.open(base_dir/\"input/hazard-landslide-arup/africa-latest/ls_eq_tiled.tif\") as src:\n",
+    "    pass\n",
+    "\n",
+    "src.width, src.height, src.transform"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "raster_height = src.height\n",
+    "raster_width = src.width\n",
+    "raster_transform = src.transform\n",
+    "data_col = 'hazard-_landslide_sum__road_damage_fraction_EAD'\n",
+    "\n",
+    "data = numpy.zeros((raster_height, raster_width))\n",
+    "\n",
+    "splits_df = splits\n",
+    "\n",
+    "def get_cell_indices_kw(geom, raster_height=0, raster_width=0, raster_transform=None):\n",
+    "    return get_cell_indices(geom, raster_height, raster_width, raster_transform)\n",
+    "\n",
+    "cell_indices_of_split_geometry = partial(get_cell_indices_kw, raster_height=src.height, raster_width=src.width, raster_transform=list(src.transform))\n",
+    "\n",
+    "splits_df = pandas.concat(\n",
+    "    [\n",
+    "        splits_df[[data_col]],\n",
+    "        splits_df.geometry.progress_apply(cell_indices_of_split_geometry)\n",
+    "    ],\n",
+    "    axis=\"columns\"\n",
+    ").rename(columns={\"geometry\":\"cell_index\"})"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "value_per_cell = splits_df[['cell_index', data_col]].groupby('cell_index').sum()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for _, item in value_per_cell.reset_index().iterrows():\n",
+    "    col, row = item.cell_index\n",
+    "    data[int(row), int(col)] = item[data_col]\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "with rasterio.open(\n",
+    "        'africa-latest_damage_EAD.tif',\n",
+    "        'w',\n",
+    "        driver='GTiff',\n",
+    "        height=data.shape[0],\n",
+    "        width=data.shape[1],\n",
+    "        count=1,\n",
+    "        dtype=data.dtype,\n",
+    "        crs='+proj=latlong',\n",
+    "        transform=raster_transform,\n",
+    "        compress='lzw'\n",
+    "    ) as dataset:\n",
+    "    dataset.write(data, 1)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "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.10.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}