batch_process.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PPRS: Week 11\n",
    "# Spatial data analysis and batch processing in Python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run this to display graphs nicely in this notebook\n",
    "\n",
    "%matplotlib inline\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "This week we will use python 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. This follows on from the introduction to lidar session and is meant as an introduction to OOSA. We will use 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",
    "We are using the DTM and DSM provided by the [Scottish government](https://remotesensingdata.gov.scot/data#/list). These were collected as point clouds and then processed to a DTM and DSM in geotiff format. After cloning these will be in the ***data*** directory. The raw data comes in 1 m resolution but the ALS point density was a little too low to support that, so we have coarsened it first to 3 m, taking the mean elevation for the DTM and the maximum elevation for the DSM, and then to 10 m taking the mean for both. Due to disk space limits on noteable, the data has already been coarsened for you.\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 rasterio package to open our data. This is built on top of GDAL. Below you will see code which uses this module and its functions to open and read GeoTIFF files. Note that *there are many many many Python modules out there, some with overlapping functionality.* Another popular geospatial package is `rasterio`.\n",
    "\n",
    "The cell below provides a function to open a single raster file. This is an example of a **User Defined Function**. 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 input ***arguments*** (a filename in this case) and return value list (the data along with the metadata in this case)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in the necessary packages\n",
    "from osgeo import gdal,osr\n",
    "import numpy as np\n",
    "from sys import exit\n",
    "import rasterio\n",
    "\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=rasterio.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",
    "    print(\"The projection is\",ds.crs)\n",
    "\n",
    "    # get the raster size\n",
    "    nX=ds.width             # number of pixels in x direction\n",
    "    nY=ds.height             # number of pixels in y direction\n",
    "    \n",
    "    # geolocation tiepoint\n",
    "    xOrigin=ds.transform[2]       # coordinate of x corner\n",
    "    yOrigin=ds.transform[5]       # coordinate of y corner\n",
    "    pixelWidth=ds.transform[0]    # resolution in x direction\n",
    "    pixelHeight=ds.transform[4]   # resolution in y direction\n",
    "\n",
    "    # read data. Returns as a 2D numpy array\n",
    "    data=ds.read(1)\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",
    "# now call the function we have defined\n",
    "\n",
    "# set a filename\n",
    "filename='data/ALS/DSM/coarsen.NN89_10M_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, using the np.min() and np.max() functions.**"
   ]
  },
  {
   "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\n"
   ]
  },
  {
   "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. This is part of the same package we used to plot data in the foundations notebook, but is set up for making 2D images.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# plot the 2D raster layer  \n",
    "plt.imshow(data)   # plot the 2D image\n",
    "plt.show()         # print to screen \n",
    "\n"
   ]
  },
  {
   "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/coarsen.NN89_10M_DSM_PHASE1.tif'\n",
    "dtmName='data/ALS/DTM/coarsen.NN89_10M_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"
   ]
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "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",
    "\n",
    "# add a loop to calculate and plot a CHM for each tile\n",
    "\n",
    "# prepare a plot space\n",
    "nFigs=len(dtmList)\n",
    "\n",
    "\n",
    "\n",
    "for i in range(0,len(dtmList)):\n",
    "    dtmName=dtmList[i]\n",
    "    dsmName=dsmList[i]\n",
    "    \n",
    "    # print the filenames to screen\n",
    "    print(dtmName,dsmName)\n",
    "    \n",
    "    # add some lines to calculate the height and plot an image\n"
   ]
  },
  {
   "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 good 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",
    "The code below loops over the tiles to make a set of CHMs and write new geotiffs to disk.\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",
    "\n",
    "# get a list of tilenames\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",
    "\n",
    "# modify this code to loop over all tiles to make a CHM\n",
    "for i in range(0,len(dtmList)):\n",
    "    dtmName=dtmList[i]\n",
    "    dsmName=dsmList[i]\n",
    "   \n",
    "    tileLabel=dtmName.split('/')[-1].split('_')[0]\n",
    "    outName=tileLabel+'_10M_CHM_PHASE1.tif'\n",
    "\n",
    "    writeCHMtiff(dtmName,dsmName,outName)\n",
    "\n"
   ]
  },
  {
   "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",
    "\n",
    "print(agbd)\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)\n"
   ]
  },
  {
   "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\n",
    "\n",
    "\n",
    "# repeat with the real 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": [
    "# read the CHM geotiff in to RAM\n",
    "\n",
    "\n",
    "# print to screen \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If there were any large outliers in the CHM, either large positive values due to noise points above ground. or negative due to holes in the data, these can stretch the contrast and make the image above hard to see. These can be removed by either filtering the data (ideally filtering the raw point cloud before making the DSM and DTM), or by adjusting the colour bar on the image. The cell below shows how to filter out all values above 10, for example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filtValue=300\n",
    "filtBiomass=biomass\n",
    "filtBiomass[biomass>filtValue]=filtValue"
   ]
  },
  {
   "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?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# find the mean and max\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Writing the biomass map to file for reuse\n",
    "\n",
    "We have displayed the biomass map to creen above. Often we will want to save it for later analysis, either in python or in GIS. We casn do this using GDAL. An example function is given below.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def writeTiff(inName,outName,data):\n",
    "    '''A function to read a tiff and then write a new tiff with the same geolocation but different properties'''\n",
    "\n",
    "    # get geolocation information from one file (the CHm file in this case)\n",
    "    inDS=gdal.Open(inName)\n",
    "    \n",
    "    # write the output to a file\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(data)  # 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",
    "\n",
    "# set a filename\n",
    "outName='biomass.tif'\n",
    "\n",
    "# write the data, matching the geolocation of the CHM\n",
    "writeTiff(chmName,outName,biomass)"
   ]
  },
  {
   "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?"
   ]
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}