Introduction

Because my project topic draws from sampling strategies in currently occupied medieval villages, I wanted to spend some time considering area excavations of medieval settlements. The interpretation of test pit data for currently occupied rural settlements is based on an analysis of the pottery from each test pit. The goal of this part of my project is to use R to calculate the average distribution of pottery within an open area excavation of a medieval rural settlement. I will create a visualization of the change in pottery distribution with increasing distance from the settlement. This visualization will be useful in considering the role of distance from the medieval street frontages in sampling currently occupied rural settlements.

The area excavation that I selected for comparison was at the London Heathrow Airport in advance of the construction of Terminal 5. I downloaded the complete spatial data archive for the excavated features and the pottery from the Archaeology Data Service

Install Packages

library(maptools)
library(rgdal)
library(rgeos)
library(leaflet)

Load Shapefiles

I downloaded the shapefiles directly from the Archaeology Data Service archive. I filtered through the finds database to create a shapefile with only the medieval pottery. I also used the medieval features shapefile and copied medieval structural features (postholes and beamslots) as a shapefile. The projection information was not provided with the original files, but I worked out that the initial downloads were in the British National Grid system. I used QGIS to transform these to WGS84: EPSG4326 when I saved the shapefiles to my desktop

heathrow <- readShapeSpatial("~/Grad Year 1/Data Mapping/Final Project/HeathrowMedStructuralwgs.shp")
heathrow.pottery <- readShapeSpatial("~/Grad Year 1/Data Mapping/Final Project/HeathrowAllPotterywgs3.shp")
heathrow.area <- readShapeSpatial("~/Grad Year 1/Data Mapping/Final Project/MedievalFeatureswgs.shp")

I could not get the readOGR commmand to load geoJSON from github to R, so my own code loads the shapefiles from my computer. I have put the shapefiles on github:

heathrow <- structural heathrow.pottery <- pottery heathrow.area <- all medieval features

Transform Coordinate System to Facilitate Area Calculations

I assigned a projection to the shapefiles and then used the spTransform function to reproject the shapefiles onto a coordinate system that uses meters. This will make meters the unit of spatial calculations for the remaining commands.

proj4string(heathrow) <- CRS("+proj=longlat +datum=WGS84")
proj4string(heathrow.pottery) <- CRS("+proj=longlat +datum=WGS84")
proj4string(heathrow.area) <- CRS("+proj=longlat +datum=WGS84")

spTransform(heathrow, CRS("+init=epsg:3035")) -> heathrow.t
spTransform(heathrow.pottery, CRS("+init=epsg:3035")) -> heathrow.pottery.t
spTransform(heathrow.area, CRS("+init=epsg:3035")) -> heathrow.area.t

gBuffer Calculations

Because the goal is to create a visualization that addresses the potential impact of increasing distance from settlement structures on pottery distribution, the first step will be to calculate a series of buffer areas from the structural features.

gBuffer(heathrow.t, width=250)-> heathrow.t.250
gBuffer(heathrow.t, width=150) -> heathrow.t.150
gBuffer(heathrow.t, width=100) -> heathrow.t.100
gBuffer(heathrow.t, width=75) -> heathrow.t.75
gBuffer(heathrow.t, width=50) -> heathrow.t.50
gBuffer(heathrow.t, width=25) -> heathrow.t.25

gDifference Calculations

The above buffer calculations result in a series of overlapping buffer areas. The area of the 25 meter buffer zone, therefore, is included as part of the area of the 250 meter buffer zone. If these are left in the current format, an attempt to overlay pottery in the 250 meter buffer zone will include all pottery within the 250 meter buffer zone rather than exclusively the pottery between the end of the 150 meter buffer zone and the end of the 250 meter buffer zone.

By using the gDifference function, the area of the 250 meter zone will only include the band of area between the end of the 150 meter zone and the end of the 250 meter zone.

gDifference(heathrow.t.250, heathrow.t.150) -> heathrow.diff.250
gDifference(heathrow.t.150, heathrow.t.100) -> heathrow.diff.150
gDifference(heathrow.t.100, heathrow.t.75) -> heathrow.diff.100
gDifference(heathrow.t.75, heathrow.t.50) -> heathrow.diff.75
gDifference(heathrow.t.50, heathrow.t.25) -> heathrow.diff.50
heathrow.t.25 -> heathrow.diff.25

Test Plot

I’ve included a section here with a simple plot of the above code to show what has been done thus far.

plot(heathrow.diff.250, col="purple")
plot(heathrow.diff.150, col="blue", add=TRUE)
plot(heathrow.diff.100, col="green", add=TRUE)
plot(heathrow.diff.75, col="yellow", add=TRUE)
plot(heathrow.diff.50, col="orange", add=TRUE)
plot(heathrow.diff.25, col="red", add=TRUE)
plot(heathrow.pottery.t, add=TRUE)
plot(heathrow.area.t, add=TRUE)

