raster_mask_bench.R

# Raster Mask Benchmark
library(sf)
library(terra)
library(dplyr)
library(elevatr)
library(microbenchmark)
# library(tictoc)
library(bench) 

# ------ get data -------------
# Get county data for England and Wales
eng_wales_counties <- read_sf("http://geoportal1-ons.opendata.arcgis.com/datasets/687f346f5023410ba86615655ff33ca9_0.geojson") %>%
  st_make_valid() %>%
  st_transform(27700)
# plot(st_geometry(eng_wales_counties), axes=T)

# Union region for raster request
eng_wales <- eng_wales_counties %>%
  st_union()

# Get Wales sf for the masking task. Save for later
wales <- eng_wales_counties %>%
  filter(grepl("W",ctyua16cd) ) %>%
  st_union()

wales_path <- tempfile(fileext = '.gpkg')
write_sf(wales, wales_path)
# plot(st_geometry(wales), axes=T)

# Download Eleavation data with Elevatr
uk_terrain <- elevatr::get_elev_raster(eng_wales, src='alos')

# get source of raster
terrain_src <- uk_terrain@file@name

# plot of what we want...
plot(uk_terrain)
plot(st_geometry(wales), add=T)

# ------ functions ---------------------
mask_crop_terra <- function(ras_path, vec_path, out_path){
  big_R <- terra::rast(ras_path)
  
  crop_R <- terra::crop(big_R, vec_path)
  
  mask_R <- terra::mask(crop_R, vect(vec_path))
  
  terra::writeRaster(mask_R, out_path, overwrite=TRUE)
}

mask_terra <- function(ras_path, vec_path, out_path){
  big_R <- terra::rast(ras_path)
  
  mask_R <- terra::mask(big_R, vect(vec_path))
  
  terra::writeRaster(mask_R, out_path, overwrite=TRUE)
}


mask_gdalutils <- function(ras_path, vec_path, out_ras){
  if (file.exists(out_ras)) (file.remove(out_ras))
  
  sf::gdal_utils('warp',
                 source=ras_path,
                 destination = out_ras,
                 options = c('-cutline', vec_path,
                             '-crop_to_cutline', ras_path,
                             '-multi'))
}

# {stars} approach not working - RAM overloaded... # is there a better way?
# stars_mask <- function(ras_path, vec_path, out_ras){
#   s_ras <- stars::read_stars(ras_path)
#   poly <- sf::read_sf(vec_path)
#   st_crs(s_ras) <- st_crs(poly)
#   stars::write_stars(s_ras[poly], out_ras)
# }
# out_temp4 <-tempfile(fileext = '.tif')
# stars_mask(ras_temp, wales_path, out_temp4)


# ---------- benchmarking -------------
# out files
out_temp1 <- tempfile(fileext = '.tif')
out_temp2 <-tempfile(fileext = '.tif')
out_temp3 <-tempfile(fileext = '.tif')

# benchmarking
microbenchmark::microbenchmark(
  mask_gdalutils(terrain_src, wales_path, out_temp1),
  mask_crop_terra(terrain_src, wales_path, out_temp2),
  mask_terra(terrain_src, wales_path, out_temp3),
  times = 3L)

# RAM usage ( I assume this is representative but not 100% sure...)
bind_rows(
mark(x <- mask_gdalutils(terrain_src, wales_path, out_temp1))[,"mem_alloc"],
mark(x <- mask_crop_terra(terrain_src, wales_path, out_temp2))[,"mem_alloc"],
mark(x <- mask_terra(terrain_src, wales_path, out_temp3))[,"mem_alloc"])  %>%
  bind_cols(data.frame(func=c('mask_gdalutils', 'mask_crop_terra', 'mask_terra')))

# final plot (what we were aiming for)
plot(rast(out_temp1))
plot(st_geometry(wales), add=T)