2e_spatial_interpolation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div>\n",
    "<table style=\"width: 100%\">\n",
    "\t<tr>\n",
    "\t\t<td>\n",
    "\t\t<table style=\"width: 100%\">\n",
    "\t\t\t<tr>\n",
    "                <td ><center><font size=\"5\"><b>Module 49</b></font><center>\n",
    "                <center><font size=\"6\">Digital Innovations for Water Challenges</font><center></td>\n",
    "\t\t\t</tr>\n",
    "\t\t\t<tr>\n",
    "                <td><center><font size=\"14\">Notebook 2.e</font><center></td>\n",
    "\t\t\t</tr>\n",
    "\t\t\t<tr>\n",
    "                <td><center><font size=\"6\"><b>Spatial Interpolation</b></font><center></td>\n",
    "\t\t\t</tr>\n",
    "\t\t</table>\n",
    "\t\t</td>\n",
    "\t\t<td><center><img src='images\\ihe-delft-institute_unesco_fc-lr.jpg'></img></td>\n",
    "\t</tr>\n",
    "</table>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Table of contents\n",
    "1. [Learning objectives](#learningobs)\n",
    "2. [Introduction](#introduction)\n",
    "3. [Convert a CSV file to PCRaster map format](#csv)\n",
    "4. [Interpolate borehole depth with IDW](#idw)\n",
    "5. [Interpolate borehole depth with Thiessen polygons](#thiessen)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Learning objectives<a name=\"learningobs\"></a>\n",
    "\n",
    "- To convert CSV files to rasters\n",
    "- To interpolate points with Inverse Distance Weighting (IDW)\n",
    "- To interpolate points with Thiessen polygons"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Introduction<a name=\"introduction\"></a>\n",
    "In this exercise you'll learn how to convert a Comma Separated Values file with borehole depth and coordinates to a PCRaster map. Next, we're going to interpolate the borehole depths to continuous raster using the Inverse Distance Weighting (IDW) method and the Thiessen method.\n",
    "\n",
    "Let's first import PCRaster and set the path to the `2e_Data` folder where we have the input data for this exercise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pcraster import *\n",
    "import os\n",
    "print(os.getcwd())\n",
    "os.chdir(\"./2e_Data\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Convert a CSV file to PCRaster map format<a name=\"csv\"></a>\n",
    "Open the `boreholesdhxy.csv` file in the `2e_Data` folder in a text editor and check the contents.\n",
    "Which columns are the coordinates and which column gives the depth of the borehole?\n",
    "\n",
    "PCRaster has a few command line applications that are useful for these tasks. You can find these applications [here](https://pcraster.geo.uu.nl/pcraster/latest/documentation/pcraster_manual/sphinx/applications.html). We will use <code>col2map</code> to convert the CSV file to a PCRaster map.\n",
    "\n",
    "The syntax of <code>col2map</code> is <code>col2map [options] columnfile PCRresult</code>.\n",
    "Here we'll use some options to indicate the output data type (Scalar) and order of the x, y and data column (v). We also need to indicate a clone map that has the raster properties (coordinate system, cell size, nr of rows/columns) of the desired output.\n",
    "\n",
    "We can run command line commands from Python using <code>os.system(stringWithCommand)</code>:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "command = \"col2map -S -x 2 -y 3 -v 1 boreholesdxy.csv boreholedepth.map --clone clone.map\"\n",
    "os.system(command)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read the `boreholedepth.map` file into the variable `BoreholeDepth`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualise `BoreholeDepth` with Aguila. Bacause we run this Jupyter Notebook locally, we can use the `aguila` Python operation. The map will popup in the background, so you might need to look for it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aguila(BoreholeDepth)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Interpolate borehole depth with IDW<a name=\"idw\"></a>\n",
    "We can now interpolate the borehole depth at the pixels. We'll first apply the Inverse Distance Weighting method. This has been implemented as the [inversedistance operation](https://pcraster.geo.uu.nl/pcraster/4.3.0/documentation/pcraster_manual/sphinx/op_inversedistance.html) in PCRaster.\n",
    "\n",
    "The syntax is: \n",
    "```Python\n",
    "Result = inversedistance(mask, points, idp, radius, maxNr)\n",
    "```\n",
    "\n",
    "Where \n",
    "* <code>mask</code> is a mask of the extent that will be interpolated. In our case we'll use the clone map. \n",
    "* <code>points</code> is the scalar raster layer with pixels that you want to interpolate. In our case <code>BoreholeDepth</code>.\n",
    "* <code>idp</code> is the exponent of the power function for the weights. An <code>idp</code> value of 2 is commonly used.\n",
    "* <code>radius</code> indicates the distance to a cell within it has to take interpolation points into account. If it's set to 0 there's no limitation.\n",
    "* <code>maxNr</code> limits the number of points used for the interpolation. If it's set to 0 there's no maximum.\n",
    "\n",
    "Let's first read the clone map that we'll use as the mask."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Mask = readmap(\"clone.map\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we interpolate the borehole depths with IDW using the following line:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "BoreholeDepthIDW = inversedistance( Mask, BoreholeDepth, 2, 0, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualise the result with Aguila."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. Interpolate borehole depth with Thiessen polygons<a name=\"thiessen\"></a>\n",
    "Another interpolation method is using so called Thiessen polygons.\n",
    "\n",
    "In PCRaster this involves three steps:\n",
    "1. Assign a unique ID to each borehole\n",
    "2. Assign to each pixel the closest borehole ID\n",
    "3. Assign the borehole depth to the corresponding borehole IDs for all cells\n",
    "\n",
    "Let's first assign a unique ID to each borehole using the <code>uniqueid</code> operation. The input needs to be a boolean raster.\n",
    "\n",
    "To create a boolean raster of the boreholes we use the <code>defined</code> operation. This will give boolean True for all cells with a value and boolean False for all nodata cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Boreholes = defined(BoreholeDepth)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the result with Aguila."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we assign a unique ID to each boolean True cell. Note that we want the output to be nominal instead of scalar. Therefore we wrap the <code>uniqueid</code> into the <code>nominal</code> operation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "BoreholeID = nominal(uniqueid(Boreholes))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the result with Aguila."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aguila(BoreholeID)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's assign to each cell its closes borehole ID. We use the <code>spreadzone</code> operation, which has the following syntax:\n",
    "\n",
    "<code>Result = spreadzone(points, initialfrictiondist, friction)</code>\n",
    "\n",
    "Where:\n",
    "\n",
    "* <code>points</code> is the raster with points to interpolate. In our case BoreholeID.\n",
    "* <code>initialfrictiondist</code> is the initial friction distance. We don't need that, so we make it 0.\n",
    "* <code>friction</code> is the friction for each cell, which we keep constant at 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ThiessenID = spreadzone(BoreholeID,0,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the result with Aguila."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aguila(ThiessenID)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final step is to assign the groundwater depth values of the boreholes to the Thiessen polygons.\n",
    "We can do that by calculating the minimum or maximum depth per class, which will just extrapolate the value to the whole class. Here we'll use <code>areamaximum</code>:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "BoreholeDepthThiessen = areamaximum(BoreholeDepth,ThiessenID)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the result and compare it with the IDW interpolation using Aguila."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aguila(BoreholeDepthIDW,BoreholeDepthThiessen)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can you write a script that calculates Thiessen polygons? Make a function that can be applied to similar data as the `BoreholeDepth` raster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}