Python-and-Gdal.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Raster data processing with Python and GDAL\n",
    "\n",
    "In this lesson we will learn basic techniques for automating raster data processing with python and the **Geospatial Data Abstraction Library GDAL**. \n",
    "\n",
    "Raster files are commonly used to store terrain models and remote sensing data and their derivate products such as vegetation indices and other environmental datasets. Raster files tend to be huge (imagine for example a raster dataset covering the globe in 30m x 30m resolution) and are often delivered and processed in smaller pieces (tiles). Efficient processing and analysis of such large datasets consequently requires automatization. During this lesson you will learn how to read and write common raster formats, and conduct basic raster data processes for a batch of files using the GDAL/OGR API in Python and GDAL command line utilities. \n",
    "\n",
    "If you want to dive deeper into raster data processing in Python you might want to check out [Rasterio](https://mapbox.github.io/rasterio/index.html) which is a relatively new Python library for efficient raster data analysis. Also remember that you can write Python scripts for using ArcGIS tools as introduced in the [demo](https://automating-gis-processes.github.io/2016/Lesson6-arcpy-script.html). \n",
    "\n",
    "For more information and examples with GDAL, please see these online resources:\n",
    "\n",
    "- http://www.gdal.org/gdal_tutorial.html\n",
    "- https://pypi.python.org/pypi/GDAL/\n",
    "- https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Extract files\n",
    "\n",
    "**Download data for this lesson from [this link](http://www.helsinki.fi/science/accessibility/opetus/autogis/LC81910182016153-SC20161208043748.tar.gz).**\n",
    "\n",
    "In this lesson, we will use a Landsat 8 satellite scene for practicing how to read raster data with Python and Gdal. Landsat images are distributed as gzipped TAR archives `*.tar.gz` when downloading from the [USGS Earth Explorer](https://earthexplorer.usgs.gov/). For the purposes of this exercise, we have pre-downloaded a scene covering South-West Finland which you can download from the above link. In order to get started, you need to extract the contents of the TAR archive either in the Terminal window or with a Python script.\n",
    "\n",
    "\n",
    "#### Option 1: Extract files in the Terminal window\n",
    "\n",
    "Navigate to the directory where you downloaded the data, and extract files into a new folder:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mkdir LandsatData\n",
    "tar -zxvf LC81910182016153-SC20161208043748.tar.gz -C LandsatData/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Option 2: Extract files using python\n",
    "\n",
    "It is also possible to extract the files in python using the `tarfile` module. First, open the file with `tarfile.open()`. The parameter `r:gz` specifies that we want to open the gzipped file in reading mode. Then, extract the files using `tarifile.extractall()` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import tarfile\n",
    "\n",
    "#Create output folder\n",
    "newFolder = \"LandsatData\"\n",
    "os.makedirs(newFolder)\n",
    "\n",
    "#Extract files\n",
    "tar = tarfile.open(\"LC81910182016153-SC20161208043748.tar.gz\", \"r:gz\")\n",
    "tar.extractall(newFolder)\n",
    "tar.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 1. Open a raster file in Python with GDAL\n",
    "\n",
    "With GDAL, you can read and write [several different raster formats](http://www.gdal.org/formats_list.html) in Python. Python automatically registers all known GDAL drivers for reading supported formats when the importing the GDAL module. Most common file formats include for example [TIFF and GeoTIFF](http://www.gdal.org/frmt_gtiff.html), [ASCII Grid](http://www.gdal.org/frmt_various.html#AAIGrid) and [Erdas Imagine .img -files](http://www.gdal.org/frmt_hfa.html).    \n",
    "\n",
    "Landsat 8 bands are stored as separate GeoTIFF -files in the original package. Each band contains information of surface reflectance from different ranges of the electromagnetic spectrum. You can read more information about Landsat 8 [here](http://landsat.gsfc.nasa.gov/landsat-8/landsat-8-bands/).\n",
    "\n",
    "Let's start with inspecting one of the files we downloaded:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from osgeo import gdal\n",
    "\n",
    "filepath = r\"LandsatData/LC81910182016153LGN00_sr_band4.tif\"\n",
    "\n",
    "# Open the file:\n",
    "raster = gdal.Open(filepath)\n",
    "\n",
    "# Check type of the variable 'raster'\n",
    "type(raster)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read raster file properties\n",
    "\n",
    "The satellite image is now stored as a GDAL Dataset object in the variable `raster`. Let's have a closer look at the properties of the file: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Projection\n",
    "raster.GetProjection()\n",
    "\n",
    "# Dimensions\n",
    "raster.RasterXSize\n",
    "raster.RasterYSize\n",
    "\n",
    "# Number of bands\n",
    "raster.RasterCount\n",
    "\n",
    "# Metadata for the raster dataset\n",
    "raster.GetMetadata()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get raster bands\n",
    "\n",
    "In our case, all bands of the Landsat 8 scene are stored as separate files. `rasterCount` is 1 as we have only opened one GeoTiff containing Landsat 8 band 4. However, different bands of a satellite images are often stacket together in one raster dataset in which case `rasterCount` would be greater than one.  \n",
    "\n",
    "In order to have a closer look at the values stoden in the band, we will take advantage of the [GDAL Band API](http://gdal.org/python/osgeo.gdal.Band-class.html). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Read the raster band as separate variable\n",
    "band = raster.GetRasterBand(1)\n",
    "\n",
    "# Check type of the variable 'band'\n",
    "type(band)\n",
    "\n",
    "# Data type of the values\n",
    "gdal.GetDataTypeName(band.DataType)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have a GDAL Raster Band object stored in the variable `band`. \n",
    "\n",
    "Data type of the band can be interpreted with the help of GDAL documentation on [Pixel data types](http://www.gdal.org/gdal_8h.html#a22e22ce0a55036a96f652765793fb7a4). Unsigned integer is always equal or greater than zero and signed integer can store also negative values. For example, an unsigned 16-bit integer can store 2^16 (=65,536) values ranging from 0 to 65,535.\n",
    "\n",
    "### Band  statistics\n",
    "Next, let's have a look at the values that are stored in the band. You might need to calculate statistics for the raster before being able to print out any information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Compute statistics if needed\n",
    "if band.GetMinimum() is None or band.GetMaximum()is None:\n",
    "    band.ComputeStatistics(0)\n",
    "    print(\"Statistics computed.\")\n",
    "    \n",
    "# Fetch metadata for the band\n",
    "band.GetMetadata()\n",
    "    \n",
    "# Print only selected metadata:\n",
    "print (\"[ NO DATA VALUE ] = \", band.GetNoDataValue()) # none\n",
    "print (\"[ MIN ] = \", band.GetMinimum())\n",
    "print (\"[ MAX ] = \", band.GetMaximum())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Open raster file as a numerical array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "GDAL is a powerful library when it comes to accessing geospatial raster data, but it does not provide many functionalities for doing calculations. For more advanced computing, we will read in the raster data as a numerical array in order to use the capabilities in the NumPy-library.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In case you want to convert and existing Gdal Dataset or a Band into a numpy array you can convert it with `ReadAsArray` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# read raster data as numeric array from Gdal Dataset\n",
    "rasterArray = raster.ReadAsArray()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "so what's the difference?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Check the datatype of variables\n",
    "type(rasterArray)\n",
    "type(raster)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, our we have now stored the same raster data into two different types of variables:\n",
    "\n",
    "-`raster` is a Gdal Dataset \n",
    "\n",
    "-`rasterArray` is a numpy array "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The GDAL library also comes with a module gdal_array that works as an interface between NumPy and GDAL. Gdal_array reads and writes raster files to and from NumPy arrays directly from file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from osgeo import gdal_array\n",
    "#from osgeo import gdalnumeric\n",
    "\n",
    "# Read raster data as numeric array from file\n",
    "rasterArray = gdal_array.LoadFile(filepath)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Take care of NoData values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# What is the minimum value of the array?\n",
    "rasterArray.min() #-9999"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the numpy array still contains the original nodata values. Any calculations will be wrong if these are not taken care of. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Get nodata value from the GDAL band object\n",
    "nodata = band.GetNoDataValue()\n",
    "\n",
    "#Create a masked array for making calculations without nodata values\n",
    "rasterArray = np.ma.masked_equal(rasterArray, nodata)\n",
    "type(rasterArray)\n",
    "\n",
    "# Check again array statistics\n",
    "rasterArray.min()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you have a two-dimensional array ready for further calculations. However, for completing our exercise for this week we will use a very simple command line tool `gdal_calc.py`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Close raster dataset\n",
    "\n",
    "It might be useful to close an exitsing GDAL object in the middle of the code if you want to free resources and remove unnecessary varables from memory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "raster = None\n",
    "band = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is however not necessary to close raster datasets or bands at the end of the python script as Pyhton will automatically take care of this."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. GDAL command line utilities \n",
    "\n",
    "We have now tested some of the basic functions from the Python [GDAL/OGR API](https://pcjericks.github.io/py-gdalogr-cookbook/index.html) for reading and inspecting raster files. However, GDAL also includes other powerful functions for data translation and processing which are not directly implemented in the library. We will have a closer look on a couple of such functions:\n",
    "\n",
    "- [gdalwarp](http://www.gdal.org/gdalwarp.html) for clipping, mosaicing, reprojection and other processes\n",
    "- [gdal_merge.py](http://www.gdal.org/gdal_merge.html) for mosaicing / stacking images\n",
    "- [gdal_calc.py](http://www.gdal.org/gdal_calc.html) for raster calculations\n",
    "\n",
    "\n",
    "These tools need to be run from the Terminal/Command Prompt or as a subprocess in Python. We will now quickly test out these tools in the Terminal window."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clipping the image with gdalwarp\n",
    "\n",
    "Among other tricks, gdalwarp is a very handy tool for quickly clipping your image. We will now practice how to clip the satellite image band based on a bounding box. Desired extent for the output file is specified using the option `-te` :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "gdalwarp -te xmin ymin xmax ymax inputfile.tif outputfile.tif"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ---------------------------------------------\n",
    "### TASK: Clip band 4 of the satellite image so that it covers cloud-free areas in the Turku archipelago (in the North-East corner of the scene). \n",
    "\n",
    "### You can open the band 4 manually in QGIS for defining the corner coordinates.\n",
    "### ---------------------------------------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's repeat the clipping for all the rest of the bands all at once. For doing this, we will use Python for generating the command for each spectral band in our scene."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import glob\n",
    "import os\n",
    "\n",
    "# List filepaths for all bands in the scence\n",
    "FileList = glob.glob(os.path.join(r'/home/geo/LandsatData','*band*.tif'))\n",
    "\n",
    "# Define clipping extent\n",
    "xmin, ymin, xmax,ymax = (0, 0, 0, 0) # INSERT HERE THE CORRECT COORDINATES\n",
    "\n",
    "# Generate gdalwarp command for each band\n",
    "command = \"\"\n",
    "\n",
    "for fp in FileList:\n",
    "    inputfile = fp\n",
    "    outputfile = inputfile[:-4] + \"_clip.tif\"\n",
    "    \n",
    "    command += \"gdalwarp -te %s %s %s %s %s %s \\n\" % (xmin, ymin, xmax, ymax, inputfile, outputfile)\n",
    "    \n",
    "# Write the commands to an .sh file\n",
    "cmd_file = \"ClipTurkufromLandsat.sh\"\n",
    "f = open(os.path.join(cmd_file), 'w')\n",
    "\n",
    "f.write(command)\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NOTE: If you are working in an windows environment, change the `.sh` extension to `.bat`. \n",
    "\n",
    "After running the above script, you should have a file \"ClipTurkufromLandsat.sh\" in your working directory. Open the file (with a text editor) and check that the commands have been written correctly to the file. \n",
    "\n",
    "Next, run the file in the Terminal window:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "bash ClipTurkufromLandsat.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you should have a bunch of clipped .tif files ready and you might want to open a few of them in QGIS to check that the process was succesful."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generating a layer stack with gdal_merge.py\n",
    "\n",
    "After clipping the image you can for example stack bands 3 (green), 4 (red), and 5 (nir) for visualizing a false-color composite. Merge the layers with `gdal_merge.py` and use the `-separate` option for indicating that you wish to save the inputs as separate bands in the output file.\n",
    "\n",
    "Let's try running the command as a subprocess in python: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "# Define input and output files\n",
    "inputfiles =  \"band3_clip.tif band4_clip.tif band5_clip.tif\"\n",
    "outputfile =  \"Landsat8_GreenRedNir.tif\"\n",
    "\n",
    "# Generate the command\n",
    "command = \"gdal_merge.py -separate %s -o %s\" % (inputfiles, outputfile)\n",
    "\n",
    "# Run the command. os.system() returns value zero if the command was executed succesfully\n",
    "os.system(command)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a result, you have three bands stacked together in the file `Landsat8_GreenRedNir.tif`. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gdal_calc.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Gdal_calc.py` is a command line raster calculator which can be useful for competing simple repetitive calculations for raster data. **Open the terminal window** and type in:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Gdal_calc.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And press enter. You should see instructions on the usage and options for the tool.The basic syntax for `gdal_calc.py` is the following:\n",
    "\n",
    "`gdal_calc.py -A input1.tif - B input2.tif [other_options] --outfile=outputfile.tif`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From other options, it is useful to notice at least the parameters `--calc` for specifying the calculation syntax and `--creation-option` (or `--co`) for controlling the output file size:\n",
    "\n",
    "- In the case of two input files `--calc=\"A+B\"` would add files A and B together. \n",
    "- By default output files tend to be huge which will quickly result in problems with disk size and memory. With `gdal_calc.py` you can add parameter `--co=\"COMPRESS=LZW\"` in order to reduce output file size.\n",
    "\n",
    "\n",
    "## Exercise 7: Delineate forest area from canopy cover data\n",
    "\n",
    "In [exercise 7](https://classroom.github.com/assignment-invitations/e59dfd42577e7fe6d98d19d70ecccf53) we will work with a piece of the [Global Forest Change Dataset](http://earthenginepartners.appspot.com/science-2013-global-forest) by Hansen et al (2013). The goal of the exercise is to delineate forest areas from an area of interest by processing raster layers containing information about canopy cover in 2000. \n",
    "\n",
    "You should be able to complete the exercise by altering the simple tools presented during this lesson. "
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}