just a placeholder for now, this will all change - notes, and see gdalio
## these belong in raadfiles, because it should create and cache them
## but here they could be data?
# rema_8m_dem_vrt_file
# rema_8m_filled_vrt_file
# rema_8m_slope_vrt_file
# rema_8m_aspect_vrt_file
# rema_8m_rugosity_vrt_file
#
# rema_100m_dem_vrt_file
# rema_100m_filled_vrt_file
# rema_100m_slope_vrt_file
# rema_100m_aspect_vrt_file
# rema_100m_rugosity_vrt_file
rema_200m_dem_vrt_file <- rema_200m_dem_files()$fullname[1]
# rema_200m_filled_vrt_file
# rema_200m_slope_vrt_file
# rema_200m_aspect_vrt_file
# rema_200m_rugosity_vrt_file
.extent_around <- function(x, radius = 1) {
if (length(x) == 2) x <- matrix(x, ncol = 2)
xx <- range(x[,1])
yy <- range(x[,2])
ex <- raster::extent(xx + c(-1, 1) * radius,
yy + c(-1, 1) * radius)
}
.polarcrs <- "+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
.extent_rema <- function(x, buffer = NULL) {
if (is.numeric(x)) {
if (!length(x) == 2 && !ncol(x) == 2) stop("grid point/s must be a single lon,lat point or a matrix of points")
if (is.null(buffer)) {
if (ncol(x) == 2) {
return(raster::extent(reproj::reproj(x, .polarcrs, source = "+proj=longlat +datum=WGS84")[,1:2, drop = FALSE]))
}
stop("if grid point supplied, 'buffer' must also be included")
}
px <- reproj::reproj(x, .polarcrs, source = "+proj=longlat +datum=WGS84")
ext <- .extent_around(px, buffer)
return(ext)
}
}
read_rema_200m_dem <- function(x, buffer = NULL, max_dim = NULL) {
## x is a lon,lat point, buffer must be included - return native
## x is an extent, check if lonlat, use max_dim - return longlat
## x is an extent, check if proj, use max_dim - return native
## x is a grid - return grid
ext <- .extent_rema(x, buffer)
## need lazyraster logic here
#actual_grid <- vapour::vapour_raster_info(rema_200m_dem_vrt_file)
lr <- lazyraster::crop(lazyraster::lazyraster(rema_200m_dem_vrt_file), ext)
actual_dim <- c(diff(lr$window$window[1:2]), diff(lr$window$window[3:4]))
## we might be asking for more than is there, so we can't use lazyraster logic (need to get more abstract, pull out of lazyr)
actual_ext <- raster::extent(lr$window$windowextent)
if (!is.null(max_dim)) {
max_dim <- rep(max_dim, length.out = 2)
actual_dim[1] <- min(actual_dim[1], max_dim[1])
actual_dim[2] <- min(actual_dim[2], max_dim[2])
}
target <- raster::raster(actual_ext, nrows = actual_dim[2], ncols = actual_dim[1])
## actual_ext might be updated because of user input
v <- vapour::vapour_warp_raster(rema_200m_dem_vrt_file, extent = raster::extent(target), dimension = actual_dim)
raster::setValues(target, v[[1L]])
}
# read_rema_200m_filled
# read_rema_200m_slope
# read_rema_200m_aspect
# read_rema_200m_rugosity
just a placeholder for now, this will all change - notes, and see gdalio