rema_longlat.md

Get rema in longlat for a given extent, and resolution.

2m resolution is some thing like 1/50000 of a degree, so 3e-5 gives cells with 4m^2 area at 70S

#' Get REMA v2 data in longlat
#' 
#' This is crafted VRT hosted on github that references the source 2m GeoTIFFs and some pre-rendered lower resolutions, 
#' so we can make data from large areas quickly as well as small areas. 
#' 
#' Note that AWS_NO_SIGN_REQUEST=YES must be set, i.e. 'Sys.setenv(AWS_NO_SIGN_REQUEST="YES")'
#' @param extent xmin,xmax,ymin,ymax in longlat values
#' @param res resolution in degrees (see Details)
#' @param method bilinear use bilinear resampling when changing resolution (set to 'near' for no resampling)
#' @param ... arguments passed to 'terra::project()' such as "filename"
#' @examples
#' rema_longlat(c(100, 120, -70, -60), res = .1, filename = "myrema.tif")
#' # rema_longlat(c(113.95, 114, -66.32, -66.3), res = c(3e-6, 3e-5), filename = "local.tif", threads = TRUE)
rema_longlat <- function(extent = c(-180, 180, -90, -65), res = c(1, 1), method = "bilinear", ...) {
  r <- terra::rast("/vsicurl/https://raw.githubusercontent.com/mdsumner/rema-ovr/main/REMA-2m_dem_ovr.vrt")
  terra::project(r, terra::rast(terra::ext(extent), res = res, crs = "EPSG:4326"), by_util = TRUE, method = method, ...)
}

## always set this for high resolution requests
Sys.setenv(AWS_NO_SIGN_REQUEST="YES")
library(terra)
plot(rema_longlat(c(100, 120, -70, -60), res = .1, filename = "myrema.tif"))


plot(rema_longlat(c(113.95, 114, -66.32, -66.3), res = c(3e-6, 3e-5), filename = "local.tif", threads = TRUE))

plot(rema_longlat(c(100, 120, -70, -60), res = .1, filename = "myrema.tif"))

plot(rema_longlat(c(113.95, 114, -66.32, -66.3), res =  3e-5, filename = "local.tif", threads = TRUE))

^{Created on 2024-12-11 with reprex v2.0.2}

Another example, Scullens monolith

plot(rema_longlat(rep(c(66.71886, -67.79353), each = 2) + c(-1, 1, -1, 1) * 2000/111120, res = 3e-5, filename = "scullen.tif", overwrite = T))

Get REMA from a KML boundary.

#' Get REMA v2 data in longlat based on a vector layer. 
#'
#' For input use 'terra::vect()' to open a 'shapfile':  KML, SHP, TAB, GPKG, FGB, Parquet etc. 
#' 
#' @param x SpatVector object 
#' @param buffer a buffer distance to add in metres
#' @param filename name of a .tif file to be created (will be a tempfile if not provided)
#' @return filename
#' @export
#'
rema_extent <- function(x, buffer = 300, filename = tempfile(fileext = ".tif")) {
  dsn <- "/vsicurl/https://raw.githubusercontent.com/mdsumner/rema-ovr/main/REMA-2m_dem_ovr.vrt"
  buffer <- rep(buffer, length.out = 2L)
  e <- as.vector(terra::project(terra::ext(x), to = "EPSG:4326", from = terra::crs(x)))
  ## buffer in metres
  bufy <- buffer[2L] / (1852 * 60)
  bufx <- bufy * cos(mean(e[3:4]) * pi/180)
  e <- unname(e + c(-bufx, bufx, -bufy, bufy))
  gdalraster::warp(dsn, filename, t_srs  = "EPSG:4326", cl_arg = c("-te", e[1], e[3], e[2], e[4]))
  filename
}
## always set this for high resolution requests
Sys.setenv(AWS_NO_SIGN_REQUEST="YES")
library(terra)

v <- vect("Ops Area (Job #0001854).kml")
tif <- rema_extent(v)

mess around with 3D a bit

library(anglr)
library(rgl)
m <- anglr::as.mesh3d(raster::raster(rast(tif)))

sv <- copy_down(anglr::PATH0(as(v, "Spatial")), raster::raster(rast(tif)))
clear3d()
plot3d(sv, col = "red", lwd = 4)
wire3d(m, col = "grey")
aspect3d(.4, 1, 1/10)
axes3d()
rgl::rglwidget()