notebooks/OtherIO_options.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 3: Other options for input\n",
    "\n",
    "Umami is a package for calculating metrics for use with for Earth surface dynamics models. This notebook is the final notebook in a three-part introduction to using umami.\n",
    "\n",
    "Umami was designed to work well with the [terrainbento](https://terrainbento.readthedocs.io/en/latest/) model package, as well as other models built using the [Landlab Toolkit](https://github.com/landlab/landlab). However, umami can be used with models built with other modeling tools and data in a variety of formats. This notebook is meant to demonstrate this capability. \n",
    "\n",
    "## Scope of this tutorial\n",
    "\n",
    "In this tutorial you will learn how to use other input options along with umami. \n",
    "\n",
    "Specifically we will use square gridded terrain stored in [ESRI ASCII](http://resources.esri.com/help/9.3/arcgisengine/java/GP_ToolRef/spatial_analyst_tools/esri_ascii_raster_format.htm) format. We will read this in as a numpy array. We will also interpolate it to an irregular grid.\n",
    "\n",
    "If you have comments or questions about the notebooks, the best place to get help is through [GitHub Issues](https://github.com/TerrainBento/umami/issues).\n",
    "\n",
    "To begin this example, we will import the required python packages. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "from io import BytesIO, StringIO\n",
    "\n",
    "import numpy as np\n",
    "from scipy.interpolate import RegularGridInterpolator\n",
    "import matplotlib.pylab as plt\n",
    "from urllib.request import urlopen\n",
    "from urllib.error import URLError\n",
    "\n",
    "import rasterio\n",
    "\n",
    "from landlab import imshow_grid, RasterModelGrid, VoronoiDelaunayGrid\n",
    "\n",
    "from umami import Metric"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Umami does not make any requirements regarding where terrain data comes from or what model or modeling package is used to construct modeled terrain. However, umami does require that modeled or observed terrain is provided to it as a Landlab grid with an at-node field called `topographic__elevation`. \n",
    "\n",
    "Using the Landlab model grid datastructure means that umami knows how large each grid cell is, and how they are connected. Landlab has five model grid classes. One of them will probably suit your needs. \n",
    "\n",
    "- [`RasterModelGrid`](https://landlab.readthedocs.io/en/release/landlab.grid.raster.html#raster)\n",
    "- [`HexModelGrid`](https://landlab.readthedocs.io/en/release/landlab.grid.hex.html#hex)\n",
    "- [`RadialModelGrid`](https://landlab.readthedocs.io/en/release/landlab.grid.radial.html#radial)\n",
    "- [`VoronoiDelaunayGrid`](https://landlab.readthedocs.io/en/release/landlab.grid.voronoi.html#voronoi)\n",
    "- [`NetworkModelGrid`](https://landlab.readthedocs.io/en/release/landlab.grid.network.html#network) \n",
    "\n",
    "In this example we will use the `RasterModelGrid` for regularly spaced square grid cells and the `VoronoiDelaunayGrid` for irregularly spaced observations. \n",
    "\n",
    "If you can read your topography into python as a numpy array, you can put it on a Landlab grid called `topographic__elevation` and use it with umami. While umami is strict about use of the Landlab grid, this quality of the grid makes it very flexible. \n",
    "\n",
    "You can use a Landlab function such as [`read_esri_ascii`](https://landlab.readthedocs.io/en/release/landlab.io.esri_ascii.html#landlab.io.esri_ascii.read_esri_ascii) or [`read_netcdf`](https://landlab.readthedocs.io/en/release/landlab.io.netcdf.html#landlab.io.netcdf.read.read_netcdf) to read your data into a numpy array. You can also create a synthetic one, or use some other package to read a file into python. The world is your oyster. \n",
    "\n",
    "In this case we will use the [rasterio](https://rasterio.readthedocs.io/en/stable/) package to read an [ESRI ASCII](http://resources.esri.com/help/9.3/arcgisengine/java/GP_ToolRef/spatial_analyst_tools/esri_ascii_raster_format.htm) format file that we will download from the [OpenTopography rest server](https://opentopography.org/developers).\n",
    "\n",
    "## Step 1: Read in a numpy array\n",
    "\n",
    "First, we download a small patch of land near Boulder, CO. You can change the values of `north`, `south`, `east` and `west` to change the location. Its not hard to download a very large file, so increase values carefully. \n",
    "\n",
    "These data are provided with horizontal units of degrees. For this example we will not convert from degrees to meters, or address the issue of changing from a geographic coordinate system (WGS84) to a projected one (e.g, UTM Zone 13 N, for Colorado). These are things you should address if you are using this sort of data in an application or research project. \n",
    "\n",
    "The code is wrapped in a `try`-`except` block because if it takes a very long time to get a response from OpenTopography (which sometimes happens on Binder), or if you don't have internet, we want you to still be able to do the tutorial. \n",
    "\n",
    "In this case, you will use some data pre-loaded into the file \"topo_data.asc\". "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    west = -105.4  # longitude (degrees)\n",
    "    east = -105.15  # longitude (degrees)\n",
    "    north = 40.1  # latitude (degrees)\n",
    "    south = 39.9  # latitue (degrees)\n",
    "\n",
    "    URL = \"http://opentopo.sdsc.edu/otr/getdem?demtype=SRTMGL3&\"\n",
    "    url = (URL + \"west=\" + str(west) + \"&\" + \"south=\" + str(south) + \"&\"\n",
    "           \"east=\" + str(east) + \"&\"\n",
    "           \"north=\" + str(north) + \"&\"\n",
    "           \"outputFormat=AAIGrid\")\n",
    "\n",
    "    f = urlopen(url)\n",
    "    file_like = BytesIO(f.read())\n",
    "    print(\"URL Sucess: Using data from OpenTopography.\")\n",
    "except URLError:\n",
    "    print(\"URL Timed out, using pre-saved file.\")\n",
    "    file_like = \"topo_data.asc\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have a variable in our python workspace called `file_like`. We can think of it like a python object that will behave like like an ESRI ASCII file on disk.\n",
    "\n",
    "Next we read `file_like` in with the rasterio package and grab important characteristics like the number of rows (`nrows`), number of columns (`ncols`), the resolution of each pixel (`dx, dy`), the coordinates of the lower left corner (`xy_lower_left`), and the actual elevation data (`elevations`). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with rasterio.open(file_like) as dataset:\n",
    "    nrows = dataset.height\n",
    "    ncols = dataset.width\n",
    "    dx, dy = dataset.res\n",
    "    xy_lower_left = (dataset.bounds.left, dataset.bounds.bottom)\n",
    "    elevations = dataset.read(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot it. As expected, it looks like the topography near Boulder, CO. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(elevations, cmap=\"terrain\", origin=\"lower\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: Create a `RasterModelGrid` to give to Umami\n",
    "\n",
    "Next we create a Landlab model grid by passing the information we got from the rasterio dataset to `RasterModelGrid`. We add the field using the name umami requires, `topographic__elevation`. \n",
    "\n",
    "One tricky step here is that the variable elevations is of datatype `int32`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "elevations.dtype"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some of the underlying tools that umami uses assume that this field is of type `float`. So when we provide the field `topographic__elevation` to the grid, we will specify the that it should be as type `float`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rmg = RasterModelGrid((nrows, ncols),\n",
    "                       xy_spacing=(dx, dy),\n",
    "                       xy_of_lower_left=xy_lower_left)\n",
    "\n",
    "z = rmg.add_field(\"topographic__elevation\", elevations.astype(float))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we use the Landlab function [`imshow_grid`](https://landlab.readthedocs.io/en/release/landlab.plot.html#landlab.plot.imshow.imshow_grid) we see that the topography is correctly represented by the grid.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imshow_grid(rmg, \"topographic__elevation\", cmap=\"terrain\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we make our Metric using the same settings we used in [Part 1](IntroductionToMetric.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metrics = {\n",
    "    \"me\": {\n",
    "        \"_func\": \"aggregate\",\n",
    "        \"method\": \"mean\",\n",
    "        \"field\": \"topographic__elevation\"\n",
    "    },\n",
    "    \"ep10\": {\n",
    "        \"_func\": \"aggregate\",\n",
    "        \"method\": \"percentile\",\n",
    "        \"field\": \"topographic__elevation\",\n",
    "        \"q\": 10\n",
    "    }\n",
    "}\n",
    "\n",
    "rmg_metric = Metric(rmg, metrics=metrics)\n",
    "rmg_metric.calculate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rmg_metric.names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rmg_metric.values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3: Use irregular data and a `VoronoiDelaunayGrid`\n",
    "\n",
    "As a final example, we will look at specifying umami with an irregular grid. We won't import any standard format of irregular data but will create some by interpolating the regular data using the scipy tool [RegularGridInterpolator](https://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.interpolate.RegularGridInterpolator.html).\n",
    "\n",
    "We use a tenth as many grid nodes as we had in the prior example. This is just for speed, feel free to adjust the value for `factor` to change this. \n",
    "\n",
    "We start by creating a set of grid node locations in x and y. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "factor = 11\n",
    "nnodes = int(rmg.x_of_node.size / factor)\n",
    "np.random.seed(27)\n",
    "\n",
    "# select a random subset of x_of_node and y_of_node and permute by a small quantity.\n",
    "# permute only in x, which allows the ordering of nodes to be maintained \n",
    "index = np.linspace(0, rmg.x_of_node.size-1, nnodes, dtype=int)\n",
    "random_x = rmg.x_of_node[index] + 0.4 * rmg.spacing[0] * np.random.randn(index.size)\n",
    "random_y = rmg.y_of_node[index] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot them in comparison with our regular grid nodes. There are a lot of nodes, so we will zoom into a corner of the plot. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(rmg.x_of_node, rmg.y_of_node, 'k.',  markersize=2, label=\"Raster Points\")\n",
    "plt.plot(random_x, random_y, 'm.', label=\"Irregular Points\")\n",
    "plt.xlim(-105.40, -105.375)\n",
    "plt.ylim(40.00, 40.025)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we create an interpolation object and interpolate to find the elevation values at our new randomly located set of model grid nodes based on the regular grid. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interp_obj = RegularGridInterpolator((rmg.y_of_node.reshape(rmg.shape)[:, 0], \n",
    "                                      rmg.x_of_node.reshape(rmg.shape)[0, :]),\n",
    "                                     z.reshape(rmg.shape), bounds_error=False, fill_value=None)\n",
    "\n",
    "interp_z = interp_obj((random_y, random_x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we create a `VoronoiDelaunayGrid` and add `topographic__elevation` to it. \n",
    "\n",
    "One nice feature of the `imshow_grid` function is that it works for both regular and irregular grids. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vdg = VoronoiDelaunayGrid(random_x, random_y)\n",
    "z = vdg.add_field(\"topographic__elevation\", interp_z, at=\"node\")\n",
    "\n",
    "imshow_grid(vdg, z, cmap=\"terrain\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected we see a slightly smoothed version of our original topography. This is expected because we decreased the number of model grid nodes by a factor of 10. \n",
    "\n",
    "The final step is to create a `Metric` and calculate values. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metrics = {\n",
    "    \"me\": {\n",
    "        \"_func\": \"aggregate\",\n",
    "        \"method\": \"mean\",\n",
    "        \"field\": \"topographic__elevation\"\n",
    "    },\n",
    "    \"ep10\": {\n",
    "        \"_func\": \"aggregate\",\n",
    "        \"method\": \"percentile\",\n",
    "        \"field\": \"topographic__elevation\",\n",
    "        \"q\": 10\n",
    "    }\n",
    "}\n",
    "\n",
    "vmg_metric = Metric(vdg, metrics=metrics)\n",
    "vmg_metric.calculate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vmg_metric.names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vmg_metric.values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparing the metric values for the two grids, we can see that the mean is slightly different in absolute value but very close based on percent change and the 10th percentile is identical. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for n in vmg_metric.names:\n",
    "    abs_change = np.abs(vmg_metric.value(n) - rmg_metric.value(n))\n",
    "    pct_change =  abs_change /( 2* (vmg_metric.value(n) + rmg_metric.value(n)))\n",
    "    print(n, \"\\n  abs_change: \", abs_change, \"\\n  pct_change: \", pct_change)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Next steps\n",
    "\n",
    "Now that you have a sense for how the `Metric` and `Residual` classes are used, try the next notebook: [Part 4: Example application](ExampleApplication.ipynb)."
   ]
  }
 ],
 "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.8.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}