{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Content:\n",
    "\n",
    "* Installation\n",
    "* ***Preparing model input parameters and dataset***\n",
    "* Runing DRYP model\n",
    "* Post processing model outputs\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preprocessing datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This section will help you to create dataset required as model parameters for runing DRYP. This section uses\n",
    "simple raster operations to store and save variables (scripts can be accessed at https://github.com/AndresQuichimbo/DRYPv2.0.1.git.).\n",
    "More complex operations can be performed using any GIS software (e.g ArcMap, QGis) o python libraries such as rasterio."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this section we will use ready-to-use datasets created for running the regional model. Raw data can be found in different\n",
    "repositories such as Hydrosheds, SoilGrids, etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following steps will be carried out in order to produce the set of model paramter datasets for a subbasin of the HAD region."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Delineate a new subbasin from the regional model domain\n",
    "* Clip extent of model parameters to reduce simulation times (reduce model domain)\n",
    "* Create a flow direction map compatible with DRYP-landlab formart\n",
    "* Create a river network\n",
    "* Create a constant head boundary condition from water table depth\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The preprocessing (also postprocessing) module is under developement therefore the list of functions available are still limited. All available functions can be found in the following link: https://github.com/AndresQuichimbo/DRYPv2.0.1.git."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import libraries from local repository\n",
    "import sys\n",
    "sys.path.append('C:/Users/Edisson/Documents/GitHub/DRYPv2.0.1') # CHANGE THIS TO YOUR MODEL PATH\n",
    "\n",
    "\n",
    "import cuwalid.dryp.components.DRYP_watershed as ppbasin\n",
    "import cuwalid.tools.DRYP_pptools as pptools\n",
    "import cuwalid.tools.DRYP_rrtools as rrtools"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**NOTE**: Before we start, please download and unzip all files required for this training by cliking in the following link: \n",
    "https://drive.google.com/drive/folders/1Lf9fdOjjZO87sU1PCsUbKTNjt5IXr9DD?usp=drive_link"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The file is organized as bellow:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Local_directory\\n|__training\\n      |__basin\\n      |   |__datasets\\n      |   |    |__csv\\n      |   |    |__shp\\n      |   |__model\\n      |   |    |__inputs\\n      |   |__outputs\\n      |__regional\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**WARNING**: Once the training files are downloaded change the paths below to access them during this training material."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "training_general_path = \"D:/HAD/training/\"\n",
    "regional_path = \"D:/HAD/training/regional/\"\n",
    "basin_path = \"D:/HAD/training/basin/\"\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1. Delineate a watershed (getting a basin mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For basin delineation, you will need to provide a Digital Elevation Map (DEM), the flow direction, and the basin outlet as csv file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, specify path of parameter files and create maps for each one of them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fname_surface = regional_path + \"/model/inputs/HAD_DEM_utm_mm.asc\"\n",
    "fname_flowdir = regional_path + \"/model/inputs/HAD_flowdir_land_utm.asc\"\n",
    "fname_outlet = regional_path + \"/model/inputs/HAD_tana_dryp_station_utm.csv\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**TASK**: Visualize the raster datasets, use the following function to quickly create a map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import rasterio\n",
    "# this function can be used for ploting raster maps, modify it as your conveniency.\n",
    "def plot_raster_file(fname, ax=None, vmin=-20.0, vmax=20.0):\n",
    "    # create plot\n",
    "    if ax is None:\n",
    "        fig, ax = plt.subplots()\n",
    "    cmap = plt.cm.get_cmap('coolwarm_r', 12)\n",
    "    data = rasterio.open(fname).read(1)\n",
    "    im = ax.imshow(data,# origin=\"lower\",#cmap=cmap, \n",
    "                   #vmin=vmin, vmax=vmax,\n",
    "                   )#extent=bounds)\t\n",
    "    \n",
    "    ax.axis('off')\n",
    "    plt.colorbar(im)\n",
    "    return im"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_raster_file(fname_surface)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_raster_file(regional_path + \"/model/inputs/HAD_flowdir_D8.asc\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_raster_file(fname_flowdir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Delineate a basin by providing a surface elevation and a the location of the basin outlet. A \n",
    "flow direction file can also be provided."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fname_out = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ppbasin.get_watershed_mask(fname_surface, fname_outlet, fname_flowDir=fname_flowdir, fname_out=fname_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_raster_file(regional_path + \"/model/inputs/HAD_DEM_utm_mm_basin.asc\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read the file created in the previous step. This file is stored in the same directory of the surface file and it will have the same name as the surface (input file) but added \"_basin\" at the end (If a name of the outout filename is provided use that name in the following step).\n",
    "\n",
    "Specify the name of the output file in the previous example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fname_mask = regional_path + \"/model/inputs/HAD_DEM_utm_mm_basin.asc\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ppbasin.get_watershed_area(fname_surface, fname_outlet, fname_flowDir=fname_flowdir, fname_mask=fname_mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**TASK** Read the csv file created in the previous state to check the catchment area"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**TASK**: Visualize the created datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2. Clip basin extent by mask"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are planing to run a subbasin of the regional model, so we are going to create input files for the new model.\n",
    "We create this files by clipling the datasets of the regional model by the extent of the subbasin.\n",
    "\n",
    "This following step will clip the raster files by identifiying the extend of the basin (positive values) provided as mask."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2.1. Clip flow accumulation file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fname = regional_path + \"model/inputs/HAD_DEM_utm_mm_flowaccum.asc\"\n",
    "fmask = regional_path + \"model/inputs/HAD_DEM_utm_mm_basin.asc\"\n",
    "fname_output = training_general_path + \"TA_HAD_DEM_utm_mm_flowaccum.asc\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rrtools.clip_raster_by_mask(fname, fmask, fname_output)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2.2. Clip flow direction map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fname_flowdirD8 = regional_path + \"model/input/HAD_flowdir_D8.asc\"\n",
    "fmask = regional_path + \"model/input/HAD_DEM_utm_mm_basin.asc\"\n",
    "fname_bflowdir = basin_path + \"model/input/TA_HAD_flowdir_D8_utm.asc\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rrtools.clip_raster_by_mask(fname_flowdirD8, fmask, fname_bflowdir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2.3. Create model parameter datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**TASK**: Create all parameters files required for the Tana basin, make sure that all files have the same grid size. A list of files required for running the Tana basin can be found in the parameter input file."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**TASK**: Visualize the created datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3. Create a flow direction map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "DRYP uses a flow direction map (optional, but recomended!) based on LandLab based notation \n",
    "(https://landlab.readthedocs.io/en/latest/reference/components/flow_director.html), so in order to use\n",
    "flow direction maps from other source that are generally in D8 format (e.g. Hydrosheds) it has to be transformed\n",
    "into LandLab format. The following function will help with this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fflowdird8 = basin_path + \"model/input/TA_HAD_flowdir_D8_utm.asc\"\n",
    "fname_bflowdirLandlab = basin_path + \"model/input/TA_HAD_flowdir_landlab_utm.asc\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rrtools.create_raster_flowdirection_dryp(fflowdird8, fname_bflowdirLandlab)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**TASK**: Visualize the created datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4. Creating river network"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "River networks can be created form the flow accumulation map. Here, we use the flow accumalation map\n",
    "created in the previous step to build the river network. The river network can be created based on a\n",
    "minimum number of cells or the minimun area."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fname_flowaccum = basin_path + \"model/input/TA_HAD_DEM_utm_mm_flowaccum.asc\"\n",
    "threshold = 10\n",
    "fname_output = basin_path + \"model/input/TA_HAD_river_network_utm.asc\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rrtools.create_raster_river_network(fname_flowaccum, threshold, fname_output, cell_area=False, fill_value=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**TASK**: Change the threshold to see how the river network varies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**TASK**: Visualize the created datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5. Create constant head boundary condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fname_wte = regional_path + \"model/inputs/HAD_long_term_wte_ini.asc\"\n",
    "fname_outlet = regional_path + \"model/inputs/HAD_tana_dryp_station_utm.csv\"\n",
    "fname_output = basin_path + \"model/input/TA_HAD_CHB_utm.asc\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rrtools.create_raster_bc_at_point(fname_wte, fname_outlet, fname_output)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**TASK**: Modify the list of points in the csv file create a new file for head boundary conditions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**6. Clip all raster parameters files for running the model**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import glob\n",
    "\n",
    "files = glob.glob(regional_path+\"/model/inputs/*.asc\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_path = \"D:/HAD_basins/Kenya/input/KE_\"\n",
    "fnames = []\n",
    "for ifile in files:\n",
    "    #fnames.append(ifile.split(\"\\\\\")[-1])\n",
    "    new_name = new_path + ifile.split(\"\\\\\")[-1]\n",
    "    fnames.append(new_name)\n",
    "    print(new_name)\n",
    "    #new_name = new_path + fnames[0]\n",
    "    #process clip raster\n",
    "    #rrtools.clip_raster_by_mask(ifile, fmask, new_name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**7. Check that all raster datasets are correct**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#### fname_base needs to be updated\n",
    "fname_base = basin_path+\"Kenya/input/KE_HAD_DEM_utm_mm.asc\"\n",
    "for ifname in fnames:\n",
    "    rrtools.check_raster_alignaments(fname_base, ifname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:cwld]",
   "language": "python",
   "name": "conda-env-cwld-py"
  },
  "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}