LiDAR in Rcheology

In this short tutrial I am scripting a workflow for working with LiDAR data downloaded from Scottish Remote Sensing Portal. The aim is to process the raw LAZ files, produce visualizations and take simple measurments solely relaying on R programming language. I am looking at a 12th century motte in Cally Forrest, Galloway, in south-west Scotland. This mound was actually surveyed by Rubicon Heritage in 2012 clik here to read about the fieldwork. Now lets see how it compares with work of an armchair archaeologist

Packages

First, I am loading all packages necessary for the following analysis.

library(lidR) # process LAZ files
library(raster)# deal with raster
library(rayshader) # 3D viz
library(rgl) # interactive plot
library(sf) # spatial class
library(mapedit) # interactive map editing
library(mapview) # interactive map viewing
library(tmap) # cartography
library(tmaptools) # little helpers
library(smoothr) # smooth contours
library(dplyr) # manipulate
library(ggplot2) # plot
library(grid) # combines plots

LAZ

Next, I am reading the LAZ files and clipping the point cloud to the extent of immidiate area around the motte.

#read laz files
las = readLAS("data/NX6055_4PPM_LAS_PHASE3.laz")

# read a shapefile of area around the motte
motte_area = st_read("data/motte_area.shp", quiet = TRUE)

#clip the laz file to the extend of aoi
motte = lasclip(las, motte_area)

#plot lidar point cloud
plot(motte, bg = "white")

rglwidget()

grab and rotate the model !

Digital Elevetion Model

DTM and DSM

In this step I am running standard alghoritms from LiDR package to compute Digital Surface Model (DSM) and Digitial Terrain Model (DTM).

rgl.clear()
# create dsm and dtm
dsm = grid_canopy(motte, 0.5, pitfree())

# assign coordinate system
crs(dsm) = CRS("+init=epsg:27700")

# create dtm
dtm = grid_terrain(motte, 0.5, algorithm = knnidw(k = 6L, p = 2))

# addign coordinate system
crs(dtm) = CRS("+init=epsg:27700")

par(mfrow = c(1,2))
plot(dtm, main = "DTM", col = height.colors(50), legend = FALSE)
plot(dsm, main = "DSM",col = height.colors(50))

Hillshade

There are dozens of alghoritms to visualize or extract characteristics from Digital Elevations Models. The most common is hillshade determined by the aspects and slope of the terrain as well as light source.

# dsm
slope_dsm = terrain(dsm, opt = 'slope')
aspect_dsm = terrain(dsm, opt = 'aspect')
hill_dsm = hillShade(slope_dsm, aspect_dsm, angle = 40, direction = 270)

# dtm
slope_dtm = terrain(dtm, opt = 'slope')
aspect_dtm = terrain(dtm, opt = 'aspect')
hill_dtm = hillShade(slope_dtm, aspect_dtm, angle = 5, direction = 270)

#plot
par(mfrow = c(1,2))
plot(hill_dtm, main = "DTM Hilshade", col = grey.colors(100, start = 0, end = 1), 
     legend = FALSE)
plot(hill_dsm, main = "DSM Hillshade", col = grey.colors(100, start = 0, end = 1))

Transect

In order to represent the height values and get idea about the profile I will draw a line across the motte.

Line

Package mapedit allows to interactively draw the line and save it as sf object.

mapviewOptions(raster.palette = grey.colors)

interactive_map <- mapview(hill_dtm) %>%
  editMap()

transect.sf = interactive_map$finished

Elevation

In the following code I use the line selector to extract the z-values from the raster. I need to give credit for the idea to the book Geocomputation with R. Although I had to amend the code to work in projected coordinate system.

transect = raster::extract(dtm, transect.sf, along = TRUE, cellnumbers = TRUE)
transect_df = purrr::map_dfr(transect, as_data_frame, .id = "ID")
transect_coords = xyFromCell(dtm, transect_df$cell)
pair_dist = pointDistance(transect_coords,lag(transect_coords), lonlat=FALSE)
transect_df$dist = c(0, cumsum(na.omit(pair_dist))) 

Transect Plot

Now I can plot the data frame with heights along the transect line. The tip is to use fixed aspect ratio in order to preserve physical representation of the data units on the axes.

profile = ggplot(transect_df) +
  geom_area(aes(x = dist, y = Z),  alpha=0.6 , size=1, colour="#AF002A", fill = "#AF002A") +
  scale_y_continuous(breaks = seq(40, 50, by = 5)) +
  scale_x_continuous(breaks = seq(0, 100, by = 10)) +
  coord_fixed(ylim = c(40,50)) +
  labs(title = "Profile",
       x = "Distance (m)",
       y = "Elevation (m)") +
  theme_bw() +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.border = element_blank())
profile

Map

Lets create a topographic map of the motte and plot the location of the transect.

#load canomre area
topo_contour <- rasterToContour(dtm, nlevels = 10) %>% 
  st_as_sf() %>% 
  st_cast("LINESTRING") %>% 
  drop_crumbs(20)

## Warning in st_cast.sf(., "LINESTRING"): repeating attributes for all sub-
## geometries for which they may not be constant

