dem_processing_with_gdal_python.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Short introduction to working with DEMs in Python GDAL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Greg Tucker, CU Boulder, Feb 2016"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Install GDAL library"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You'll need to install the GDAL library. If you have Anaconda installed, you can do this from the command line by:\n",
    "\n",
    "``conda install gdal``"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Download some DEM data to work with"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Navigate a browser to ``http://criticalzone.org/boulder/data/dataset/2915/``\n",
    "\n",
    "Select *Betasso (Snow off - filtered) - 1m Filtered DSM*\n",
    "\n",
    "Save the zip file, and double-click to unzip it. Inside the folder ``img`` you will see a file called ``czo_1m_bt1.img``. This is a 1 m resolution lidar-derived DEM of a stretch of Boulder Creek Canyon, with the small Betasso tributary catchment located roughly in the center."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import the GDAL library"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we import GDAL and NumPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from osgeo import gdal\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Open and read data from the DEM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Change the path name below to reflect your particular computer, then run the cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "betasso_dem_name = '/Users/gtucker/Dev/dem_analysis_with_gdal/czo_1m_bt1.img'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "geo = gdal.Open(betasso_dem_name)\n",
    "zb = geo.ReadAsArray()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the previous two lines worked, ``zb`` should be a 2D numpy array that contains the DEM elevations. There are some cells along the edge of the grid with invalid data. Let's set their elevations to zero, using the numpy ``where`` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "zb[np.where(zb<0.0)[0],np.where(zb<0.0)[1]] = 0.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's make a color image of the data. To do this, we'll need Pylab and a little \"magic\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.imshow(zb, vmin=1600.0, vmax=2350.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Questions:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Note: to answer the following, open Google Earth and enter Betasso Preserve in the search bar. Zoom out a bit to view the area around Betasso)\n",
    "\n",
    "(1) Use a screen shot to place a copy of this image in your lab document. Label Boulder Creek Canyon and draw an arrow to show its flow direction.\n",
    "\n",
    "(2) Indicate and label the confluence of Fourmile Creek and Boulder Canyon.\n",
    "\n",
    "(3) What is the mean altitude? What is the maximum altitude? (Hint: see numpy functions ``mean`` and ``amax``)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make a slope map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the numpy ``gradient`` function to make an image of absolute maximum slope angle at each cell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def slope_gradient(z):\n",
    "    \"\"\"\n",
    "    Calculate absolute slope gradient elevation array.\n",
    "    \"\"\"\n",
    "    x, y = np.gradient(z)  \n",
    "    #slope = (np.pi/2. - np.arctan(np.sqrt(x*x + y*y)))\n",
    "    slope = np.sqrt(x*x + y*y)\n",
    "    return slope"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sb = slope_gradient(zb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see what it looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.imshow(sb, vmin=0.0, vmax=1.0, cmap='pink')\n",
    "print np.median(sb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Questions:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(1) Place a copy of this image in your lab document. Identify and label the Betasso Water Treatment plant.\n",
    "\n",
    "(2) How many degrees are in a slope gradient of 1.0 (or 100%)?\n",
    "\n",
    "(3) What areas have the steepest slopes? What areas have the gentlest slopes? What do you think the distribution of slopes might indicate about the distribution of erosion rates within this area?\n",
    "\n",
    "(4) What is the median slope gradient? What is this gradient in degrees? (Hint: numpy has a ``median`` function)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make a map of slope aspect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def aspect(z):\n",
    "    \"\"\"Calculate aspect from DEM.\"\"\"\n",
    "    x, y = np.gradient(z)\n",
    "    return np.arctan2(-x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ab = aspect(zb)\n",
    "plt.imshow(ab)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can make a histogram (frequency diagram) of aspect. Here 0 degrees is east-facing, 90 is north-facing, 180 is west-facing, and 270 is south-facing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "abdeg = (180./np.pi)*ab # convert to degrees\n",
    "n, bins, patches = plt.hist(abdeg.flatten(), 50, normed=1, facecolor='green', alpha=0.75)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Questions:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(1) Place a copy of this image in your lab notes.\n",
    "\n",
    "(2) Compare the aspect map to imagery in Google Earth. Is there any correlation aspect and vegetation? If so, what does it look like?\n",
    "\n",
    "(3) What is the most common aspect? (N, NE, E, SE, S, SW, W, or NW)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Shaded relief"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a shaded relief image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def hillshade(z, azimuth=315.0, angle_altitude=45.0):  \n",
    "    \"\"\"Generate a hillshade image from DEM.\n",
    "    \n",
    "    Notes: adapted from example on GeoExamples blog,\n",
    "    published March 24, 2014, by Roger Veciana i Rovira.\n",
    "    \"\"\"\n",
    "    x, y = np.gradient(z)  \n",
    "    slope = np.pi/2. - np.arctan(np.sqrt(x*x + y*y))  \n",
    "    aspect = np.arctan2(-x, y)  \n",
    "    azimuthrad = azimuth*np.pi / 180.  \n",
    "    altituderad = angle_altitude*np.pi / 180.\n",
    "   \n",
    "    shaded = np.sin(altituderad) * np.sin(slope)\\\n",
    "     + np.cos(altituderad) * np.cos(slope)\\\n",
    "     * np.cos(azimuthrad - aspect)  \n",
    "    return 255*(shaded + 1)/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "hb = hillshade(zb)\n",
    "plt.imshow(hb, cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Questions:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(1) Place a copy of this image in your lab document. Label at least one area of relatively smooth terrain, and one area of relatively rough terrain.\n",
    "\n",
    "(2) Do the areas of smooth and rough topography bear any relation to other geomorphic features?\n",
    "\n",
    "(3) Note in your document any other observations, comments, or questions that have occurred to you while examining these data."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}