Landscape modeling
Thomas Huet
March 19, 2021
Abstract
A data-driven presentation linked to the ANTHEA ERC projectPresentation
Concepts of landscape, space, and place, came from the research of Lefebvre (1974) and Tilley (1994). According to Tilley, the landscape is a “serie of named locales, a set of relational places linked by paths, movements and narratives. It is a ‘natural’ topography perspectivally linked to the existential being of the body in social space.” For Anschuetz et al., “Landscape archaeology can be said to be the archaeology of place” (Anschuetz, Wilshusen, and Scheick 2001). Place “is a particular kind of space” (Low 2016). It results from space social valuation (Ingold 1993; Meskell and Preucel 2004). Landscape archaeology aims to render places social uses through models.
The Saint Michel church on the top of a Neolithic tumulus (Carnac, Bretagne, France). A clear example of time persistence between a past structure and its effects in the present landscape
Places are different from archaeological sites. Archaeological sites have a delimited duration and extinguished geographical effects. Instead, places are durable constructions, resilient structures determining future uses (path dependency) and creating palimpsest over the long-term, between past forms and present usages – what Chouquer (2000) calls “hystérochronie.”
Tilley (1994, 2004) and Ingold (1993) also introduce the concepts of “phenomenological landscape.” The landscape is perceived by the human vision (visual scape), audition (sound scape), social geography (task scape), etc. Landscape and quantitative archaeology have already taken over viewsheds and inter-visibilities (Wheatley 1995) and social geography is studied since the very beginning of the archaeological discipline
IT settings
This HTML web page presentation is hosted on this GitHub repository. It results from an R report generation. The source document is a Rmarkdown file using mainly packages like Leaflet and Plotly to create interactive and enriched maps and graphics. We will mostly manage raster files, so the R package raster is at the root of most data management.
This presentation attempts to illustrate the reuse of existing data for landscape modeling. All used data are open source, directly read by their URL or downloaded and stored in the GitHub folder:
https://github.com/zoometh/popland/tree/main/doc/data
To simplify the presentation, part of the R code used to generate the graphics and maps is shown and the analysis mostly focuses on Northern Europe between 2800 BC and 1000 AD
Figure 1: The four dimensions of the study with the list of studied values for the Northern Europe landscape
Lines of evidences
Lines of evidences of space anthropogenic management are controlled fires, fencing systems, palynological evidence, etc.
Palynological evidences
Some pollen ratios, like phosphate concentrations, are clear indexes of land anthropogenic management (manuring, etc.) (Odgaard 1985; Nielsen et al. 2019). For a better understanding of the sampling context (e.g. in- or outside the agrarian structure), sampling locations and – eventual OSL dates – should be precisely mapped, for example on a LiDAR digital surface model, with ground-truthing (Nielsen et al. 2019)
Controlled fires
Since its domestication, the fire helps human societies to manage their environment (niche construction) (Scott 2017). In Northern Europe, more than 4,000 years ago, agropastoral communities began the first fire-based management (“Anthropogenic Heathlands,” n.d.). The Global Charcoal Database (GCD) gather thousand of sedimentary records of fire. This database can be read with the R package paleofire (Blarquez et al. 2014). The key field ‘id_site’ allows connecting each database record to the root URL of the GDC database (with the package kableExtra). Here we subset the data coming from Northwest Europe
id.url <- "https://paleofire.org/index.php?p=CDA/site_view&gcd_menu=CDA&site_id="
ID <- pfSiteSel(lat > 49 & lat < 59,
long > -6 & long < 18)
sumID <- summary(ID)
samleID <- sumID[sample(nrow(sumID), 10), ]
samleID$url <- paste0(id.url, samleID$id_site)
samleID %>%
mutate(id_site = cell_spec(id_site, "html", link = url)) %>%
dplyr::select (-url) %>%
kable("html", escape = FALSE,
caption = "Sample of mtDNA metatadata with hyperlinks") %>%
kable_styling(bootstrap_options = c("hover", "condensed"), font_size=11)| id_site | lat | long | elevation | min_est_age | max_est_age | num_dating | date_int | num_samp | l12 | rf99 | methods | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Black Loch | 242 | 56.00 | -3.00 | -9999 | -45.00 | 10760.27 | NA | NA | 56 | 0 | 0 | POLS |
| Ballachrink | 479 | 54.37 | -4.47 | 20 | 6596.00 | 9379.00 | 5 | 556.6000 | 52 | 5 | 0 | POLS |
| Cerna Horna | 942 | 50.66 | 15.75 | 1190 | 2866.40 | 10048.00 | NA | NA | 53 | 5 | 5 | NA |
| West Lomond | 969 | 56.25 | -3.28 | 522 | 494.00 | 10012.00 | NA | NA | 66 | 5 | 5 | NA |
| Bolton Fell Moss | 656 | 55.02 | -2.85 | 110 | 6008.19 | 9780.65 | 10 | 377.2460 | 87 | 5 | 5 | POLS |
| Migneint | 660 | 52.82 | -3.82 | 440 | -31.94 | 4586.07 | NA | NA | 72 | 5 | 5 | SIEV |
| Cothill Fen | 654 | 51.68 | -1.32 | -9999 | 3986.61 | 11139.69 | 4 | 1788.2700 | 49 | 5 | 5 | POLS |
| Kragso | 547 | 56.21 | 9.02 | 47 | 11.00 | 10233.00 | NA | NA | 59 | 5 | 5 | POLS |
| Stavsakra Bog | 470 | 56.91 | 14.90 | 187 | 25.00 | 10535.00 | 18 | 583.8889 | 77 | 5 | 4 | POLS |
| Skanso | 555 | 56.51 | 8.83 | 8 | 15.00 | 6215.00 | NA | NA | 72 | 5 | 5 | POLS |
The GCD web site already offers a web map interface (e.g. Jackfish archaeological site) but one can use the Leaflet package for further data integration
sumID$lbl <- paste0("<a href=",shQuote(paste0(id.url,sumID$id_site)),
"\ target=\"_blank\"", ">",
sumID$id_site, "</a><b> ", rownames(sumID), "</b> <br>",
"dating: ", sumID$max_est_age, " to ", sumID$min_est_age)
leaflet(data = sumID, width = "70%", height = "500px") %>%
addTiles(group = 'OSM') %>%
addCircleMarkers(layerId = ~id_site,
lng = ~long,
lat = ~lat,
weight = 1,
radius = 3,
popup = ~lbl,
opacity = 0.7,
fillOpacity = 0.7)Enclosure systems
Rectangular plots appear around 1900 BC in Northwestern Europe with a peak around 1600 BC (e.g. Wessex Celtic fields) (Løvschal 2014). Geographical data come from raster (satellite image, LiDAR, etc.) and shapefiles
Raster approach
A large part of landscape data is raster. The R package raster allows reading raster layers from .DEM, .tif and .asc files. To illustrate R capacities, we choose the environment of the archaeological site of Sidbury Hill (Wiltshire, England), an Iron Age hill fort close to Celtic fields. To set the geographical context of Sidbury Hill, one can use a Leaflet with a layer button to select between different backgrounds, here: OSM, Esri.WorldImagery and a WMS tile, the CORINE European Land Cover (this later is not avalaible at all scales)
Sidbury.long <- -1.691734639258629
Sidbury.lat <- 51.253526272001714
corine.url = "https://image.discomap.eea.europa.eu/arcgis/services/Corine/CLC2018_WM/MapServer/WmsServer"
leaflet(width = "60%", height = "400px") %>%
addTiles(group = 'OSM') %>%
addWMSTiles(corine.url, layers = "12",
options = WMSTileOptions(format = "image/png", transparent = T),
group = 'CORINE') %>%
addProviderTiles(providers$Esri.WorldImagery, group='Esri.WorldImagery') %>%
addMarkers(Sidbury.long,
Sidbury.lat,
label = "Sidbury Hill") %>%
addLayersControl(
baseGroups = c('OSM', 'Esri.WorldImagery', 'CORINE')) %>%
addScaleBar(position = "bottomleft")
Some LiDAR tiles are freely available online like those of the British Departement for Environment Food & Rural Affairs. The projection system is the British National Grid (aka BNG, Ordnance Survey National Grid, EPSG: 27700). The LiDAR tile of Sidbury Hill is su2150_DSM_2M.asc (BNG grid reference: SU216505). Once downloaded, we stored it into GitHub. After reading it with the raster package, its conversion into a data frame, it can be plotted with the RasterVis package
Figure 2: Sidbury Hill (Wiltshire, England) LiDAR tile
The site is characterized by a large double ditch and rampart enclosure. One can apply a sharpening filter to accentuate the changes in the slopes. The imager package, allows applying such a filter with the function imsharpen
par(mar = c(0, 0, 2, 0))
rdem.img <- as.cimg(rdem)
rdem.img %>%
plot(main="Normal", axes = F)
imsharpen(rdem.img, 1000) %>%
plot(main="Sharpened", axes = F)
From there, the geometry of the enclosure can be extracted by shapefile digitalization or raster segmentation
Shapefiles approach
Another common data format is the shapefiles. The open-access Landscape Dynamics (landDX) database (Amoke et al. 2020) offers a training dataset, by the mean of shapefiles, to analyze the settlement distribution. Here we sampled some polygons and polylines from Nakuru county (South Kenya)
Spatial distribution of settlements can be achieved in different ways, both at the inter-site scale, as well at the intra-site scale. Here, we will only consider the shape of the buildings (plots and fences). These latter are considered POLYGONS. One can use the R package stampr for shape analysis (Long, Robertson, and Nelson 2018) calculating their minimum bounding rectangle (MBR) and the minimum bounding circle (MBC). After this calculation, as an example, the first polygon can be plotted
par(mar = c(0, 0, 2, 0))
rd.i.cp <- data.frame(idf = numeric(0),
compactness = numeric(0),
roundness = numeric(0))
mbrs.list <- mbcs.list <- list()
for(p in 1:nrow(landDX.polygons)){
a.polyg <- landDX.polygons[p, ]
a.polyg.sf <- st_as_sf(a.polyg)
idf <- rownames(a.polyg@data)
# MBR
a.mbr <- st_geometry(st_as_sfc(st_bbox(a.polyg)))
a.mbr.sp <- as_Spatial(a.mbr)
a.mbr.spdf <- SpatialPolygonsDataFrame(a.mbr.sp,
data.frame(id=idf),
match.ID = F)
mbrs.list[[length(mbrs.list)+1]] <- a.mbr.spdf
# MBC
a.mbc <- st_minimum_bounding_circle(a.polyg.sf)
a.mbc.spdf <- as(a.mbc, "Spatial")
mbcs.list[[length(mbcs.list)+1]] <- a.mbc.spdf
}
mbrs.list.sp <- do.call(bind, mbrs.list) # list to dfs
mbcs.list.sp <- do.call(bind, mbcs.list) # list to dfs
plot(mbcs.list.sp[1, ], lwd = 1.5, border="blue",
main = row.names(landDX.polygons@data)[1])
plot(mbrs.list.sp[1, ], add=T, border="red", lwd = 1.5)
plot(landDX.polygons[1, ], add=T, lwd = 1.5, border="black")
Figure 3: Calculation of minimum bounding rectangle (MBR) and the minimum bounding circle (MBC), anf plot of the first polygon (black line) with its MBR (red line) and MBC (blue line)
The sf package offers many spatial functions, quite close to PostGIS ones. Here we extract the centroids of all features and select the centroid of the first one to set the Leaflet view (see below)
# mbrs.list.sp@polygons[[1]]@Polygons[[1]]@coords
sf_cent <- st_centroid(st_as_sf(mbrs.list.sp))
xy <- st_coordinates(sf_cent)[1,]
xy## X Y
## 35.727343 -1.688538
Now, we can spatialize these MBRs and MBCs and add them to the Leaflet map
Equivalently, the shape index compactness (i.e. how the shape fits a rectangle) and roundness (i.e. how the shape fits a circle) can be calculated and stored in a data frame
rd.i.cp <- data.frame(idf = numeric(0),
compactness = numeric(0),
roundness = numeric(0))
mbrs.list <- mbcs.list <- list()
for(p in 1:nrow(landDX.polygons)){
a.polyg <- landDX.polygons[p, ]
a.polyg.sf <- st_as_sf(a.polyg)
idf <- rownames(a.polyg@data)
a.compactness <- fs_compact(a.polyg.sf)
a.roundness <- as.numeric(fs_area(a.polyg.sf) / fs_area(fs_mbc(a.polyg.sf)))
rd.i.cp[nrow(rd.i.cp)+1, ] <- c(idf,
round(a.compactness, 2),
round(a.roundness, 2))
}
The rd.i.cp data frame list all compactness and roundness values. One can plot their distributions with a scatter plot. the Plotly package allows enriching the graphic by offering dynamic labeling tools like ‘hoverinfo’
Figure 4: Scatterplot of compactness and roundness values
Data integration
Landscape archaeology is strongly interdisciplinary and needs to put together different lines of evidence. Here, we reuse the Sidbury Hill LiDAR file and plot at its top identified gates. This data are stored in a shapefile of POINTS (.zip file) saved in GitHub:
https://raw.github.com/zoometh/popland/master/doc/data/sidbury/sidbury_gates.zip
temp <- tempfile()
download.file(paste0(root.data.sidbury, "sidbury_gates.zip"),
temp, quiet = TRUE)
db_path <- unzip(temp, exdir = tempdir())
SHP <- st_read(dsn = dirname(db_path)[1], layer = "features")## Reading layer `features' from data source `C:\Users\supernova\AppData\Local\Temp\Rtmp4WjMGM' using driver `ESRI Shapefile'
## Simple feature collection with 3 features and 2 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 421494.9 ymin: 150488.7 xmax: 421741.4 ymax: 150701
## projected CRS: OSGB 1936 / British National GridSHP.df <- as.data.frame(SHP)
SHP.df <- cbind(SHP.df, st_coordinates(SHP))
plot(rdem)
plot(SHP, pch=16, col="black", add=TRUE)
text(SHP.df[,"X"], SHP.df[,"Y"], SHP.df[,"text"], pos = 2)
Figure 5: Sidubury LiDAR DEM with the location of the main entrances
Data visualization
Today, new means in scripting languages, like R, allow going further in data visualization. The package rayshader allows one to create mutliscalar 3D views and files from a raster file, or an x, y, z matrix, and to enrich them with other data (POINTS, LINES, POLYGONS, etc.). As an example, we can reuse the Sidbury LiDAR file and the shapefile of the entrances within an enriched 3D plot
# Not run in Rmarkdown, see https://stackoverflow.com/questions/62417149/rayshader-in-rmarkdown-file
attr(rdem.mat, "extent") <- rdem@extent
rdem.mat <- raster_to_matrix(rdem)
rdem.mat %>%
sphere_shade(zscale = 1, texture = "imhof1") %>%
plot_3d(rdem.mat, zscale = 1, fov = 0, theta = 135, zoom = 0.6, phi = 45,
baseshape = "hex", windowsize = c(800, 800))
render_label(rdem.mat, long = SHP.df[1,"X"], lat = SHP.df[1,"Y"],
extent = attr(rdem.mat, "extent"), altitude = 50,
zscale = 1, text = SHP.df[1,"text"], textsize = 2, linewidth = 5)
render_label(rdem.mat, long = SHP.df[2,"X"], lat = SHP.df[2,"Y"],
extent = attr(rdem.mat, "extent"), altitude = 50,
zscale = 1, text = SHP.df[2,"text"], textsize = 2, linewidth = 5)
render_label(rdem.mat, long = SHP.df[3,"X"], lat = SHP.df[3,"Y"],
extent = attr(rdem.mat, "extent"), altitude = 50,
zscale = 1, text = SHP.df[3,"text"], textsize = 2, linewidth = 5)
render_scalebar(limits=c(0, 50, 100, 250), label_unit = "m",
position = "W", y = 150, scale_length = c(0.5, 1))
render_compass(position = "E", compass_radius = 40)
Sys.sleep(0.2)
render_snapshot(clear=TRUE)The result is a 3D png of the archaeological site:
3D view of the Sidbury Hill archaeological site with the R package rayshader