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 [ ]:
%matplotlib inline

from osgeo import gdal
import gdal
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 [ ]:
# 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

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 [ ]:
#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 = rasterio.open('http://ropitz.github.io/digitalantiquity/data/red.tif')
GREEN = rasterio.open('http://ropitz.github.io/digitalantiquity/data/green.tif')
BLUE = rasterio.open('http://ropitz.github.io/digitalantiquity/data/blue.tif')
IR = rasterio.open('http://ropitz.github.io/digitalantiquity/data/IR.tif')
ALL = rasterio.open('http://ropitz.github.io/digitalantiquity/data/carnuntum.tif')

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 [ ]:
#Let's start by calculating NDVI

# Read out band 16 and 19
band16 = dataset.read((16))
band19 = dataset.read((19))
band19 = band19.astype(np.float64)
band16 = band16.astype(np.float64)
def ratio1(band19, band16):
    """Calculate custom ratio"""
    return (band19 - band16) / (band19 + band16)

fig = plt.figure(figsize=(15,15))
ratio1 = ratio1(band19, band16)
plt.imshow(ratio1, cmap='RdYlGn')
plt.colorbar()

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 [ ]:
# Read out and plot band 19
band19 = dataset.read((19))
plt.imshow(band19, cmap='RdYlGn')
plt.colorbar()
In [ ]:
# Read out and plot band 16
band16 = dataset.read((16))
plt.imshow(band16, cmap='RdYlGn')
plt.colorbar()
In [ ]:
# Read out band 44
band44 = dataset.read((44))
band44
In [ ]:
# Read out band 16
band16 = dataset.read((16))
band16
In [ ]:
# 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()
In [ ]:
#Let's try it with some other bands
band60 = dataset.read((60))
band60
In [ ]:
band30 = dataset.read((30))
band30
In [ ]:
def ratio2(band60, band30):
    """Calculate custom ratio"""
    return (band60 - band30) / (band60 + band30)


band60 = band60.astype(np.float64)
band30 = band30.astype(np.float64)
fig = plt.figure(figsize=(15,15))
ratio2 = ratio2(band60, band30)
plt.imshow(ratio2, cmap='RdYlGn')
plt.colorbar()

Compare the results. Are some band combinations more effective at showing cropmarks than others?

In [ ]:
#REIP is another popular index for highlighting cropmarks. Calculate REIP
b30 = dataset.read((30))
b35 = dataset.read((35))
b39 = dataset.read((39))
b42 = dataset.read((42))
b43 = dataset.read((43))
b34 = dataset.read((34))
b30 = b30.astype(np.float64)
b35 = b35.astype(np.float64)
b39 = b39.astype(np.float64)
b42 = b42.astype(np.float64)
b43 = b43.astype(np.float64)
b34 = b34.astype(np.float64)
In [ ]:
REIP = 700 + 40 * ((((b30 + b42)/2)-b35) / (b39 - b35))



REIPplot = 700 + (40 * ((((b30 + b43)/2)-b34) / (b39 - b34)))
fig = plt.figure(figsize=(15,15))
plt.imshow(REIPplot, cmap='RdYlGn')
plt.colorbar()

I don't think that index shows very much in our data. Let's look up Another REIP Index and try that instead.

In [ ]:
REIPCL = (b39/b35)-1
fig = plt.figure(figsize=(15,15))
plt.imshow(REIPCL, cmap='RdYlGn')
plt.colorbar()

That shows rather more!

That concludes this tutorial.

Hopefully you have:

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

We'll be talking more about remote sensing methods in archaeology throughout the course.