gIntersection

I eventually plan to calculate the pottery weight density per square meter. This calculation will make the furthest distances look especially skewed to be smaller if I use the current calculations. In order to have an accurate understanding of pottery distribution, it is important to only divide the total pottery sherds by the area excavated rather than by the total area that happens to fall within a given distance of the structural features. The gIntersection function will return only areas excavated within each buffer zone.

gIntersection(heathrow.diff.250, heathrow.area.t, byid=T) -> heathrow.area.250
gIntersection(heathrow.diff.150, heathrow.area.t, byid=T) -> heathrow.area.150
gIntersection(heathrow.diff.100, heathrow.area.t, byid=T) -> heathrow.area.100
gIntersection(heathrow.diff.75, heathrow.area.t, byid=T) -> heathrow.area.75
gIntersection(heathrow.diff.50, heathrow.area.t, byid=T) -> heathrow.area.50
gIntersection(heathrow.diff.25, heathrow.area.t, byid=T) -> heathrow.area.25

Test Plot 2

I’ve included another rough plot of the newly calculated areas with an overlay of the 250 meter buffer to show the difference that accounting for excavation area will make on the calculations that will follow.

plot(heathrow.area.250, col="purple")
plot(heathrow.area.150, col="blue", add=TRUE)
plot(heathrow.area.100, col="green", add=TRUE)
plot(heathrow.area.75, col="yellow", add=TRUE)
plot(heathrow.area.50, col="orange", add=TRUE)
plot(heathrow.area.25, col="red", add=TRUE)
plot(heathrow.diff.250, add=TRUE)

Excavated Area Calculations

The excavated areas have now been plotted by buffer zone based on distance to the structural features. It is now a simple calculation to find the total excavated area within each buffer zone.

gArea(heathrow.area.250) -> area.250.calc
gArea(heathrow.area.150) -> area.150.calc
gArea(heathrow.area.100) -> area.100.calc
gArea(heathrow.area.75) -> area.75.calc
gArea(heathrow.area.50) -> area.50.calc
gArea(heathrow.area.25) -> area.25.calc

Spatial Overlay with Pottery Data

The sp::over function will create a vector for each buffer zone with “NA” values if the pottery was outside the buffer area and numeric values if the pottery was found within the given area.

over(heathrow.pottery.t, heathrow.area.250) -> pottery.250
over(heathrow.pottery.t, heathrow.area.150) -> pottery.150
over(heathrow.pottery.t, heathrow.area.100) -> pottery.100
over(heathrow.pottery.t, heathrow.area.75) -> pottery.75
over(heathrow.pottery.t, heathrow.area.50) -> pottery.50
over(heathrow.pottery.t, heathrow.area.25) -> pottery.25

Find the Total Pottery by Buffer Zone

This section of code combines the sherd count and weight columns with the vectors created above based on whether the pottery fell within each individual buffer zone. The subset function will make it possible to omit rows with “NA” in any specified column in the calculation of the total pottery weight. The total pottery weight for the 50 meter zone, for instance, will only be calculated with the subset of rows for which the value of pottery.50 is numeric.

weight.df <- data.frame(heathrow.pottery[,22:23], pottery.25, pottery.50, pottery.75, pottery.100, pottery.150, pottery.250)
sum(weight.df$WEIGHT[!is.na(weight.df$pottery.250)]) -> sum.250
sum(weight.df$WEIGHT[!is.na(weight.df$pottery.150)]) -> sum.150
sum(weight.df$WEIGHT[!is.na(weight.df$pottery.100)]) -> sum.100
sum(weight.df$WEIGHT[!is.na(weight.df$pottery.75)]) -> sum.75
sum(weight.df$WEIGHT[!is.na(weight.df$pottery.50)]) -> sum.50
sum(weight.df$WEIGHT[!is.na(weight.df$pottery.25)]) -> sum.25

Pottery per Square Meter

Using the calculations from above, it is now possible to calculate the amount of pottery per square meter for each respective buffer distance from the structural features.

sum.250/area.250.calc -> ppm.250
sum.150/area.150.calc -> ppm.150
sum.100/area.150.calc -> ppm.100
sum.75/area.75.calc -> ppm.75
sum.50/area.50.calc -> ppm.50
sum.25/area.25.calc -> ppm.25

#create a data frame with the pottery per square meter measurements
dist <- c(25, 50, 75, 100, 150, 250)
ppm <- c(ppm.25, ppm.50, ppm.75, ppm.100, ppm.150, ppm.250)
ppm.df <- data.frame(dist, ppm)
ppm.df
##   dist       ppm
## 1   25 5.6145285
## 2   50 1.9256627
## 3   75 1.0798820
## 4  100 0.1675445
## 5  150 0.0000000
## 6  250 1.6713298

