Learning about spatial data and maps for archaeology (and other things)

Spatial Thinking and Skills Exercise 1 for Theory and Practice

Made by Rachel Opitz, Archaeology, University of Glasgow

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.

The aim of this exercise is for you to:

  • learn to work with spectral images, including extracting individual bands and calculating common band ratios
  • learn which bands and ratios are frequently used to spot vegetation stress aka cropmarks
  • start thinking about the use of remote sensing in archaeological survey and landscape studies

You'll do this using data collected over Carnuntum, a Roman site in Austira, made available by the LBI.

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.

Let's get started... Hit 'Ctrl'+'Enter' to run the code in any cell in the page.

In [60]:
%matplotlib inline

from osgeo import gdal
import gdal; gdal.UseExceptions()
import rasterio as rasterio
from matplotlib import pyplot
import numpy as np
import pygeoprocessing
from skimage import data, io, filters
import matplotlib.pyplot as plt
import cv2
import warnings
warnings.filterwarnings('ignore')
In [7]:
# First we want to know a little about our dataset. 
# How many bands does our image have? Let's count them.

dataset = rasterio.open('http://ropitz.github.io/digitalantiquity/data/carnuntum.tif')

dataset.count
Out[7]:
65

65 bands is a lot! Now we want to know which wavelength is stored in each band.

Get the information about which wavelength is in each band from the files header.

The header info is here. You may want to keep it open in another tab for easy reference.

In [56]:
#I pre-extracted the most commonly used and important bands into their own raster images. 
# Rasters are just fancy images.
# We'll assign RGB + IR to their own variables so we can work with them.
RED = 'data/red.tif'
GREEN = 'data/green.tif'
BLUE = 'data/blue.tif'
IR = 'data/IR.tif'
ALL = 'data/carnuntum.tif'
In [49]:
nir_ds = gdal.Open('data/IR.tif')

nir_ds.GetRasterBand(1)
Out[49]:
<osgeo.gdal.Band; proxy of <Swig Object of type 'GDALRasterBandShadow *' at 0x7fba33eb2f30> >
In [65]:
 # You can open and plot each raster.
nir_ds = gdal.Open(IR)
nir_band = nir_ds.GetRasterBand(1)
nir = nir_band.ReadAsArray()
plt.imshow(nir)
plt.colorbar()
Out[65]:
<matplotlib.colorbar.Colorbar at 0x7fb9fc063198>
In [63]:
red_ds = gdal.Open(RED)
red_band = red_ds.GetRasterBand(1)
red = red_band.ReadAsArray()
plt.imshow(red)
plt.colorbar()
Out[63]:
<matplotlib.colorbar.Colorbar at 0x7fb9fc1a3ba8>

Calculating Spectral Indices

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. Spectral Indices that enhance the appearance of cropmarks (and other indices) can be found here.

In [66]:
#Let's start by calculating NDVI
def ndvi(red, nir):
    """Calculate NDVI."""
    return (nir - red) / (nir + red)


red = red.astype(np.float64)
nir = nir.astype(np.float64)
fig = plt.figure(figsize=(15,15))
ndvi = ndvi(red, nir)
plt.imshow(ndvi, cmap='RdYlGn')
plt.colorbar()
Out[66]:
<matplotlib.colorbar.Colorbar at 0x7fb9f5bb3ac8>

Compare the NDVI image to the images of individual bands. Can you see more features in one than another? Different features?

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.

In [67]:
# Read out band 44
band44 = dataset.read((44))
band44
Out[67]:
array([[36315, 36161, 36731, ..., 37971, 38070, 38012],
       [36194, 36181, 36967, ..., 37979, 37873, 37758],
       [36475, 36557, 36298, ..., 37850, 37909, 37946],
       ...,
       [37321, 37244, 37151, ..., 37739, 37549, 37460],
       [37136, 37299, 37395, ..., 37545, 37722, 37635],
       [37056, 37071, 37259, ..., 37462, 37652, 37615]], dtype=uint16)
In [68]:
# Read out band 16
band16 = dataset.read((16))
band16
Out[68]:
array([[34206, 34129, 34047, ..., 33682, 33704, 33784],
       [34222, 34178, 34037, ..., 33708, 33727, 33744],
       [34038, 34011, 34337, ..., 33659, 33732, 33748],
       ...,
       [33852, 33844, 33793, ..., 33299, 33326, 33330],
       [33810, 33839, 33848, ..., 33293, 33317, 33321],
       [33835, 33819, 33834, ..., 33302, 33308, 33313]], dtype=uint16)
In [69]:
# Now we make a custom ratio, like NDVI but with different bands.
def ratio1(band44, band16):
    """Calculate custom ratio"""
    return (band44 - band16) / (band44 + band16)


band44 = band44.astype(np.float64)
band16 = band16.astype(np.float64)
fig = plt.figure(figsize=(15,15))
ratio1 = ratio1(band44, band16)
plt.imshow(ratio1, cmap='RdYlGn')
plt.colorbar()
Out[69]:
<matplotlib.colorbar.Colorbar at 0x7fb9f52ceac8>