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:
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.
%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')
# 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
The header info is here. You may want to keep it open in another tab for easy reference.
#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 = 'http://ropitz.github.io/digitalantiquity/data/red.tif'
GREEN = 'http://ropitz.github.io/digitalantiquity/data/green.tif'
BLUE = 'http://ropitz.github.io/digitalantiquity/data/blue.tif'
IR = 'http://ropitz.github.io/digitalantiquity/data/IR.tif'
ALL = 'http://ropitz.github.io/digitalantiquity/data/carnuntum.tif'
# You can open and plot each raster.
nir_ds = gdal.Open(IR)
nir_band = nir_ds.GetRasterBand(1)
nir = nir_band.ReadAsArray()
fig = plt.figure(figsize=(15,15))
plt.imshow(nir)
plt.colorbar()
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.
#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()
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.
# Read out band 44
band44 = dataset.read((44))
band44
# Read out band 16
band16 = dataset.read((16))
band16
# 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()
#Let's try it with some other bands
band60 = dataset.read((60))
band60
band30 = dataset.read((30))
band30
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?
#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))
b30 = b30.astype(np.float64)
b35 = b35.astype(np.float64)
b39 = b39.astype(np.float64)
b42 = b42.astype(np.float64)
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.
REIPCL = (b39/b35)-1
fig = plt.figure(figsize=(15,15))
plt.imshow(REIPCL, cmap='RdYlGn')
plt.colorbar()
That shows rather more!
Hopefully you have:
We'll be talking more about remote sensing methods in archaeology throughout the course.