Top-down/notebooks/Heat_stress.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7ba77d7d",
   "metadata": {},
   "source": [
    "CCDR Hazard Analysis Notebook\n",
    "Developed by M. Amadio and T. Iwanaga"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "27ad9a5b-154d-4f7a-8579-59487c6a7d76",
   "metadata": {},
   "outputs": [],
   "source": [
    "from common import *  # import necessary packages\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6738bd3b",
   "metadata": {},
   "outputs": [],
   "source": [
    "country_code_map = {\n",
    "    \"NPL\": 175,\n",
    "    \"PAK\": 188,\n",
    "    \"BGD\": 23,\n",
    "    \"GHA\": 94,\n",
    "    \"ETH\": 79,\n",
    "    \"BFA\": 42,\n",
    "    \"MLI\": 155,\n",
    "    \"NER\": 181,\n",
    "    \"TCD\": 50,\n",
    "    \"MRT\": 159\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "53ebd7a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_analysis(rb):\n",
    "    with output:\n",
    "        output.clear_output()\n",
    "        print(\"Running analysis...\")\n",
    "        rb.disabled = True\n",
    "\n",
    "    # Ensure values are valid\n",
    "    bin_seq = [w.value for w in class_edges.values()]\n",
    "\n",
    "    seq = np.all([True if bin_seq[i] < bin_seq[i+1] else False for i in range(0, len(bin_seq)-1)])\n",
    "\n",
    "    if not seq:\n",
    "        ValueError(\"Class thresholds are not sequential. Lower classes must be less than class thresholds above.\")\n",
    "        rb.disabled = False\n",
    "        return\n",
    "\n",
    "    max_bin_value = np.max(bin_seq)\n",
    "    max_haz_threshold = max_bin_value + 1e-4\n",
    "    bin_seq = bin_seq + [np.inf]\n",
    "    num_bins = len(bin_seq)\n",
    "\n",
    "    # Get user input\n",
    "    country = country_dd.value\n",
    "    exp_cat = exp_cat_dd.value\n",
    "    target_ADM = adm_dd.value\n",
    "    adm_name = target_ADM.replace('_', '')\n",
    "\n",
    "    min_haz_threshold = class_edges['class_1'].value\n",
    "\n",
    "    valid_RPs = [5, 20, 100]\n",
    "\n",
    "    # Testing data file locations\n",
    "    # TODO: Temp data store, to be replaced with a config spec (.env file?) before deployment\n",
    "        \n",
    "    if exp_cat_dd.value == 'pop':\n",
    "        try:\n",
    "            exp_ras = f\"{DATA_DIR}/EXP/{country}_WPOP20.tif\"\n",
    "        except ValueError:\n",
    "            # TODO: Download from website if layer not available\n",
    "            exp_ras = f\"{DATA_DIR}/EXP/{country}_WPOP20.tif\"\n",
    "    elif exp_cat_dd.value == 'builtup':\n",
    "        exp_ras = f\"{DATA_DIR}/EXP/{country}_WSF19.tif\"\n",
    "    elif exp_cat_dd.value == 'agri':\n",
    "        exp_ras = f\"{DATA_DIR}/EXP/{country}_ESA20_agri.tif\"\n",
    "    else:\n",
    "        ValueError(\"Missing data layer\")\n",
    "\n",
    "    # Hazard data location\n",
    "    hazard_RP_data_loc = f\"{DATA_DIR}/HZD\"\n",
    "\n",
    "    # Run analysis\n",
    "\n",
    "    # Open exposure dataset\n",
    "    exp_data = rxr.open_rasterio(exp_ras)\n",
    "\n",
    "    # Indicate -1 values as representing no data.\n",
    "    exp_data.rio.write_nodata(-1, inplace=True)\n",
    "\n",
    "    # Load ADM based on country code value\n",
    "    try:\n",
    "        adm_dataset = gpd.read_file(os.path.join(DATA_DIR, f\"ADM/{country}_ADM.gpkg\"), layer=f\"{country}_{adm_name}\")\n",
    "    except ValueError:\n",
    "        print(\"Missing ADM layer!\")\n",
    "\n",
    "    adm_data = adm_dataset.loc[adm_dataset.ADM0_CODE == country_code_map[country], :]\n",
    "    \n",
    "    # Get all ADM code/name columns to save with results\n",
    "    adm_cols = adm_data.columns\n",
    "    all_adm_codes = adm_data.columns.str.contains(\"_CODE\")\n",
    "    all_adm_names = adm_data.columns.str.contains(\"_NAME\")\n",
    "    \n",
    "    all_adm_name_tmp = adm_cols[all_adm_names].tolist()\n",
    "    all_adm_code_tmp = adm_cols[all_adm_codes].to_list()\n",
    "\n",
    "    result_df = adm_data.loc[:, all_adm_code_tmp + all_adm_name_tmp + [\"geometry\"]]\n",
    "    \n",
    "    # Prep result structure     \n",
    "    for rp in valid_RPs:\n",
    "        \n",
    "        # Get total exposure for each ADM region\n",
    "        exp_per_ADM = gen_zonal_stats(vectors=adm_data[\"geometry\"], raster=exp_ras, stats=[\"sum\"])\n",
    "        \n",
    "        result_df[f\"{adm_name}_{exp_cat}\"] = [x['sum'] for x in exp_per_ADM]\n",
    "\n",
    "        # Load corresponding hazard dataset\n",
    "        hazard_rst = rxr.open_rasterio(os.path.join(hazard_RP_data_loc, f\"GLB_HS_RP{rp}.tif\"))\n",
    "\n",
    "        # Reproject and clip raster to same bounds as exposure data\n",
    "        hazard_rst = hazard_rst.rio.reproject_match(exp_data)\n",
    "\n",
    "        # Get raw array values for exposure and hazard layer\n",
    "        hazard_arr = hazard_rst[0].values\n",
    "\n",
    "        hazard_arr[hazard_arr < min_haz_threshold] = 0  # Set values below min threshold to nan\n",
    "        hazard_arr[hazard_arr > max_bin_value] = max_haz_threshold  # Cap large values to maximum threshold value\n",
    "        \n",
    "        # Calculate affected exposure in ADM        \n",
    "        # Filter down to valid areas\n",
    "        valid_impact_areas = hazard_rst.values > 0\n",
    "        affected_exp = exp_data.where(valid_impact_areas)  # Get total exposure in affected areas\n",
    "        affected_exp = affected_exp.where(affected_exp > 0)  # Out of the above, get areas that have people\n",
    "        \n",
    "        if save_inter_rst_chk.value:\n",
    "            affected_exp.rio.to_raster(os.path.join(OUTPUT_DIR, f\"{country}_HS_{rp}_{exp_cat}_hazard_affected_.tif\"))\n",
    "            \n",
    "        # Conduct analyses for classes\n",
    "\n",
    "        # Assign bin values to raster data\n",
    "        # Follows: x_{i-1} <= x_{i} < x_{i+1}\n",
    "        bin_idx = np.digitize(hazard_arr, bin_seq)\n",
    "\n",
    "        for bin_x in range(1, num_bins):\n",
    "            impact_class = gen_zonal_stats(vectors=adm_data[\"geometry\"], raster=np.array(bin_idx == bin_x).astype(int) * affected_exp.data[0],\n",
    "                                stats=[\"sum\"], affine=affected_exp.rio.transform(), nodata=np.nan)\n",
    "            result_df[f\"RP{rp}_{exp_cat}_C{bin_x}\"] = [x['sum'] for x in impact_class]\n",
    "        # end\n",
    "\n",
    "        C_cols = result_df.columns.str.contains(f\"RP{rp}_{exp_cat}_C\")\n",
    "        result_df[f\"RP{rp}_{exp_cat}_tot_exposed\"] = result_df.loc[:, C_cols].sum(axis=1)\n",
    "\n",
    "    # End RP loop\n",
    "    \n",
    "    # Round to three decimal places to avoid giving the impression of high precision\n",
    "    result_df = result_df.round(3)\n",
    "\n",
    "    # Write table of total population in each class, in each ADM2\n",
    "    df_cols = result_df.columns\n",
    "\n",
    "    no_geom = result_df.loc[:, df_cols[~df_cols.isin(['geometry'])]].fillna(0)\n",
    "    no_geom.to_csv(os.path.join(OUTPUT_DIR, f\"{country}_HS_{adm_name}_{exp_cat}_class.csv\"), index=False)\n",
    "    result_df.to_file(os.path.join(OUTPUT_DIR, f\"{country}_HS_{adm_name}_{exp_cat}_class.gpkg\"))\n",
    "\n",
    "    with output:\n",
    "        print(\"Finished analysis.\")\n",
    "        rb.disabled = False\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "526bfa3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data option widgets\n",
    "country_dd = widgets.Dropdown(\n",
    "    options=[('Nepal', 'NPL'), ('Pakistan', 'PAK'),('Bangladesh', 'BGD'),('Ghana', 'GHA'),\n",
    "            ('Ethiopia', 'ETH'), ('Burkina Faso', 'BFA'), ('Mali', 'MLI'), ('Niger', 'NER'),\n",
    "            ('Chad', 'TCD'), ('Mauritania', 'MRT')],\n",
    "    value='NPL',\n",
    "    description='Country:',\n",
    "    style={'description_width': 'initial'}\n",
    ")\n",
    "\n",
    "exp_cat_dd = widgets.Dropdown(\n",
    "    options=[(\"Population\", \"pop\"), (\"Built-up\", \"builtup\"), (\"Agriculture\", \"agri\")],\n",
    "    value='pop',\n",
    "    description='Exposure Category:',\n",
    "    style={'description_width': 'initial'}\n",
    ")\n",
    "\n",
    "adm_dd = widgets.Dropdown(\n",
    "    options=['ADM1', 'ADM2', 'ADM3'],\n",
    "    value='ADM2',\n",
    "    description='Administrative Unit Level:',\n",
    "    style={'description_width': 'initial'}\n",
    ")\n",
    "\n",
    "analysis_app_dd = widgets.Dropdown(\n",
    "    options=[\"Classes\"],\n",
    "    value=\"Classes\",\n",
    "    description='Analysis Approach:',\n",
    "    style={'description_width': 'initial'}\n",
    ")\n",
    "\n",
    "# Class value inputs\n",
    "class_edges = OrderedDict({\n",
    "    f'class_{i+1}': widgets.BoundedFloatText(\n",
    "        value=k,\n",
    "        min=18.0,\n",
    "        max=50.0,\n",
    "        step=1.0,\n",
    "        description=f'Class {i+1}:',\n",
    "        tooltip=f'Minimum value of class {i+1}. Value must be less than the next entry.',\n",
    "        disabled=False\n",
    "    ) for (i,k) in enumerate([18, 23, 28, 30])\n",
    "})\n",
    "\n",
    "# User action widgets\n",
    "save_inter_rst_chk = widgets.Checkbox(\n",
    "    value=False,\n",
    "    description='Export Intermediate Rasters',\n",
    "    tooltip='Save rasters generated between each step (saves to nominated output directory)',\n",
    "    disabled=False,\n",
    "    indent=False\n",
    ")\n",
    "\n",
    "# Run button to perform analysis\n",
    "run_button = widgets.Button(\n",
    "    description='Run Analysis',\n",
    "    button_style='', # 'success', 'info', 'warning', 'danger' or ''\n",
    "    tooltip='Click to run analysis with selected options',\n",
    "    # icon='check' # (FontAwesome names without the `fa-` prefix)\n",
    ")\n",
    "\n",
    "reset_display_button = widgets.Button(\n",
    "    description='Reset',\n",
    "    button_style='', # 'success', 'info', 'warning', 'danger' or ''\n",
    "    tooltip='Reset display',\n",
    "    # icon='check' # (FontAwesome names without the `fa-` prefix)\n",
    ")\n",
    "\n",
    "def reset_display(bt):\n",
    "    output.clear_output()\n",
    "    run_button.disabled = False\n",
    "\n",
    "run_button.on_click(run_analysis)\n",
    "reset_display_button.on_click(reset_display)\n",
    "\n",
    "def analysis_change(change):\n",
    "    if change.new == \"Classes\":\n",
    "        disable_class = False\n",
    "        \n",
    "    for w in class_edges.values():\n",
    "        w.disabled = disable_class\n",
    "\n",
    "analysis_app_dd.observe(analysis_change, names='value')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "59a0428c",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "621cfc3843884e16b6d40ebb94bdb79d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Dropdown(description='Country:', options=(('Nepal', 'NPL'), ('Pakistan', 'PAK'), ('Bangladesh', 'BGD'), ('Ghan…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4d503463751049f2917a4ced7c374b82",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Dropdown(description='Exposure Category:', options=(('Population', 'pop'), ('Built-up', 'builtup'), ('Agricult…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2d34fc8e4e064e82b17b7f9c6438dd17",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Dropdown(description='Administrative Unit Level:', index=1, options=('ADM1', 'ADM2', 'ADM3'), style=Descriptio…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "426fd651df934bca95fd74825a88b8a8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Dropdown(description='Analysis Approach:', options=('Classes',), style=DescriptionStyle(description_width='ini…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5db295a62aaf46a395782e820d4a466b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "BoundedFloatText(value=18.0, description='Class 1:', max=50.0, min=18.0, step=1.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5b278b5c4ffa4009acad793f2b1423a9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "BoundedFloatText(value=23.0, description='Class 2:', max=50.0, min=18.0, step=1.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d9b98fe48c3648038791d0ae5ce9e995",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "BoundedFloatText(value=28.0, description='Class 3:', max=50.0, min=18.0, step=1.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d85f432e3bb34beaacc507f112be2cd9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "BoundedFloatText(value=30.0, description='Class 4:', max=50.0, min=18.0, step=1.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2f0bc21898554d7dbee7ed29cb8a21e0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Button(description='Run Analysis', style=ButtonStyle(), tooltip='Click to run analysis with sel…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "01702f991f8f4b129e1e968318e38808",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Button(description='Reset', style=ButtonStyle(), tooltip='Reset display')"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d7fe56dc444c4bc2960d1277cd548572",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(country_dd)\n",
    "display(exp_cat_dd)\n",
    "display(adm_dd)\n",
    "display(analysis_app_dd)\n",
    "[display(w) for w in class_edges.values()]\n",
    "\n",
    "display(HBox([run_button, save_inter_rst_chk]), \n",
    "        reset_display_button)\n",
    "\n",
    "output = widgets.Output()\n",
    "display(output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32340823",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "b6d8b26c4f649b8a27f8996ac75dd87b5987470d4d00d1b2c02ba594e155f4b5"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}