Transform Spatial Polygons to WGS84

In order to visualize the results of the buffer zone and calculations on an interactive map, it is first necessary to change the coordinate system to WGS84: EPSG4326.

spTransform(heathrow.area.250, CRS("+proj=longlat +datum=WGS84")) -> heathrow.250.wg
spTransform(heathrow.area.150, CRS("+proj=longlat +datum=WGS84")) -> heathrow.150.wg
spTransform(heathrow.area.100, CRS("+proj=longlat +datum=WGS84")) -> heathrow.100.wg
spTransform(heathrow.area.75, CRS("+proj=longlat +datum=WGS84")) -> heathrow.75.wg
spTransform(heathrow.area.50, CRS("+proj=longlat +datum=WGS84")) -> heathrow.50.wg
spTransform(heathrow.area.25, CRS("+proj=longlat +datum=WGS84")) -> heathrow.25.wg

Leaflet Maps

The final result of the above calculations will be shown on two leaflet maps. I decided to have one leaflet map without the OSM basemap and one with the OSM basemap because the OSM basemap of the excavated area is busy and potentially distracting but I also wanted to have a sense of scale on the visualization.

I plotted the excavation areas within each buffer zone in the order of the rainbow with the closest features in red and furthest in purple. I included a layer of the structural features in gray. There is also a popup feature for each buffer zone with the pottery per square meter calculation.

leaflet() %>% addPolygons(data=heathrow.250.wg, color="purple", popup="1.67g per square meter") %>%
  addPolygons(data=heathrow.150.wg, color="blue", popup="0.00g per square meter") %>%
  addPolygons(data=heathrow.100.wg, color="green", popup="0.17g per square meter") %>%
  addPolygons(data=heathrow.75.wg, color="yellow", popup="1.08g per square meter") %>% 
  addPolygons(data=heathrow.50.wg, color="orange", popup="1.93g per square meter") %>%
  addPolygons(data=heathrow.25.wg, color="red", popup="5.61g per square meter") %>%
  addPolygons(data=heathrow, color="grey") %>% addTiles()

leaflet() %>% addPolygons(data=heathrow.250.wg, color="purple", popup="1.67g per square meter") %>%
  addPolygons(data=heathrow.150.wg, color="blue", popup="0.00g per square meter") %>%
  addPolygons(data=heathrow.100.wg, color="green", popup="0.17g per square meter") %>%
  addPolygons(data=heathrow.75.wg, color="yellow", popup="1.08g per square meter") %>% 
  addPolygons(data=heathrow.50.wg, color="orange", popup="1.93g per square meter") %>%
  addPolygons(data=heathrow.25.wg, color="red", popup="5.61g per square meter") %>%
  addPolygons(data=heathrow, color="grey")

Conclusions

This spatial analysis of pottery distribution at a dispersed rural settlement at the present site of London Heathrow Terminal 5 provided helpful results on the spatial distribution of pottery. Unsurprisingly, the highest pottery distribution density was in the area immediately surrounding the medieval structural remains. There was still a substantial amount of pottery recovered within 50m of the structural remains. The excavated features within the 100m and 150m areas yielded 0.168g per square meter and 0.000 g per square meter respectively, showing that excavations at these distances from the occupation structures produced little pottery. The higher pottery density represented at 250m away from the structural features could represent disposal of waste at the settlement boundaries (since many of these features are ditches) or disposal of waste in the fields through manuring.

The London Heathrow Terminal 5 pottery distribution densities provide interesting data, but these should not be taken as an absolute measure of expected pottery distribution in rural medieval settlement space. The excavators only provided finds data for pottery recovered within features. Medieval pottery could have survived at higher rates within these features since it would be less likely to be subjected to disturbance through trampling or yard maintenance. Archaeologists sampling currently occupied rural settlements with 1 meter x 1 meter test pits cannot strip the entire surface and excavate only within the visible features. As a result, archaeologists should not expect to be able to replicate these pottery distribution densities through test pits.

Another cause for caution is that the amount of medieval structural features lost to post-medieval and contemporary construction is unknown. The fact that there are multiple clusters of medieval structural evidence suggests that there were probably additional structures within this area.

Although the numbers derived from this spatial analysis should not be taken as directly replicable, they provide potentially helpful guidelines in currently occupied rural settlement sampling. The minimal pottery recovered from features between 75 and 150 meters away from structural features suggests that samples at these distances from the street frontages of currently occupied rural settlements will not be very helpful to inferring occupation. The steady decrease of pottery distribution density with increased distance from the settlement structures suggests that test pits within 50 meters of street frontages will be most informative of occupation history.