landslide_puertorico.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "522c0574",
   "metadata": {},
   "source": [
    "# Data Component Use Case for Landslide Susceptibility Calculation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f6caae7",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "Landslide susceptibility is the likelihood of a landslide occurring in an area on the basis of local terrain condition to estimate “where” landslides are likely to occur. This Jupyter notebook demonstrates how to use several [CSDMS Data Components](https://csdms.colorado.edu/wiki/DataComponents) to download topography and soil datasets to calculate the landslide susceptibility for a study area in Puerto Rico. \n",
    "\n",
    "In this notebook, it includes the following sections:\n",
    "- [Step 0: Initial Setup](#setup)\n",
    "  \n",
    "  Install API key files and create the input/output folders.\n",
    "  <br>\n",
    "- [Step 1: Download Datasets](#step1) \n",
    "\n",
    "  Download the topography and soil datasets.\n",
    "  <br>\n",
    "- [Step 2: Regrid Datasets](#step2)\n",
    "\n",
    "  Regrid the datasets in the same grid resolution. \n",
    "  <br>\n",
    "- [Step 3: Calculate Susceptibility](#step3)\n",
    "\n",
    "  Use the datasets to calculate the landslide susceptibility.\n",
    "  <br>\n",
    "- [Step 4: Visualize Results](#step4)\n",
    "\n",
    "  Visualize the final results as a short video.\n",
    "  <br>\n",
    "\n",
    "\n",
    "**Suggested Citation**: Gan, T., Campforts, B., Tucker, G. E., Overeem, I. (2023). Data Component Use Case for Landslide Susceptibility Calculation, HydroShare, https://www.hydroshare.org/resource/df5fa2f5d1b74be4bf0a049e1e59889c/\n",
    "\n",
    "**Run this notebook**: Please follow the instructions [here](https://github.com/gantian127/landslide_usecase#data-component-use-case-for-landslide-susceptibility-calculation) to run this notebook on the local PC or the online platform."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48a9f425",
   "metadata": {},
   "source": [
    "<a id='setup'></a>\n",
    "## Initial Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9f8d9bd",
   "metadata": {},
   "source": [
    "### Install API key files\n",
    "For the ERA5 and Topography data components, there is a need to create API key files to download the datasets. The install_api_key( ) function will ask for your [CDS API Key](https://cds.climate.copernicus.eu/api-how-to) and [Open Topography API Key](https://opentopography.org/blog/introducing-api-keys-access-opentopography-global-datasets) to create API key files. Please make sure you have already obtained those API Keys before you run this helper function. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96ddd865",
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils import install_api_key\n",
    "install_api_key()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90f2bfe5",
   "metadata": {},
   "source": [
    "### Create folders\n",
    "We will first import all the python packages and then create three folders for this notebook:\n",
    "- **configuration file folder**: this folder includes several configuration files which will be used by the data components. In this example, we have prepared these configuration files ('dem_config.yaml' and 'era5_config.yaml') and put them in this folder. \n",
    "- **cache folder**: this folder stores the downloaded data files.\n",
    "- **results folder**: this folder stores the final results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1f07eb7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import packages \n",
    "import os\n",
    "import warnings\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray\n",
    "import xesmf as xe\n",
    "import rioxarray \n",
    "import cftime\n",
    "from datetime import datetime\n",
    "from tqdm import trange\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import colors\n",
    "import imageio.v2 as imageio\n",
    "from IPython.display import Video\n",
    "\n",
    "from landlab import RasterModelGrid\n",
    "from pymt.models import Topography, Era5\n",
    "\n",
    "from utils import regrid_data, cal_subsurface_flow_depth, cal_safety_factor\n",
    "\n",
    "warnings.simplefilter(action='ignore', category=FutureWarning)\n",
    "plt.rcParams.update({'font.size': 14})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82f450b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create folders\n",
    "study_area = 'puerto_rico'\n",
    "\n",
    "config_dir = os.path.join(os.getcwd(), 'config_files_{}'.format(study_area))\n",
    "results_dir = os.path.join(os.getcwd(), 'results_{}'.format(study_area)) \n",
    "cache_dir = os.path.join(os.getcwd(),'cache_{}'.format(study_area))\n",
    "\n",
    "\n",
    "for folder in [config_dir, results_dir, cache_dir]:\n",
    "    if not os.path.isdir(folder):\n",
    "        os.mkdir(folder)\n",
    "        print(folder)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16eca92d",
   "metadata": {},
   "source": [
    "<a id='step1'></a>\n",
    "## Step 1 Download Datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f60dfec",
   "metadata": {},
   "source": [
    "### Background \n",
    "\n",
    "Landslides are frequent hazards in Puerto Rico which are mainly caused by the steep terrain and heavy rainfall from hurricanes and other tropical weather systems. For example, Hurricane Maria hit the island of Puerto Rico on September 20th, 2017 and triggered more than 40,000 landslides in Puerto Rico (see details [here](https://www.usgs.gov/supplemental-appropriations-for-disaster-recovery-activities/landslides-triggered-hurricane-maria)). The map below shows concentration and distribution of landslides generated by rainfall associated with Hurricane Maria in Puerto Rico (map source from USGS).\n",
    "\n",
    "In this example, we will calculate the hourly landslide susceptibility for the area that has high concentration of landslide of Puerto Rico during Hurricane Maria. We will prepare the following datasets to calculate the landslide susceptibility. Details for how to calculate the landslide susceptibility is described in [Step 3](#step3)\n",
    "- OpenTopography DEM \n",
    "- ERA5 volumetric soil water and precipitation \n",
    "- Soil depth  \n",
    "- Slope angle\n",
    "\n",
    "<img src=\"https://github.com/gantian127/landslide_usecase/blob/master/Landslide_map.jpeg?raw=true\" width=\"900\">"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "573e5e1d",
   "metadata": {},
   "source": [
    "### OpenTopography DEM \n",
    "\n",
    "We will use the [Topography data component](https://csdms.colorado.edu/wiki/Model:Topography_Data_Component) to download the Digital Elevation Model (DEM) data with 90m resolution. The figure below shows the bounding box of the study area. The 'dem_config.yaml' file includes the parameter settings of this data component. The following cells demonstrate how to use the configuration file to initialize a data component and how to use the variable and grid related methods of this data component to get the metadata as well as the data values.\n",
    "\n",
    "\n",
    "<img src=\"https://github.com/gantian127/landslide_usecase/blob/master/study_area.png?raw=true\" width=\"850\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a264880b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialize Topography data component\n",
    "dem = Topography()\n",
    "dem.initialize(os.path.join(config_dir, 'dem_config.yaml'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "239ce76e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get DEM variable info\n",
    "var_name = dem.output_var_names[0]\n",
    "var_unit = dem.var_units(var_name)\n",
    "var_location = dem.var_location(var_name)\n",
    "var_type = dem.var_type(var_name)\n",
    "var_grid = dem.var_grid(var_name)\n",
    "var_itemsize = dem.var_itemsize(var_name)\n",
    "var_nbytes = dem.var_nbytes(var_name)\n",
    "print('variable_name: {} \\nvar_unit: {} \\nvar_location: {} \\nvar_type: {} \\nvar_grid: {} \\nvar_itemsize: {}' \n",
    "            '\\nvar_nbytes: {} \\n'. format(var_name, var_unit, var_location, var_type, var_grid, var_itemsize, var_nbytes))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "985d5bea",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get DEM grid info \n",
    "dem_grid_ndim = dem.grid_ndim(var_grid) \n",
    "dem_grid_type = dem.grid_type(var_grid)\n",
    "dem_grid_shape = dem.grid_shape(var_grid)\n",
    "dem_grid_spacing = dem.grid_spacing(var_grid)\n",
    "dem_grid_origin = dem.grid_origin(var_grid)\n",
    "\n",
    "print('grid_ndim: {} \\ngrid_type: {} \\ngrid_shape: {} \\ngrid_spacing: {} \\ngrid_origin: {}'.format(\n",
    "    dem_grid_ndim, dem_grid_type, dem_grid_shape, dem_grid_spacing, dem_grid_origin))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a9816f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get DEM variable data\n",
    "dem_data = dem.get_value(var_name)\n",
    "dem_data_2D = dem_data.reshape(dem_grid_shape)\n",
    "\n",
    "# get X, Y extent for plot\n",
    "min_y, min_x = dem_grid_origin\n",
    "max_y = min_y + dem_grid_spacing[0]*(dem_grid_shape[0]-1)\n",
    "max_x = min_x + dem_grid_spacing[1]*(dem_grid_shape[1]-1)\n",
    "dy = dem_grid_spacing[0]/2\n",
    "dx = dem_grid_spacing[1]/2\n",
    "dem_extent = [min_x - dx, max_x + dx, min_y - dy, max_y + dy]\n",
    "\n",
    "# plot DEM data\n",
    "fig, ax = plt.subplots(1,1,figsize=(10,5))\n",
    "im = ax.imshow(dem_data_2D, extent=dem_extent)\n",
    "ax.title.set_text('Topography Data')\n",
    "ax.set_xlabel('longitude [degree_east]')\n",
    "ax.set_ylabel('latitude [degree_north]')\n",
    "fig.colorbar(im,label='elevation(m)')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d34b5f63",
   "metadata": {},
   "source": [
    "### ERA5 Volumetric Soil Water & Precipitation \n",
    "We will use the [ERA5 data component](https://csdms.colorado.edu/wiki/Model:ERA5_Data_Component) to download the hourly volumetric soil water data and the precipitation data of the study area with 0.25 degrees (27-28km) resolution. The volumetric soil water data will be used for calculating the susceptibility, while the precipitation data is mainly used for results visualization. \n",
    "\n",
    "The 'era5_config.yaml' file includes the parameter settings of this data component. The following cells demonstrate how to use the configuration file to initialize an ERA5 data component and how to use the variable, grid and time related methods to get the metadata as well as the data values. You'll notice that although the ERA5 and Topography data components download the datasets from different sources, they are using the same methods to get information from the datasets. Please note that sometimes the request for ERA5 data may be queued which may take a while (>10min) to get the data downloaded."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c730ff9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialize ERA5 data component\n",
    "era5 = Era5()\n",
    "era5.initialize(os.path.join(config_dir,'era5_config.yaml'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cece7801",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get ERA5 variable info\n",
    "for var_name in era5.output_var_names:\n",
    "    var_unit = era5.var_units(var_name)\n",
    "    var_location = era5.var_location(var_name)\n",
    "    var_type = era5.var_type(var_name)\n",
    "    var_grid = era5.var_grid(var_name)\n",
    "    var_itemsize = era5.var_itemsize(var_name)\n",
    "    var_nbytes = era5.var_nbytes(var_name)\n",
    "    print('variable_name: {} \\nvar_unit: {} \\nvar_location: {} \\nvar_type: {} \\nvar_grid: {} \\nvar_itemsize: {}' \n",
    "            '\\nvar_nbytes: {} \\n'. format(var_name, var_unit, var_location, var_type, var_grid, var_itemsize, var_nbytes))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1bc1e2e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get ERA5 grid info\n",
    "era5_grid_ndim = era5.grid_ndim(var_grid) \n",
    "era5_grid_type = era5.grid_type(var_grid)\n",
    "era5_grid_shape = era5.grid_shape(var_grid)\n",
    "era5_grid_spacing = era5.grid_spacing(var_grid)\n",
    "era5_grid_origin = era5.grid_origin(var_grid)\n",
    "\n",
    "print('grid_ndim: {} \\ngrid_type: {} \\ngrid_shape: {} \\ngrid_spacing: {} \\ngrid_origin: {}'.format(\n",
    "    era5_grid_ndim, era5_grid_type, era5_grid_shape, era5_grid_spacing, era5_grid_origin))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "894f48bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get ERA5 time info\n",
    "era5_start_time = era5.start_time\n",
    "era5_end_time = era5.end_time\n",
    "era5_time_step = era5.time_step\n",
    "era5_time_unit = era5.time_units\n",
    "era5_time_steps = int((era5_end_time - era5_start_time)/era5_time_step) + 1\n",
    "\n",
    "print('start_time:{} \\nend_time:{} \\ntime_step:{} \\ntime_unit:{} \\ntime_steps:{}'.format(\n",
    "    era5_start_time, era5_end_time, era5_time_step, era5_time_unit, era5_time_steps))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f5583aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get ERA5 variables data and plot (at the first time step)\n",
    "fig = plt.figure(figsize=(18,16)) \n",
    "nrows, ncols = 3, 2\n",
    "i = 1\n",
    "\n",
    "for var_name in era5.output_var_names:\n",
    "    ax = fig.add_subplot(nrows, ncols, i)\n",
    "    var_unit = era5.var_units(var_name)\n",
    "    \n",
    "    # get variable data    \n",
    "    era5_data = era5.get_value(var_name)\n",
    "    era5_data_2D = era5_data.reshape(era5_grid_shape)\n",
    "    \n",
    "    # get X, Y extent for plot\n",
    "    min_y, min_x = era5_grid_origin\n",
    "    max_y = min_y + era5_grid_spacing[0]*(era5_grid_shape[0]-1)\n",
    "    max_x = min_x + era5_grid_spacing[1]*(era5_grid_shape[1]-1)\n",
    "    dy = era5_grid_spacing[0]/2\n",
    "    dx = era5_grid_spacing[1]/2\n",
    "    era5_extent = [min_x - dx, max_x + dx, min_y - dy, max_y + dy]\n",
    "\n",
    "    # plot data\n",
    "    im = ax.imshow(era5_data_2D, extent=era5_extent, cmap='Blues')\n",
    "    ax.title.set_text('{} ({})'.format(var_name,var_unit ))\n",
    "    ax.set_xlabel('longitude [degree_east]')\n",
    "    ax.set_ylabel('latitude [degree_north]')\n",
    "    cbar = plt.colorbar(im, ax=ax)\n",
    "    \n",
    "    i += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b2b0ded",
   "metadata": {},
   "source": [
    "### Soil Depth data\n",
    "\n",
    "Since the data component for soil depth data is not available yet, we will use the xarray and rioxarray to subset and download the soil depth data with 250m resolution from [SoilGrids](https://www.isric.org/explore/soilgrids) system. In this dataset, the maximum soil depth value is 200cm. Grid with values larger than 200cm represents open water area. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa492bb7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# download soil depth data\n",
    "soil_raster = rioxarray.open_rasterio(\"https://files.isric.org/soilgrids/former/2017-03-10/data/BDRICM_M_250m_ll.tif\")\n",
    "soil_depth_data = soil_raster.rio.clip_box(\n",
    "    minx=-66.8,\n",
    "    miny=18.0,\n",
    "    maxx=-66.2,\n",
    "    maxy=18.45,\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e196652",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot soil depth data\n",
    "soil_depth_data.plot(figsize=(10,5),cbar_kwargs={'label': 'depth(m)'})\n",
    "soil_depth_data.rio.to_raster(os.path.join(cache_dir,'soil_depth.tif'))\n",
    "plt.title('Soil Depth')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e12ca486",
   "metadata": {},
   "source": [
    "### Slope angle\n",
    "\n",
    "Slope angle is one of the input for calculating the landslide susceptibility. In this example, we will use the Topography data and the [RasterModelGrid](https://landlab.readthedocs.io/en/master/reference/grid/raster.html) from [Landlab](https://landlab.readthedocs.io/en/master/index.html) to calculate the slope angle for the study area. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b0163f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate slope angle using Topography data\n",
    "model_grid = RasterModelGrid(dem_data_2D.shape,xy_spacing=(90,90))\n",
    "slope = model_grid.calc_slope_at_node(elevs=dem_data) # slope in radians, 1D array\n",
    "slope_angle = slope.reshape(dem_data_2D.shape) # reshape as 2D array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8cf0fcf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot slope angle\n",
    "fig, ax = plt.subplots(figsize=(10,5))\n",
    "im=ax.imshow(slope_angle, extent=dem_extent)\n",
    "cbar = fig.colorbar(im, label= 'radians')\n",
    "ax.title.set_text('Slope Angle')\n",
    "ax.set_xlabel('longitude [degree_east]')\n",
    "ax.set_ylabel('latitude [degree_north]')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6222e476",
   "metadata": {},
   "source": [
    "<a id='step2'></a>\n",
    "## Step 2 Regrid datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d310a739",
   "metadata": {},
   "source": [
    "Before calculating the landslide susceptibility, we need to regrid all the datasets in the same grid resolution. We will take the Topography dataset as the template to regrid the soil depth and the ERA5 datasets. This will make these datasets to have the same resolution as the Topography dataset. Since all the original datasets are in the WGS84 geographic coordinate system, there is no need for data reprojection. Because ERA5 datasets are multidimensional space time data, we will regrid ERA5 datasets in a for loop in [Step3](#step3). In this example, we will use the regrid_data( ) function for regridding purpose. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0365ba3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define destination grid coordinate using Topography data\n",
    "dem_y = np.flip(np.arange(dem_grid_shape[0])*dem_grid_spacing[0] + dem_grid_origin[0])\n",
    "dem_x = np.arange(dem_grid_shape[1])*dem_grid_spacing[1] + dem_grid_origin[1]\n",
    "dest_coor = {'lon': dem_x, \n",
    "             'lat': dem_y}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "498b9f1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# regrid soil depth data\n",
    "soil_depth_coor = {'lon': soil_depth_data.x.values,\n",
    "                   'lat': soil_depth_data.y.values}\n",
    "\n",
    "soil_depth = regrid_data(soil_depth_data.values[0], soil_depth_coor, dest_coor) \n",
    "soil_depth = soil_depth/100 # units conversion as m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19e2f66d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot regridded soil depth data\n",
    "fig, ax = plt.subplots(1,1,figsize=(10,5))\n",
    "im = ax.imshow(soil_depth, extent=dem_extent)\n",
    "ax.title.set_text('Regridded Soil Depth Data')\n",
    "ax.set_xlabel('longitude [degree_east]')\n",
    "ax.set_ylabel('latitude [degree_north]')\n",
    "cbar = plt.colorbar(im, label='depth(m)')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96e9588c",
   "metadata": {},
   "source": [
    "<a id='step3'></a>\n",
    "## Step 3  Calculate Susceptibility"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26df085c",
   "metadata": {},
   "source": [
    "This section will loop through each time step to calculate the hourly landslide susceptibility, which mainly includes the following tasks:\n",
    "- Regrid ERA5 datasets\n",
    "- Calculate subsurface flow depth\n",
    "- Calculate factor of safety\n",
    "- Calculate susceptibility\n",
    "- Plot results\n",
    "\n",
    "\n",
    "**Regrid ERA5 datasets**\n",
    "\n",
    "The ERA5 precipitation and volumetric soil water datasets will be regridded in a higher resolution by using the Topography dataset as the template. The regrid_data( ) function will be used which is similar as the code shown in [Step2](#step2). As mentioned before, the volumetric soil water data is used for calculation and the precipitation data is mainly used to visualize the rainfall change during the time period. \n",
    "\n",
    "**Subsurface flow depth** \n",
    "\n",
    "The subsurface flow depth is calculated using the soil depth data and the ERA5 volumetric soil water datasets. There are four layers of the ERA5 volumetric soil water data (Layer 1: 0 - 7cm, Layer 2: 7 - 28cm, Layer 3: 28 - 100cm, Layer 4: 100 - 289cm). We first calculate the product of the soil depth at each layer and its corresponding volumetric soil water data. Then this product is divided by the porosity, which represents the water depth at each soil layer. Then we add up the water depth values of all four layers, which will be an approximation of the subsurface flow depth. The cal_subsurface_flow_depth( ) function is used for calculation.\n",
    "\n",
    "**Factor of safety**\n",
    "\n",
    "In geological engineering, it is common to take the ratio of the resisting stresses to driving stresses. This ratio is called the factor of safety (FS). When FS is larger than 1, the slope should be stable, while if it is below 1, the driving stress exceeds the resistance and the slope is likely to fail. FS can be calculated with the following function, and cal_safety_factor( ) is implemented based on this function.\n",
    "\n",
    "$$\n",
    "FS = \\frac{(C_r + C_s)/h_s\\rho_sg}{\\sin\\theta} + \\frac{\\cos\\theta \\tan\\phi (1-\\frac{h_w}{h_s}\\rho_w / \\rho_s)}{\\sin\\theta}\n",
    "$$\n",
    "\n",
    "where, \n",
    "- Cr: root cohesion (Pa kg/ms^2)\n",
    "- Cs: soil cohesion (Pa kg/ms^2)\n",
    "- hs: soil depth (m)\n",
    "- hw: subsurface flow depth (m)\n",
    "- ρs: soil density (kg/m^3)\n",
    "- ρw: water density (kg/m^3)\n",
    "- g: gravity acceleration (m/s^2)\n",
    "- θ: slope angle \n",
    "- φ: soil internal friction angle \n",
    "\n",
    "**Susceptibility**\n",
    "\n",
    "Susceptibility is the inverse of FS. When the susceptibility is larger than 1, it means that the slope of the area is not stable and susceptible to landslide.  \n",
    "\n",
    "$$\n",
    "\\text{Susceptibility} = \\frac{1}{FS}\n",
    "$$\n",
    "\n",
    "\n",
    "**Plot Results**\n",
    "\n",
    "There are four subplots created at each time step. The top left plot shows the landslide susceptibility at first time step and then it shows the difference of the landslide susceptibility between the current and the first time step. Similarly, the lower left plot shows the subsurface flow depth at the first time step and then the difference of the subsurface flow depth. The two plots on the right show the time series of mean total precipitation and the mean volumetric soil water content for the study area respectively. Please note that these two plots are not shown at the first time step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5efc9d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set time steps for the loop\n",
    "time_steps = 48\n",
    "\n",
    "# define mask for non-data area\n",
    "mask = (slope_angle==0)&(soil_depth>2.0)\n",
    "\n",
    "# define a dataframe for ERA5 time series data\n",
    "df_columns = ['Time', 'Prcp', 'VSW_L1', 'VSW_L2', 'VSW_L3','VSW_L4']\n",
    "era5_ts_df = pd.DataFrame(columns = df_columns)\n",
    "\n",
    "# define ERA5 input coordinate and an empty array for regridding results\n",
    "era5_y = np.flip(np.arange(era5_grid_shape[0])*era5_grid_spacing[0] + era5_grid_origin[0])\n",
    "era5_x = np.arange(era5_grid_shape[1])*era5_grid_spacing[1] + era5_grid_origin[1]\n",
    "era5_coor = {'lon': era5_x, \n",
    "             'lat': era5_y}\n",
    "\n",
    "soil_water_layer = np.empty([4,*dem_grid_shape], dtype=era5.var_type('Volumetric soil water layer 1'))\n",
    "\n",
    "\n",
    "# calculation and plot\n",
    "for time_step in trange(0, time_steps): \n",
    "\n",
    "    # get time value\n",
    "    time = cftime.num2date(era5.time, era5_time_unit, only_use_cftime_datetimes=False, only_use_python_datetimes=True)\n",
    "    \n",
    "    # regrid ERA5 precipitation \n",
    "    era5_prec = era5.get_value('Total precipitation')\n",
    "    era5_prec_2D = era5_prec.reshape(era5_grid_shape)\n",
    "    era5_prec_data = regrid_data(era5_prec_2D, era5_coor, dest_coor)\n",
    "    \n",
    "    \n",
    "    # regrid ERA5 volumetric soil water (4 layers)\n",
    "    era5_soil_ts = []\n",
    "    for i in range(0,4):\n",
    "        var_name = 'Volumetric soil water layer {}'.format(i+1)\n",
    "\n",
    "        # get original data\n",
    "        era5_soil_water = era5.get_value(var_name)\n",
    "        era5_soil_water_2D = era5_soil_water.reshape(era5_grid_shape)\n",
    "        \n",
    "        # regrid\n",
    "        soil_water_layer[i] = regrid_data(era5_soil_water_2D, era5_coor, dest_coor)\n",
    "        era5_soil_ts.append(np.mean(soil_water_layer[i]))\n",
    "        \n",
    "    # store ERA5 time series data\n",
    "    era5_ts_df.loc[len(era5_ts_df)] = [time, np.mean(era5_prec_data)] + era5_soil_ts\n",
    "       \n",
    "    # calculate subsurface flow depth\n",
    "    subsurface_flow_depth = cal_subsurface_flow_depth(soil_depth, soil_water_layer)\n",
    "    \n",
    "    # calculate FS\n",
    "    safety_factor = cal_safety_factor(slope_angle, subsurface_flow_depth, soil_depth, \n",
    "                                      root_cohesion=5000, soil_cohesion=5000, soil_bulk_density=1300,\n",
    "                                      soil_internal_friction_angle=35)\n",
    "    \n",
    "    # calculate susceptibility\n",
    "    susceptibility = 1.0 / safety_factor\n",
    "    susceptibility = np.where(~mask, susceptibility, np.nan)\n",
    "    \n",
    "    # plot susceptibility\n",
    "    fig = plt.figure(figsize=(16,10))\n",
    "    \n",
    "    fig.suptitle(\"Results at {}\".format(time))\n",
    "    plt.rcParams[\"figure.autolayout\"] = True\n",
    "    \n",
    "    nrows, ncols = 2,2\n",
    "    ax_sus = fig.add_subplot(nrows, ncols, 2)\n",
    "    ax_subflow = fig.add_subplot(nrows, ncols, 4)\n",
    "    ax_prcp = fig.add_subplot(nrows, ncols, 1)\n",
    "    ax_soil = fig.add_subplot(nrows, ncols, 3)\n",
    "    \n",
    "    \n",
    "    # susceptibility plot \n",
    "    if time_step==0:\n",
    "        im_data = susceptibility\n",
    "        reference_sus = np.copy(susceptibility)\n",
    "        \n",
    "        divnorm = None\n",
    "        ax_sus.title.set_text('Susceptibility of landslide')\n",
    "    else:\n",
    "        im_data = susceptibility - reference_sus \n",
    "        \n",
    "        divnorm = colors.TwoSlopeNorm(vcenter=0.0)       \n",
    "        ax_sus.title.set_text('Difference of susceptibility of landslide')\n",
    "    \n",
    "    ax_sus.set_xlabel('longitude [degree_east]')\n",
    "    ax_sus.set_ylabel('latitude [degree_north]')\n",
    "    \n",
    "    im_sus = ax_sus.imshow(im_data, cmap='BrBG_r', norm=divnorm, extent=dem_extent)\n",
    "    fig.colorbar(im_sus, ax=ax_sus)\n",
    "    \n",
    "    \n",
    "    # subsurface flow depth plot\n",
    "    if time_step==0:\n",
    "        im_data = subsurface_flow_depth\n",
    "        ax_subflow.title.set_text('Subsurface flow depth (m)')\n",
    "        reference_subflow = np.copy(subsurface_flow_depth) \n",
    "    else:\n",
    "        im_data = subsurface_flow_depth - reference_subflow\n",
    "        ax_subflow.title.set_text('Difference of subsurface flow depth (m)')\n",
    "    \n",
    "    im_data = np.where(~mask, im_data, np.nan)\n",
    "    ax_subflow.set_xlabel('longitude [degree_east]')\n",
    "    ax_subflow.set_ylabel('latitude [degree_north]')\n",
    "    \n",
    "    divnorm=colors.TwoSlopeNorm(vcenter=0.0) \n",
    "    im_subflow = ax_subflow.imshow(im_data, cmap='RdBu', norm=divnorm,  extent=dem_extent)\n",
    "    fig.colorbar(im_subflow, ax=ax_subflow)\n",
    "    \n",
    "  \n",
    "    if time_step==0:\n",
    "        ax_prcp.axis('off')\n",
    "        ax_soil.axis('off')\n",
    "    \n",
    "    else:\n",
    "        # mean precipitation plot\n",
    "        era5_ts_df.plot('Time','Prcp', ax=ax_prcp).legend(loc='upper right')\n",
    "        ax_prcp.set_ylabel('total precipitation (m/hr)')\n",
    "        ax_prcp.set_ylim([0.0, 0.02])\n",
    "        ax_prcp.title.set_text('Mean total precipitation')\n",
    "\n",
    "        # mean volumetric soil water plot\n",
    "        era5_ts_df.plot('Time',df_columns[2:], ax=ax_soil).legend(loc='upper right')\n",
    "        ax_soil.set_ylabel('volumetric soil water (m^3/m^3)')\n",
    "        ax_soil.set_ylim([0.24, 0.40])\n",
    "        ax_soil.title.set_text('Mean volumetric soil water content')\n",
    "    \n",
    "    # save plot    \n",
    "    plt.close(fig)\n",
    "    fig.savefig(os.path.join(results_dir, 'sus_{}.png'.format(time_step)))\n",
    "    \n",
    "    # update time step\n",
    "    era5.update()\n",
    "\n",
    "print('calculation is done')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "763b528e",
   "metadata": {},
   "source": [
    "<a id='step4'></a>\n",
    "## Step 4 Visualize Results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2425217",
   "metadata": {},
   "source": [
    "Run the cells below and it will show a short video. In the mean total precipitation plot, you will find that there was heavy rain during Sept 20th and the volumetric soil water content at 1st and 2nd layer responds to this rainfall event quickly to reach its capacity. There is a bit time lag for the 3rd layer to increase the volumetric soil water content and there is little change for 4th layer when compared with other layers. You will also notice that the increase of the landslide susceptibility has a positive relationship with the increase of the subsurface flow depth and the slope angle of the study area. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "947a1385",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Make a short video\n",
    "img_files = [os.path.join(results_dir, file) for file in os.listdir(results_dir) if '.png' in file]\n",
    "img_files.sort(key=lambda x: os.path.getmtime(x))\n",
    "\n",
    "with imageio.get_writer(os.path.join(results_dir,'landslide.mp4'), mode='I', fps=1) as writer:\n",
    "    for f in img_files:\n",
    "        im = imageio.imread(os.path.join(results_dir, f))\n",
    "        writer.append_data(im)\n",
    "\n",
    "writer.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83cbbab0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display the video\n",
    "Video(\"./results_puerto_rico/landslide.mp4\", embed=True, width=900, height=700)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2af04b5b",
   "metadata": {},
   "source": [
    "## References\n",
    "- Strauch, R., Istanbulluoglu, E., Nudurupati, S. S., Bandaragoda, C., Gasparini, N. M., and Tucker, G. E. (2018), A hydroclimatological approach to predicting regional landslide probability using Landlab, Earth Surf. Dynam., 6, 49–75, https://doi.org/10.5194/esurf-6-49-2018\n",
    "- Montgomery, D. R., and Dietrich, W. E. (1994), A physically based model for the topographic control on shallow landsliding, Water Resour. Res., 30( 4), 1153– 1171, https://doi.org/10.1029/93WR02979."
   ]
  }
 ],
 "metadata": {
  "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.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}