Batch_process.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Key Methods: Week 10\n",
    "\n",
    "# Spatial data analysis and batch processing in Python and GIS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run this to display graphs nicely in this notebook\n",
    "\n",
    "%matplotlib inline\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "This week we will use a mixture of python and GIS to process a large set of data. GIS allows us to visually assess and interact with each step, whilst python allows us to easily scale up to processing large datasets with a single button press, making repeat processing easy. We will answer the same question as week 2; using lidar data to map biomass along the Spey river valley.\n",
    "\n",
    "Tasks for you to complete are written in **bold font**.\n",
    "\n",
    "\n",
    "## Data\n",
    "\n",
    "Again we are using the DTM and DSM provided by the [Scottish government](https://remotesensingdata.gov.scot/data#/list). These have been provided through Noteable in the ***data*** directory. The raw data comes in 1 m resolution, but as we saw in week 2, the ALS point density was a little too low to support that, so we coarsened it to 3 m, taking the mean elevation for the DTM and the maximum elevation for the DSM. Due to disk space limits on noteable, the data has already been coarsened to 3 m resolution to remove the gaps between ALS points, and then a 10 m mean product calculated.\n",
    "\n",
    "The 10 m resolution DTM and DSMs are available in the folders below, in geotiff raster format:\n",
    "\n",
    "    data/ALS/DTM\n",
    "    data/ALS/DSM\n",
    "\n",
    "Later on a python cell will be given to show you how that coarsening was done, and a single full-res DSM file has been provided in in case you want to run this process from first principles. Note that you can install Python on your own laptop, or use it on the University's Windows or Linux computers, where you will get access to move more disk space.\n",
    "\n",
    "\n",
    "You can download some of these tiles to your computer and view them in QGIS if you would like. Or you can use some of the plotting functions within python to view them here.\n",
    "\n",
    "This week we will use the GDAL package to open our data. Below you will see code which uses this module and its functions to open and read GeoTIFF files. Note that `rasterio`, which we used last week, provides functions which do the same. *There are many many many Python modules out there, some with overlapping functionality.* However, `gdal` makes it easier to *write GeoTIFFs to file*, which is more difficult with `rasterio`.\n",
    "\n",
    "The cell below provides a function to open a single raster file. This is an example of a **User Defined Function**, something we did not cover in previous practicals due to lack of time. After you run the cell, the function can be called just like any of the other functions we use (without needing to import it). See it called toward the bottom of the cell to understand its arguments (a filename) and return value list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gdal,osr\n",
    "import numpy as np\n",
    "from sys import exit\n",
    "\n",
    "def readTiff(filename):\n",
    "    '''A function to read a geotiff in to a numpy array'''\n",
    "\n",
    "    # open the input raster file\n",
    "    ds=gdal.Open(filename)\n",
    "\n",
    "    # check that the file is open\n",
    "    if ds is None:\n",
    "        print('Failed to open',filename)\n",
    "        exit\n",
    "    \n",
    "    # read all geolocation information.\n",
    "    # This tells the computer how to map a 2D coordinate on to the Earth\n",
    "    proj=osr.SpatialReference(wkt=ds.GetProjection())\n",
    "    epsg=int(proj.GetAttrValue('AUTHORITY',1))\n",
    "\n",
    "    # get the raster size\n",
    "    nX=ds.RasterXSize             # number of pixels in x direction\n",
    "    nY=ds.RasterYSize             # number of pixels in y direction\n",
    "    \n",
    "    # geolocation tiepoint\n",
    "    transform_ds = ds.GetGeoTransform()# extract geolocation information\n",
    "    xOrigin=transform_ds[0]       # coordinate of x corner\n",
    "    yOrigin=transform_ds[3]       # coordinate of y corner\n",
    "    pixelWidth=transform_ds[1]    # resolution in x direction\n",
    "    pixelHeight=transform_ds[5]   # resolution in y direction\n",
    "\n",
    "    # read data. Returns as a 2D numpy array\n",
    "    data=ds.GetRasterBand(1).ReadAsArray(0,0,nX,nY)\n",
    "    \n",
    "    # set no data values to 0, otherwise they will cause issues.\n",
    "    missingInd=np.where(data<-100.0)\n",
    "    if(len(missingInd)>0):\n",
    "        data[missingInd[0],missingInd[1]]=0.0\n",
    "        \n",
    "    # pass data back\n",
    "    return(data,xOrigin,yOrigin,pixelWidth,pixelHeight,nX,nY)\n",
    "\n",
    "\n",
    "\n",
    "# set a filename\n",
    "filename='data/ALS/DSM_1m/NN89_1M_DSM_PHASE1.tif'\n",
    "\n",
    "# use the function above to read the data in to RAM\n",
    "data,xOrigin,yOrigin,pixelWidth,pixelHeight,nX,nY=readTiff(filename)\n",
    "\n",
    "# here data is a numpy array containing the DSM data\n",
    "# xOrigin tells us the x coordinate of the left side of the data\n",
    "# yOrigin tells us the y coordinate of the top side of the data (geotiffs read from top to bottom)\n",
    "# pixelWidth tells use the size of each pixel in the x direction\n",
    "# pixelHeight tells use the size of each pixel in the y direction. Will be negative if from top to bottom\n",
    "# nX tells us how many x columns there are\n",
    "# nY tells us how many y rows there are\n",
    "\n",
    "# let us pring some of these to screen to see them\n",
    "print('pixelHeight',pixelHeight)\n",
    "print('The raster has',nX,'by',nY,'pixels')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now our data is in RAM, we can manipulate it in anyway we like. We can use the numpy functions to calculate some summary statistics of the data, as below.\n",
    "\n",
    "**Add some lines to the end of the cell to calculate and print the minimum and maximum elevation for this tile.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mean elevation\n",
    "meanElev=np.mean(data)\n",
    "print('Mean elevation is',meanElev,'m')\n",
    "\n",
    "# some measure of roughness\n",
    "stdevElev=np.std(data)\n",
    "print('The standard deviation of elevation is',stdevElev,'m')\n",
    "\n",
    "# calculate and print the min and max elevation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Image maps\n",
    "\n",
    "It is very useful to be able to plot images of our data to quicly asses it and check for artefacts. Simply calling `plt.imshow(data)` will now make an image map, but there are several things that we will want to change about this default map.\n",
    "\n",
    "First, note that the numbers on the y axis are in reverse order (and, if you check this tile in QGIS you will see that the map is inverted). This follows the convention for storing images (i.e. a digital photograph will be shown correctly) but can be corrected for map axes with an origin in the lower left corner by adding optional (keyword) argument `origin='lower'`. \n",
    "\n",
    "<span style=\"color:blue\">*Note*: </span>In Week 9, we did not need to provide this optional argument, because the data was already aligned in the correct order for plotting with `imshow()` (i.e. the last row in the array corresponded to the southernmost location). All we needed to do was swap \"bottom\" and \"top\" in the *extent* optional (keyword) argument so that the axes displayed georeferencing information correctly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# plot inverted  \n",
    "plt.imshow(data)   # plot the 2D image\n",
    "plt.show()         # print to screen \n",
    "plt.clf()          # clear out plt ready for next plot \n",
    "\n",
    "\n",
    "# flip the y axis to put it the right way around\n",
    "plt.imshow(data,origin='lower')\n",
    "plt.show()         # print to screen \n",
    "plt.clf()          # clear out plt ready for next plot "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ***Optional:*** Pre-processing data\n",
    "\n",
    "This part is included for your own interest. You can skip to the \"Mapping biomass\" part if you like. The data has already been processed to a 10 m resolution DTM and DSM ready for you to use to map biomass. This is included here if you would like to see how to process the data from first principles, but it goes in to the fine details of the geolocation information within the raster file.\n",
    "\n",
    "Remember that in week 2, we visually inspected the data and noticed that the DSMs had stripes between the lidar scan lines, due to the DSM being produced at too high a resolution. We need to coarsen the DSM to around 3 m resolution, taking the maximum so that the false zero pixels do not bias the canopy height. This can be done with the script below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is taken from:\n",
    "# https://gis.stackexchange.com/questions/110769/gdal-python-aggregate-raster-into-lower-resolution\n",
    "\n",
    "# libraries needed to handle the data\n",
    "import gdal,osr\n",
    "\n",
    "\n",
    "def coarsenRaster(filename,outRes,outName):\n",
    "    '''Coarsen the resolution of a raster and write to file'''\n",
    "\n",
    "    inDS=gdal.Open(filename)\n",
    "\n",
    "\n",
    "    # reproject that data to a coarser resolution\n",
    "    # first read the geolocation information, which includes the resolution\n",
    "    outGeo=list(inDS.GetGeoTransform())  # puts the geolocation information in to a python list\n",
    "\n",
    "    # read each element of the list in to variables we can use to rescale\n",
    "    x0=outGeo[0]   # the x corner\n",
    "    xRes=outGeo[1] # x resolution\n",
    "    i0=outGeo[2]   # x pixel number of the corner, normally 0\n",
    "    y0=outGeo[3]   # the y corner\n",
    "    j0=outGeo[4]   # y pixel number of the corner, normally 0\n",
    "    yRes=outGeo[5] # y resolution\n",
    "    nX=inDS.RasterXSize  # number of x pixels\n",
    "    ny=inDS.RasterYSize  # number of y pixels\n",
    "\n",
    "    # coarsen the above by the factor you want\n",
    "    outNx=int(nX*xRes/outRes)  # number of x pixels in the output\n",
    "    outNy=int(nY*abs(yRes)/outRes)  # number of y pixels in the output. Abs makes the negative positive\n",
    "    outXres=outRes\n",
    "    outYres=-1*outRes\n",
    "    \n",
    "    # load the new geoloction information in to a python tuple\n",
    "    OutGT = tuple([x0,outXres,i0,j0,y0,outYres])\n",
    "    \n",
    "    # set up a new geotiff dataset\n",
    "    OutDriver = gdal.GetDriverByName('GTiff')\n",
    "\n",
    "    # load the geolocation in to the new file\n",
    "    outDS = OutDriver.Create(outName,outNx,outNy,1,gdal.GDT_Float32 )\n",
    "    outDS.SetGeoTransform(OutGT)\n",
    "    outDS.SetProjection(inDS.GetProjection()) # copy projection info\n",
    "\n",
    "    # read the data\n",
    "    inData=inDS.GetRasterBand(1).ReadAsArray(0,0,nX,nY)\n",
    "    \n",
    "    # coarsen it\n",
    "    outData=np.zeros((outNx,outNy),dtype=np.float)\n",
    "    xWindow=int(outXres/xRes)/2\n",
    "    yWindow=int(outYres/yRes)/2\n",
    "    nXratio=nX/outNx\n",
    "    nYratio=nY/outNy\n",
    "    \n",
    "    # loop over the x and y pixels of the coarse image and read the original image pixels\n",
    "    for i in range(0,outNx):\n",
    "        \n",
    "        # set x slicing bounds\n",
    "        minI=int(i*nXratio-xWindow)\n",
    "        if(minI<0):\n",
    "            minI=0\n",
    "        maxI=int(i*nXratio+xWindow+1)\n",
    "        if(maxI>=nX):\n",
    "            maxI=nX-1\n",
    "            \n",
    "        for j in range(0,outNy):\n",
    "            \n",
    "            # set slicing bounds\n",
    "            minJ=int(j*nYratio-yWindow)\n",
    "            if(minJ<0):\n",
    "                minJ=0\n",
    "            maxJ=int(j*nYratio+yWindow+1)\n",
    "            if(maxJ>=nY):\n",
    "                maxJ=nY-1\n",
    "            \n",
    "            # calculate the maximum pixel within the sliced bounds\n",
    "            outData[i,j]=np.max(inData[minI:maxI,minJ:maxJ])\n",
    "\n",
    "    # write the data to the geotiff\n",
    "    outDS.GetRasterBand(1).WriteArray(outData)  # write image to the raster\n",
    "    outDS.GetRasterBand(1).SetNoDataValue(0)  # set no data value\n",
    "    outDS.FlushCache()                     # write to disk\n",
    "    outDS = None\n",
    "\n",
    "    # tell the user where the data has been written to\n",
    "    print('New file written to',outName)\n",
    "    \n",
    "    \n",
    "#####################################\n",
    "# now call the above function\n",
    "    \n",
    "# define an output resolution\n",
    "outRes=3   # 3 m resolution\n",
    "\n",
    "# set an input and output filename\n",
    "inName='data/ALS/DSM_1m/NN89_1M_DSM_PHASE1.tif'\n",
    "outName='data/ALS/DSM_1m/NN89_3M_DSM_PHASE1.tif'\n",
    "\n",
    "# call the above function to coarsen and write\n",
    "coarsenRaster(inName,outRes,outName)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By changing the ***np.max()*** command to ***np.mean()***, the same script can be used to coarsen both the 3 m resolution DSM we have just created and the 1 m resolution DTM to create a 10 m resoltion DSM and DTM. **You can use the cell below to do that if you would like.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mapping biomas\n",
    "\n",
    "Now we can have some 10 m resolution DTM and DSM tiles and functions to read them in to RAM, we can to repeat the analysis from week 2 using an automatic python workflow.\n",
    "\n",
    "### Canopy height map\n",
    "\n",
    "The first step was to make a canopy height map which we can use to first calibrate a model between our lidar measurement and some ground biomass data, and then predict biomass across our site. To do this we need to read the DTM and DSM in to RAM as numpy arrays and then subtract the DTM from the DSM to make a new CHM array. A function to perform this, making use of our earlier ***readTiff()*** function, is provided below.\n",
    "\n",
    "**Add some lines to the end of this cell to plot the CHM you have just made.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# make sure you have run the readTiff() cell above to load that function.\n",
    "\n",
    "def heightAboveGround(dtmName,dsmName):\n",
    "    '''\n",
    "    A function to calculate height above ground from a DTM and DSM\n",
    "    Note that this assumes the two datasets are aligned and the same resolution\n",
    "    '''\n",
    "    \n",
    "    # open the DTM and DSM and read data\n",
    "    tData,tX0,tY0,tXres,tYres,nX,nY=readTiff(dtmName)  # data plus all metadata\n",
    "    sData,sX0,sY0,sXres,sYres,nX,nY=readTiff(dsmName)  # data plus all metadata\n",
    "    \n",
    "    # Subtract the two to get height, as two are aligned\n",
    "    hData=sData-tData\n",
    "    \n",
    "    # pass back to the calling function\n",
    "    return(hData,tX0,tY0,tXres,tYres,nX,nY)\n",
    "\n",
    "\n",
    "\n",
    "# call those functions\n",
    "\n",
    "dsmName='data/ALS/DSM_1m/NN89_1M_DSM_PHASE1.tif'\n",
    "dtmName='data/ALS/DTM_1m/NN89_1M_DTM_PHASE1.tif'\n",
    "\n",
    "\n",
    "# calculate the height array and metadata and read in to RAM\n",
    "hData,hX0,hY0,hXres,hYres,nX,nY=heightAboveGround(dtmName,dsmName)\n",
    "\n",
    "\n",
    "# Add some lines below here to make an image of the CHM, with the y axis correctly aligned."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have a function to calculate canopy height from a single DSM and DTM. We then want to process all of the data. As a first step, let us practice processing all of the data to height and making an image of each tile to see how a loop can be used to batch process.\n",
    "\n",
    "The code below will create a list of all of the DTM and DSM filenames. **Add a loop to step through these filenames and call the ***heightAboveGround()*** function defined above to calculate height, then plot an image of each tile.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['data/ALS/DTM/coarsen.NH70_10M_DTM_PHASE1.tif', 'data/ALS/DTM/coarsen.NH80_10M_DTM_PHASE1.tif', 'data/ALS/DTM/coarsen.NH81_10M_DTM_PHASE1.tif', 'data/ALS/DTM/coarsen.NH90_10M_DTM_PHASE1.tif', 'data/ALS/DTM/coarsen.NH91_10M_DTM_PHASE1.tif', 'data/ALS/DTM/coarsen.NH92_10M_DTM_PHASE1.tif', 'data/ALS/DTM/coarsen.NJ01_10M_DTM_PHASE1.tif', 'data/ALS/DTM/coarsen.NJ02_10M_DTM_PHASE1.tif', 'data/ALS/DTM/coarsen.NJ03_10M_DTM_PHASE1.tif', 'data/ALS/DTM/coarsen.NJ13_10M_DTM_PHASE1.tif', 'data/ALS/DTM/coarsen.NJ14_10M_DTM_PHASE1.tif', 'data/ALS/DTM/coarsen.NJ23_10M_DTM_PHASE1.tif', 'data/ALS/DTM/coarsen.NJ24_10M_DTM_PHASE1.tif', 'data/ALS/DTM/coarsen.NN69_10M_DTM_PHASE1.tif', 'data/ALS/DTM/coarsen.NN79_10M_DTM_PHASE1.tif', 'data/ALS/DTM/coarsen.NN89_10M_DTM_PHASE1.tif']\n"
     ]
    }
   ],
   "source": [
    "from glob import glob\n",
    "\n",
    "dtmDir='data/ALS/DTM'\n",
    "dsmDir='data/ALS/DSM'\n",
    "\n",
    "dtmList=sorted(glob(dtmDir+'/*.tif'))\n",
    "dsmList=sorted(glob(dsmDir+'/*.tif'))\n",
    "\n",
    "# check it has worked by printing the list to screen\n",
    "print(dtmList)\n",
    "\n",
    "# add a loop to calculate and plot a CHM for each tile"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Making the merged CHM\n",
    "\n",
    "In the cell above, you have looped over each tile and calculated the CHM. At the moment, you have plotted this on the screen, then the computer has deleted the data. We now want to save all of that data to disk so that we can use it in our later analysis. We have an irregular set of tiles that we want to intersect with our field data and then map biomass.\n",
    "\n",
    "There are lots of different ways to do this, but to simplify the intersection of the CHM raster with our field data, we shall create a set of CHM raster tiles and then merge all of these in to a single CHM raster tile. This is less RAM efficient than keeping all of the tiles separate, but avoids having to loop when intersecting the field data with the tiles.\n",
    "\n",
    "For the first step, we need to create a set of geotiffs for the CHM of each tile. The code below defines a function to write a single CHM to a new geotiff. **This is a great example of code reuse**. `writeCHMtiff` makes use of `readTiff`, and the eventual call to `writeCHMtiff` is one line, replacing about 50 lines of code that would be required to read two GeoTiffs (the DTM and DSM), subtract them, and then write to a new GeoTIFF. We can be even more efficient by creating a loop which creates and saves CHM rasters for **all** tiles -- saving hundreds of lines of code!!!\n",
    "\n",
    "**Add a loop to write a CHM geotiff for all tiles.** (*Hint: it will look similar to the loop you created above.*)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "def writeCHMtiff(dtmName,dsmName,outName):\n",
    "    '''A function to write a CHM geotiff, based on heightAboveGround()'''\n",
    "    \n",
    "    # open the DTM and DSM and read data\n",
    "    tData,tX0,tY0,tXres,tYres,nX,nY=readTiff(dtmName)  # data plus all metadata\n",
    "    sData,sX0,sY0,sXres,sYres,nX,nY=readTiff(dsmName)  # data plus all metadata\n",
    "    \n",
    "    # Subtract the two to get height, as two are aligned\n",
    "    hData=sData-tData\n",
    "\n",
    "    # get geolocation information from one file\n",
    "    inDS=gdal.Open(dtmName)\n",
    "    OutDriver = gdal.GetDriverByName('GTiff')\n",
    "    outDS = OutDriver.Create(outName,nX,nY,1,gdal.GDT_Float32 )\n",
    "    outDS.SetGeoTransform(inDS.GetGeoTransform())\n",
    "    outDS.SetProjection(inDS.GetProjection()) # copy projection info\n",
    "    outDS.GetRasterBand(1).WriteArray(hData)  # write image to the raster\n",
    "    outDS.FlushCache()                     # write to disk\n",
    "    outDS = None\n",
    "    \n",
    "    print('Written to',outName)\n",
    "    return\n",
    "\n",
    "\n",
    "# pick a tile\n",
    "dsmName='data/ALS/DSM/coarsen.NN89_10M_DSM_PHASE1.tif'\n",
    "dtmName='data/ALS/DTM/coarsen.NN89_10M_DTM_PHASE1.tif'\n",
    "    \n",
    "# make a new CHM raster geotiff\n",
    "outName='NN89_10M_CHM_PHASE1.tif'\n",
    "writeCHMtiff(dtmName,dsmName,outName)\n",
    "\n",
    "\n",
    "# modify this code to loop over all tiles to make a CHM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you have a set of CHM geotiff files. The code below will read a list of your CHM geotiffs and merge them all in to a single raster geotiff, which will then be stored in:\n",
    "\n",
    "    data/ALS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from glob import glob\n",
    "import os\n",
    "\n",
    "chmDir='.'\n",
    "chmList=glob(chmDir+'/*CHM_PHASE1.tif')\n",
    "\n",
    "outName='data/ALS/merged_CHM_10m.tif'\n",
    "\n",
    "\n",
    "# turn the above in to a string so that the command below works\n",
    "chmStr=\"\"\n",
    "for f in chmList:\n",
    "    chmStr=chmStr+\" \"+f\n",
    "\n",
    "\n",
    "command = \"gdal_merge.py -o \"+outName+\" -of gtiff \" + chmStr\n",
    "os.system(command)\n",
    "\n",
    "print('Written to',outName)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calibrating the biomass model\n",
    "\n",
    "Now you have a 10 m resolution CHM. You also have some field data in a csv file stored in:\n",
    "\n",
    "    data/ground/ground_data.csv\n",
    "    \n",
    "You can calibrate the biomass model in python. The ground data can be read in using the numpy.loadtxt() function to store the biomass (AGBD) and coordinates (x and y) in RAM.\n",
    "\n",
    "A few more optional arguments are used than in the week 9 practical -- a noteable difference is that the call returns not a single, multidimensional array, but 3 arrays corresponding to different columns in the spreadsheet. (We could alternatively use `pandas`.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "groundname='data/ground/ground_data.csv'\n",
    "\n",
    "# load data in to RAM\n",
    "agbd,x,y=np.loadtxt(groundname, usecols=(1,2,3), unpack=True, dtype=float,skiprows=1,delimiter=',')\n",
    "\n",
    "# Here agbd is in kg/ha. Let's use the more common Mg/ha\n",
    "agbd=agbd/1000\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we want to intersect this with our CHM raster to create an array with the biomass for each plot along with the mean CHM. We can then plot that data to see what type of mathematical model would be appropriate to predict biomass from lidar data.\n",
    "\n",
    "The script below will read in the raster layer and then return an array of CHM values corresponding to the coordinates you read in from the field data file.\n",
    "\n",
    "In this code, we are not actually *interpolating* the CHM to the point locations given by `x` and `y`, like we did in Week 9. Rather, for each location given by `x` and `y` we are finding a nearby grid point in the CHM, and using that value. The difference is generally small, as long as we are working with a closely spaced raster, but differences could be large with a coarse raster, or a raster in which values change rapidly between neighboring grid points.\n",
    "\n",
    "**Add some code to make a scatterplot of biomass against mean canopy height. What type of relationship do you see?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# read the CHM geotiff in to RAM\n",
    "chmName='data/ALS/merged_CHM_10m.tif'\n",
    "data,xOrigin,yOrigin,pixelWidth,pixelHeight,nX,nY=readTiff(chmName)\n",
    "\n",
    "\n",
    "# get a list of pixel coordinates for the ground plots\n",
    "xInds=np.array((x-xOrigin)/pixelWidth,dtype=int)\n",
    "yInds=np.array((y-yOrigin)/pixelHeight,dtype=int)\n",
    "\n",
    "# extract the values to make a new array of mean canopy height\n",
    "meanCH=data[xInds,yInds]\n",
    "\n",
    "# add some code to make a scatterplot of biomass (agbd) against mean canopy height (meanCH)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you have viewed the data and chosen an appopriate mathamatical model (was it linear, loigarithmic etc.?), you can perform linear regression within python to get the equation for your biomass model. \n",
    "\n",
    "The code below shows you can example of fitting a line of best fit to two random datasets, and calculates the correlation (Pearson's $r$ value). We saw `linregress` last week, but not `pearsonr`. However, Pearson's $r$ is the correlation between $x$ and $y$, which is also returned by `linregress`. The difference is that `pearsonr` returns its value *and significance*. You may recall $R^2$, the *coefficient of determination*, from learning about regression in high school. Despite the upper/lower case difference, in linear regression, $R^2$ is the square of Pearson's $r$.\n",
    "\n",
    "**Modify that code to find the equation needed to predict biomass from mean canopy height and report the strength of the correlation.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# packages to do linear regression and Pearson's correlation\n",
    "from scipy.stats import linregress, pearsonr\n",
    "\n",
    "# make some fake data\n",
    "x=np.random.random((20))*100\n",
    "y=x*3+20+40*(np.random.random((20))-0.5)\n",
    "\n",
    "# fit a line of best fit. See manual for scipy.stats.linregress for more details\n",
    "m, c, r, _, _ = linregress(x,y)\n",
    "\n",
    "print(\"Line of best fit parameters are m\",m,\"c\",c)\n",
    "\n",
    "# find the correlation. See manual for scipy.stats.pearsonr\n",
    "print('Correlation is',pearsonr(x,y))\n",
    "\n",
    "# the first number returned by pearsonr is the correlation, 99.28% for this fake data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Making the biomass map\n",
    "\n",
    "\n",
    "Now you have the parameters needed to parameterise a biomass model. **In the cell below, write some code to read in the CHM again (you can reuse the code above) and make a new array of biomass using the model you have just set.** If this was a linear model, the biomass will be:\n",
    "\n",
    "$biomass = m \\times meanCH + c $\n",
    "\n",
    "You can use ***plt.imshow()*** to print it as an image on your screen, or modify the CHM writing function above to write the biomass raster layer as a new geotiff. Use whichever approach you prefer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Using numpy operations, what is the mean biomass for this valley? What is the maximum biomass for any 10 m pixel?\n",
    "\n",
    "Try using the pyplot histogram function to look at the distribution of biomass values in this area.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A comment on automation\n",
    "\n",
    "\n",
    "As you have seen above, making a biomass map with Python requires quite a few functions, some using some unintuitive aspects of the language. However, you now have a workflow to process a set of geotiff tiles through to biomass. You can rerun the whole workflow by pressing the double arrow button at the top of the screen.\n",
    "\n",
    "You can do this for any data by downloading new tiles and putting them in the directories above. You can make adjustments to the workflow instructions (eg. change the resolution you do the calculation at, or add in an extra bit of processing) and then reurn everything with a single button."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visually assessing the results\n",
    "\n",
    "You have now produced a biomass map entirely in Python. Whilst python does allow you to make plots of maps and to calculate summary statistics, it can be clunky to assess one map against other. For that task, GIS allows you to more interactively work with the data.\n",
    "\n",
    "Download the biomass map you have made (go to the tab to the left of this, navigate to the map and tick and download it). Load this in to QGIS (either on your own computer or Apps.ed) and overlay it with a Bing aerial image.\n",
    "\n",
    "Can you see any obvious errors in the biomass map? Can you see a way to limit your analysis to a particular estate or a particular land type?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Further work\n",
    "\n",
    "There are lots of different options to colour your maps. Some can make very attractive and effective visualisations. You can look through the notes below to see some examples.\n",
    "\n",
    "Keyword argument `cmap` in `plt.imshow()` selects the colour map. The default is called `'viridis'`, which is a perceptually uniform sequential colour map (one which gives uniform steps in perceived brightness for uniform steps in data values), but `cmap='gray'` gives a greyscale which works perfectly well for this example. Experiment with the matplotlib colour maps listed at  http://matplotlib.org/users/colormaps.html to find ones that you think are particularly good or bad for showing this data. Until recently, the matplotlib default was the rainbow scale called `'jet'`, and you will still see many examples of plots with rainbow scales in research literature, even though they are now considered to be a poor choice.\n",
    "\n",
    "The axes are labelled in grid coordinates (numbers of points) rather than geographical coordinates. The correct extents of the axes (800 × 2 = 1600 m) can be set with keyword argument `extent=(0,1600,0,1600)`. Use `plt.xticks()` and `plt.yticks()` if you want to control where the axis tick marks are placed and how they are labelled.\n",
    "\n",
    "To have a colour bar beside the map, add `cbar = plt.colorbar()`. \n",
    "\n",
    "Having assigned a name to the colour bar object, there are several methods that can modify it. Use the `set_label()` method with a text string to label the colour bar, e.g. `cbar.set_label('Elevation (m)')`, and use `cbar.set_ticks(list)` with a list of numbers to set the locations of tick marks on the bar.\n",
    "\n",
    "Here we have set an elevation parameter equal to the variable z, and the resolutions equal to dx and dy. That can be done by reading a DTM geotiff as z, rather than data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set a filename\n",
    "filename='data/ALS/DSM_1m/NN89_1M_DSM_PHASE1.tif'\n",
    "\n",
    "# use the function above to read the data in to RAM\n",
    "z,xOrigin,yOrigin,dz,dy,nX,nY=readTiff(filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradients and slopes\n",
    "\n",
    "The components of the gradient vector for array `z` can be estimated by central differences using \n",
    "```\n",
    "dzdy, dzdx = np.gradient(z, dx)\n",
    "```\n",
    "where `dx` is the grid spacing (2 m for the Arthur’s Seat DEM). The slope is calculated and converted to slope angle in degrees by \n",
    "```\n",
    "slope = np.sqrt(dzdx**2 + dzdy**2)\n",
    "angle = np.rad2deg(np.arctan(slope))\n",
    "```\n",
    "\n",
    "Make a figure with three images of `dzdx`, `dzdy` and `angle` side by side. In this layout, you might think that horizontal colour bars under the plots will fit better; use `plt.colorbar(orientation='horizontal')` to do this. If you want to shrink or expand the bar to a fraction of its default size, use keyword argument `shrink`.\n",
    "\n",
    "Because `dzdx` and `dzdy` have both positive and negative values, use a diverging scale from the list of colour maps. Ensure that zero is in the middle of the scale by setting symmetrical negative and positive values for keyword arguments `vmin` and `vmax` in `plt.imshow()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "\n",
    "Find the maximum slope angle. The places with the steepest slopes should be clearly visible in the slope angle image; check that they are where you expect them to be.\n",
    "\n",
    "\n",
    "## Topographic shading\n",
    "\n",
    "A shaded surface map can be created with the lines\n",
    "```\n",
    "from matplotlib.colors import LightSource\n",
    "ls = LightSource(azdeg=315, altdeg=45)\n",
    "shadsurf = ls.shade(z, cmap=plt.cm.gray)\n",
    "plt.imshow(shadsurf, origin='lower')\n",
    "```\n",
    "for an array of elevations `z` that has already been loaded. `azdeg=315` and `altdeg=45` are the default azimuth and elevation angles in degrees for the light source.\n",
    "\n",
    "Make a shaded surface plot of Arthur’s Seat and experiment with different azimuth angles (0º to 360º) and elevation angles (0º to 90º) for the light source to see how the appearance of the plot changes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Contours can be added to a shaded map simply by putting `plt.contour()` after `plt.imshow()`. Pleasing effects can be produced by shading a coloured image map, but this requires making the shading semi-transparent; add keyword argument `alpha` with a value between 0 (fully transparent) and 1 (opaque)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}