{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Learning about spatial data and maps for archaeology (and other things)\n",
    "\n",
    "### Spatial Thinking and Skills Exercise 1 for Theory and Practice\n",
    "\n",
    "#### Made by Rachel Opitz, Archaeology, University of Glasgow\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remote sensing data, from satellite images, to lidar, to multispectral data, is used in archaeology to identify sites and features throughout the landscape. Recently, hyperspectral data - which is like multispectral data but with even more and narrower bands - has been used to spot cropmarks before they appear in the visible spectrum and are apparent to the naked eye. \n",
    "\n",
    "The aim of this exercise is for you to:\n",
    "* learn to work with spectral images, including extracting individual bands and calculating common band ratios\n",
    "* learn which bands and ratios are frequently used to spot vegetation stress aka cropmarks\n",
    "* start thinking about the use of remote sensing in archaeological survey and landscape studies\n",
    "\n",
    "You'll do this using data collected over Carnuntum, a Roman site in Austira, made available by the LBI. \n",
    "\n",
    "As you may recall from Archaeology of Scotland, to start working with spatial data and imagery, you need to put together your toolkit. You're currently working inside something called a jupyter notebook. It's a place to keep notes, pictures, code and maps together. You can add tools and data into your jupyter notebook and then use them to ask spatial questions and make maps and visualisations that help answer those questions. \n",
    "\n",
    "\n",
    "### Let's get started... Hit 'Ctrl'+'Enter' to run the code in any cell in the page."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from osgeo import gdal\n",
    "import gdal\n",
    "import rasterio as rasterio\n",
    "from matplotlib import pyplot\n",
    "import numpy as np\n",
    "import pygeoprocessing\n",
    "from skimage import data, io, filters\n",
    "import matplotlib.pyplot as plt\n",
    "import cv2\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First we want to know a little about our dataset. \n",
    "# How many bands does our image have? Let's count them.\n",
    "\n",
    "dataset = rasterio.open('http://ropitz.github.io/digitalantiquity/data/carnuntum.tif')\n",
    "\n",
    "dataset.count"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 65 bands is a lot! Now we want to know which wavelength is stored in each band.\n",
    "### Get the information about which wavelength is in each band from the files header.\n",
    "\n",
    "[The header info is here.](http://ropitz.github.io/digitalantiquity/data/carnuntum.txt) You may want to keep it open in another tab for easy reference.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#I pre-extracted the most commonly used and important bands into their own raster images. \n",
    "# Rasters are just fancy images.\n",
    "# We'll assign RGB + IR to their own variables so we can work with them.\n",
    "RED = rasterio.open('http://ropitz.github.io/digitalantiquity/data/red.tif')\n",
    "GREEN = rasterio.open('http://ropitz.github.io/digitalantiquity/data/green.tif')\n",
    "BLUE = rasterio.open('http://ropitz.github.io/digitalantiquity/data/blue.tif')\n",
    "IR = rasterio.open('http://ropitz.github.io/digitalantiquity/data/IR.tif')\n",
    "ALL = rasterio.open('http://ropitz.github.io/digitalantiquity/data/carnuntum.tif')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating Spectral Indices\n",
    "\n",
    "To make cropmarks more visible we calculate different indices. Indices are combinations of spectral bands that highlight certain properties of the vegetation or soil. NDVI is one of the most common indices, and highlights the amount of clorophyll, which roughly corresponds to how healthy or stressed a plant is. \n",
    "[Spectral Indices that enhance the appearance of cropmarks (and other indices) can be found here.](https://www.indexdatabase.de/db/i-single.php?id=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Let's start by calculating NDVI\n",
    "\n",
     "# Read out band 16 and 19\n",
    "band16 = dataset.read((16))\n",
    "band19 = dataset.read((19))\n",
    "band19 = band19.astype(np.float64)\n",
    "band16 = band16.astype(np.float64)\n",
    "def ratio1(band19, band16):\n",
    "    \"\"\"Calculate custom ratio\"\"\"\n",
    "    return (band19 - band16) / (band19 + band16)\n",
    "\n",
    "fig = plt.figure(figsize=(15,15))\n",
    "ratio1 = ratio1(band19, band16)\n",
    "plt.imshow(ratio1, cmap='RdYlGn')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare the NDVI image to the images of individual bands. Can you see more features in one than another? Different features? \n",
    "\n",
    "We can also create our own custom ratios by reading individual bands out of our big Carnuntum dataset and combining them in different ways. This is a good way to experiment and explore our data."
   ]
  },  
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read out and plot band 19\n",
    "band19 = dataset.read((19))\n",
    "plt.imshow(band19, cmap='RdYlGn')\n",
    "plt.colorbar()"
   ]
  }, 
    {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read out and plot band 16\n",
    "band16 = dataset.read((16))\n",
    "plt.imshow(band16, cmap='RdYlGn')\n",
    "plt.colorbar()"
   ]
  }, 
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read out band 44\n",
    "band44 = dataset.read((44))\n",
    "band44"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read out band 16\n",
    "band16 = dataset.read((16))\n",
    "band16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we make a custom ratio, like NDVI but with different bands.\n",
    "def ratio1(band44, band16):\n",
    "    \"\"\"Calculate custom ratio\"\"\"\n",
    "    return (band44 - band16) / (band44 + band16)\n",
    "\n",
    "\n",
    "band44 = band44.astype(np.float64)\n",
    "band16 = band16.astype(np.float64)\n",
    "fig = plt.figure(figsize=(15,15))\n",
    "ratio1 = ratio1(band44, band16)\n",
    "plt.imshow(ratio1, cmap='RdYlGn')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Let's try it with some other bands\n",
    "band60 = dataset.read((60))\n",
    "band60"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "band30 = dataset.read((30))\n",
    "band30"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ratio2(band60, band30):\n",
    "    \"\"\"Calculate custom ratio\"\"\"\n",
    "    return (band60 - band30) / (band60 + band30)\n",
    "\n",
    "\n",
    "band60 = band60.astype(np.float64)\n",
    "band30 = band30.astype(np.float64)\n",
    "fig = plt.figure(figsize=(15,15))\n",
    "ratio2 = ratio2(band60, band30)\n",
    "plt.imshow(ratio2, cmap='RdYlGn')\n",
    "plt.colorbar()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare the results. Are some band combinations more effective at showing cropmarks than others?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#REIP is another popular index for highlighting cropmarks. Calculate REIP\n",
    "b30 = dataset.read((30))\n",
    "b35 = dataset.read((35))\n",
    "b39 = dataset.read((39))\n",
    "b42 = dataset.read((42))\n",
    "b43 = dataset.read((43))\n",   
    "b34 = dataset.read((34))\n",  
    "b30 = b30.astype(np.float64)\n",
    "b35 = b35.astype(np.float64)\n",
    "b39 = b39.astype(np.float64)\n",
    "b42 = b42.astype(np.float64)\n",
    "b43 = b43.astype(np.float64)\n",
    "b34 = b34.astype(np.float64)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "REIP = 700 + 40 * ((((b30 + b42)/2)-b35) / (b39 - b35))\n",
    "\n",
    "\n",
    "\n",
    "REIPplot = 700 + (40 * ((((b30 + b43)/2)-b34) / (b39 - b34)))\n",
    "fig = plt.figure(figsize=(15,15))\n",
    "plt.imshow(REIPplot, cmap='RdYlGn')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I don't think that index shows very much in our data. Let's look up [Another REIP Index](https://www.indexdatabase.de/db/i-single.php?id=190) and try that instead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "REIPCL = (b39/b35)-1\n",
    "fig = plt.figure(figsize=(15,15))\n",
    "plt.imshow(REIPCL, cmap='RdYlGn')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That shows rather more!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### That concludes this tutorial.\n",
    "\n",
    "Hopefully you have:\n",
    "* learned to work with spectral images, including extracting individual bands and calculating common band ratios\n",
    "* learned which bands and ratios are frequently used to spot vegetation stress aka cropmarks\n",
    "* started thinking about the use of remote sensing in archaeological survey and landscape studies\n",
    "\n",
    "We'll be talking more about remote sensing methods in archaeology throughout the course."
   ]
  }
 ],
 "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.5"
  },
  "widgets": {
   "state": {
    "3c0b9dc3198943dc810a7ba3a4c04249": {
     "views": [
      {
       "cell_index": 18
      }
     ]
    },
    "614fe957a514497bbb4c2b7059c09274": {
     "views": [
      {
       "cell_index": 27
      }
     ]
    },
    "908bc6883c354c2699a617c27ace9de7": {
     "views": [
      {
       "cell_index": 19
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}