knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, comment = "#>") knitr::opts_knit$set(root.dir = "/Users/gnj/Dropbox/R-Projects/map.displ.r/inst/extdata/")
The BUnwarpJ plugin in the open-source image processing package ImageJ/Fiji can be a very useful tool for high-resolution surface deformation modelling from digital elevation models (DEMs). The map.disp.r
package was built for importing bUnwarpJ image registration based on DEMs and computing displacements in the x,y and z directions.
This vignette
As a case study, we will map the deformation patterns of a rock glaicer using DEMs that were derived from lidar and structure-from-motion topographic surveying (Goetz et al 2019).
To start, we will load our DEMs for the years 2012 and 2017 that were used for performing image registration in BUnwarpJ for displcament mapping.
library(raster) library(sp) dem2012 <- raster("rock_glacier_dem_2012_50cm.tif") dem2017 <- raster("rock_glacier_dem_2017_50cm.tif")
The raw transformation file from BUnwarpJ describes the registration by providing the coordinates for each pixel of the source image and the corresponding row and column position of the target image. To apply this transformation to georeferenced images, we need to convert these positions to a local CRS. The dem.displacement.mapping()
function does this for us by using the source and target DEMs used for performing image registration. Additionally, it computes the vertical displacement (in the z direction), as well as the direction (aspect) of the 2D (xy) displacements and the slope of the z displacement.
library(map.displ.r) f_tx <- "RAW_rock_glacier_hillshade_2012_50cm_direct_transf.txt" d_tx <- dem.displacement.mapping(tx_file = f_tx, r_source = dem2012, r_target = dem2017) head(d_tx)
Using raster
we can assign the values of the displacement mapping to a grid. This is useful for exporting.
# Assign displacement values to an emptry raster to export r_na <- setValues(dem2012, NA) disp_3d <- setValues(r_na, d_tx$xyz_disp) plot(disp_3d)
For visualization of the displacements in R, we can use ggplot2
combined with the metR
and ggnewscale
packages. metR
can create arrows indicating the magnitude and direction of displacement vectors, and ggnewscale
lets us plot multiple geom_raster()
's each with their own fill scale and legend.
library(ggplot2) library(metR) library(ggnewscale) # Load hillshade for our map of displacements hs2017 <- raster("rock_glacier_hillshade_2017_50cm.tif") hs2017_df <- as.data.frame(hs2017, xy=TRUE) names(hs2017_df) <- c("x", "y", "hs") # Estimate mean annual surface velocity (m/yr) over the five year period d_tx$mean_disp <- d_tx$xyz_disp/5 # Make displacement map map <- ggplot(d_tx, aes(x_source,y_source)) + geom_raster(data=hs2017_df, aes(x=x, y=y, fill = hs), show.legend = FALSE) + scale_fill_gradient(high = "white", low = "black", na.value = "#FFFFFF") + # Allow for multiple scale fills using ggnewscale package new_scale("fill") + geom_raster(data=d_tx, alpha = 0.4, aes(x=x_source, y=y_source, fill = mean_disp)) + scale_fill_viridis_c(name = "Mean annual\nsurface velocity\n(m/yr)", direction = 1) + # Create arrows using metR package geom_arrow(data=d_tx, aes(dx = x_disp, dy = y_disp), skip = 30, show.legend = FALSE) + xlab("Easting (m)") + ylab("Northing (m)") + coord_fixed() + theme(text = element_text(size = 9), axis.title = element_text(size = 9), axis.text = element_text(size = 6), axis.text.y = element_text(angle = 90)) map
Arganda-Carreras, I., Sorzano, C. O., Marabini, R., Carazo, J. M., Ortiz-de-Solorzano, C., & Kybic, J. (2006, May). Consistent and elastic registration of histological sections using vector-spline regularization. In International Workshop on Computer Vision Approaches to Medical Image Analysis (pp. 85-95). Springer, Berlin, Heidelberg.
Goetz, J., Fieguth, P., Kasiri, K., Bodin, X., Marcer, M., & Brenning, A. (2019). Accounting for permafrost creep in high-resolution snow depth mapping by modelling sub-snow ground deformation. Remote Sensing of Environment, 231, 111275.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.