# color palette
brewer_grey = get_brewer_pal("Greys", n = 10, contrast = c(0.22, 0.82), plot = FALSE)

map = tm_shape(dtm, unit = "m") +
  tm_raster(style= "pretty", n = 10, palette = brewer_grey, legend.show=F, alpha = 0.75) +
tm_shape(topo_contour) +
  tm_lines(col = "level",
           palette = "Greys",
           legend.col.show = F) +
  tm_text(text = "level",
          size = .5,
          col = "level",
          palette = "-Greys",
          along.lines = F, 
          overwrite.lines = T,
          legend.col.show = F) +
tm_shape(transect.sf) +
  tm_lines(col = "#AF002A",
           lty = "dashed",
           lwd = 2) +
tm_compass(position = c("left", "bottom")) +
tm_scale_bar(breaks = c(0,5,10),
             position = c("left", "bottom")) +
tm_credits("@topographos2") +
tm_layout(
    main.title = "Cally Forrest Motte",
    frame = F)

map

Final Graphic

Finally, I can combine the two graphics.

# plot using grid
library(grid)

grid.newpage()
# set up the layout
pushViewport(viewport(layout=grid.layout(1,2)))
# plot using the print command
print(map, vp=viewport(layout.pos.col = 1, height  = 20))
print(profile, vp=viewport(layout.pos.col = 2, height = 20))

Going 3D

And a cherry on the top by amazing rayshader package

#And convert it to a matrix:
elmat = raster_to_matrix(dtm)

elmat %>%
  sphere_shade(texture = "imhof1") %>%
  add_shadow(ray_shade(elmat, zscale = 0.5, sunaltitude = 30,lambert = TRUE),
             max_darken = 1) %>%
  add_shadow(lamb_shade(elmat,zscale = 0.5,sunaltitude = 30), max_darken = 0.2) %>%
  add_shadow(ambient_shade(elmat, zscale = 0.5), max_darken = 0.2) %>%
  plot_3d(elmat, zscale = 0.5,windowsize = c(1000,1000),
          phi = 40, theta = 180, zoom = 0.8, fov = 1)
rglwidget()

grab and rotate the model !

sessionInfo()

## R version 4.0.0 (2020-04-24)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] shiny_1.4.0.2    ggplot2_3.3.0    dplyr_0.8.5      smoothr_0.1.2   
##  [5] tmaptools_3.0    tmap_3.0         mapview_2.7.8    mapedit_0.6.0   
##  [9] sf_0.9-2         rgl_0.100.54     rayshader_0.15.1 lidR_2.2.3      
## [13] raster_3.1-5     sp_1.4-1        
## 
## loaded via a namespace (and not attached):
##  [1] jsonlite_1.6.1          viridisLite_0.3.0       foreach_1.5.0          
##  [4] assertthat_0.2.1        stats4_4.0.0            yaml_2.2.1             
##  [7] progress_1.2.2          pillar_1.4.3            lattice_0.20-41        
## [10] glue_1.4.0              digest_0.6.25           manipulateWidget_0.10.1
## [13] RColorBrewer_1.1-2      promises_1.1.0          leaflet.providers_1.9.0
## [16] colorspace_1.4-1        htmltools_0.4.0         httpuv_1.5.2           
## [19] XML_3.99-0.3            pkgconfig_2.0.3         magick_2.3             
## [22] stars_0.4-1             purrr_0.3.4             xtable_1.8-4           
## [25] scales_1.1.0            webshot_0.5.2           satellite_1.0.2        
## [28] later_1.0.0             leaflet.extras_1.0.0    tibble_3.0.1           
## [31] farver_2.0.3            ellipsis_0.3.0          withr_2.2.0            
## [34] lazyeval_0.2.2          leafsync_0.1.0          magrittr_1.5           
## [37] crayon_1.3.4            mime_0.9                evaluate_0.14          
## [40] doParallel_1.0.15       lwgeom_0.2-3            class_7.3-16           
## [43] tools_4.0.0             data.table_1.12.8       prettyunits_1.1.1      
## [46] hms_0.5.3               rlas_1.3.5              lifecycle_0.2.0        
## [49] stringr_1.4.0           munsell_0.5.0           compiler_4.0.0         
## [52] e1071_1.7-3             rlang_0.4.5             classInt_0.4-3         
## [55] units_0.6-6             dichromat_2.0-0         iterators_1.0.12       
## [58] htmlwidgets_1.5.1       crosstalk_1.1.0.1       miniUI_0.1.1.1         
## [61] leafem_0.1.1            base64enc_0.1-3         rmarkdown_2.1          
## [64] gtable_0.3.0            codetools_0.2-16        abind_1.4-5            
## [67] DBI_1.1.0               markdown_1.1            R6_2.4.1               
## [70] rgdal_1.4-8             knitr_1.28              rgeos_0.5-2            
## [73] fastmap_1.0.1           KernSmooth_2.23-16      stringi_1.4.6          
## [76] parallel_4.0.0          Rcpp_1.0.4.6            vctrs_0.2.4            
## [79] png_0.1-7               tidyselect_1.0.0        leaflet_2.0.3          
## [82] xfun_0.13

 

---

A work by Michal Michalski

michal.m.michalski@durham.ac.uk