.r_to_gdal_datatype <- function(x) {
if (nchar(x) == 0 || is.na(x)) return("")
xout <- toupper(x[1L])
xout <- c("RAW" = "Byte", "INTEGER" = "Int32", "DOUBLE" = "Float64", "NUMERIC" = "Float64",
"BYTE" = "Byte",
"UINT16" = "Int32", "INT16" = "Int32",
"UINT32" = "Int32", "INT32" = "Int32",
"FLOAT32" = "Float32", "FLOAT64" = "Float64")[xout]
if (is.na(xout)) {
message(sprintf("unknown 'band_output_type = \'%s\'', ignoring", x))
xout <- ""
}
xout
}
#' Raster IO (read)
#'
#' Read a window of data from a GDAL raster source. The first argument is the source
#' name and the second is a 6-element `window` of offset, source dimension, and output dimension.
#'
#' The value of `window` may be input as only 4 elements, in which case the source dimension
#' Will be used as the output dimension.
#'
#' This is analogous to the `rgdal` function `readGDAL` with its arguments `offset`, `region.dim`
#' and `output.dim`. There's no semantic wrapper for this in vapour, but see
#' `https://github.com/hypertidy/lazyraster` for one approach.
#'
#' Resampling options will depend on GDAL version, but currently 'NearestNeighbour' (default),
#' 'Average', 'Bilinear', 'Cubic', 'CubicSpline', 'Gauss', 'Lanczos', 'Mode' are potentially
#' available. These are compared internally by converting to lower-case. Detailed use of this is barely tried or tested with vapour, but is
#' a standard facility used in GDAL. Easiest way to compare results is with gdal_translate.
#'
#' There is no write support in vapour.
#'
#' Currently the `window` argument is required. If this argument unspecified and `native = TRUE` then
#' the default window specification will be used, the entire extent at native resolution. If 'window'
#' is specified and `native = TRUE` then the window is used as-is, with a warning (native is ignored).
#'
#' 'band_output_type' can be 'raw', 'integer', 'double', or case-insensitive versions of the GDAL types
#' 'Byte', 'UInt16', 'Int16', 'UInt32', 'Int32', 'Float32', or 'Float64'. These are mapped to one of the
#' supported types 'Byte' ('== raw'), 'Int32' ('== integer'), or 'Float64' ('== double').
#'
#' @param x data source
#' @param band index of which band to read (1-based)
#' @param window src_offset, src_dim, out_dim
#' @param resample resampling method used (see details)
#' @param ... reserved
#' @param native apply the full native window for read, `FALSE` by default
#' @param sds index of subdataset to read (usually 1)
#' @param set_na specify whether NA values should be set for the NODATA
#' @param band_output_type numeric type of band to apply (else the native type if ''), is mapped to one of 'Byte', 'Int32', or 'Float64'
#' @param unscale default is `TRUE` so native values will be converted by offset and scale to floating point
#' @param nara logical whether to return a (scaled) nativeRaster
#' @export
#' @return list of numeric vectors (only one for 'band')
#' @examples
#' f <- system.file("extdata", "sst.tif", package = "vapour")
#' ## a 5*5 window from a 10*10 region
#' vapour_read_raster(f, window = c(0, 0, 10, 10, 5, 5))
#' vapour_read_raster(f, window = c(0, 0, 10, 10, 5, 5), resample = "Lanczos")
#' ## find the information first
#' ri <- vapour_raster_info(f)
#' str(matrix(vapour_read_raster(f, window = c(0, 0, ri$dimXY, ri$dimXY)), ri$dimXY[1]))
#' ## the method can be used to up-sample as well
#' str(matrix(vapour_read_raster(f, window = c(0, 0, 10, 10, 15, 25)), 15))
#'
vapour_read_raster <- function(x, band = 1, window, resample = "nearestneighbour", ..., sds = NULL, native = FALSE, set_na = TRUE, band_output_type = "", unscale = TRUE, nara = FALSE) {
x <- .check_dsn_single(x)
if (!length(native) == 1L || is.na(native[1]) || !is.logical(native)) {
stop("'native' must be a single 'TRUE' or 'FALSE'")
}
band_output_type <- .r_to_gdal_datatype(band_output_type)
datasourcename <- sds_boilerplate_checks(x, sds = sds)
resample <- tolower(resample) ## ensure check internally is lower case
if (!resample %in% c("nearestneighbour", "average", "bilinear", "cubic", "cubicspline",
"gauss", "lanczos", "mode")) {
warning(sprintf("resample mode '%s' is unknown", resample))
}
ri <- vapour_raster_info(x, sds = sds)
if (native && !missing(window)) warning("'window' is specified, so 'native = TRUE' is ignored")
if (native && missing(window)) window <- c(0, 0, ri$dimXY, ri$dimXY)
if (!is.numeric(band) || band < 1 || length(band) < 1 || anyNA(band)) {
stop("'band' must be an integer of length 1, and be greater than 0")
}
if (any(band > ri$bands)) stop(sprintf("specified 'band = %i', but maximum band number is %i", band, ri$bands))
## turn these warning cases into errors here, + tests
## rationale is that dev can still call the internal R wrapper function to
## get these errors, but not the R user
stopifnot(length(window) %in% c(4L, 6L))
## use src dim as out dim by default
if (length(window) == 4L) window <- c(window, window[3:4])
## these error at the GDAL level
if (any(window[1:2] < 0)) stop("window cannot index lower than 0")
if (any(window[1:2] > (ri$dimXY-1))) stop("window offset cannot index higher than grid dimension")
## this does not error in GDAL, gives an out of bounds value
if (any(window[3:4] < 1)) stop("window size cannot be less than 1")
## GDAL error
if (any((window[1:2] + window[3:4]) > ri$dimXY)) stop("window size cannot exceed grid dimension")
## GDAL error
if (any(window[5:6] < 1)) stop("requested output dimension cannot be less than 1")
## pull a swifty here with [[ to return numeric or integer
##vals <- raster_io_gdal_cpp(dsn = datasourcename, window = window, band = band, resample = resample[1L], band_output_type = band_output_type)
vals <- lapply(band, function(iband) {
raster_io_gdal_cpp(dsn = datasourcename, window = window, band = iband, resample = resample[1L], band_output_type = band_output_type, unscale = unscale, nara = nara)[[1L]]
})
if (set_na && !is.raw(vals[[1L]][1L])) {
for (i in seq_along(vals)) {
vals[[i]][vals[[i]] == ri$nodata_value] <- NA ## hardcode to 1 for now
}
}
names(vals) <- sprintf("Band%i",band)
vals
}
#' type safe(r) raster read
#'
#' These wrappers around [vapour_read_raster()] guarantee single vector output of the nominated type.
#'
#' `*_hex` and `*_chr` are aliases of each other.
#' @inheritParams vapour_read_raster
#' @aliases vapour_read_raster_raw vapour_read_raster_int vapour_read_raster_dbl vapour_read_raster_chr vapour_read_raster_hex
#' @export
#' @return atomic vector of the nominated type raw, int, dbl, or character (hex)
#' @examples
#' f <- system.file("extdata", "sst.tif", package = "vapour")
#' vapour_read_raster_int(f, window = c(0, 0, 5, 4))
#' vapour_read_raster_raw(f, window = c(0, 0, 5, 4))
#' vapour_read_raster_chr(f, window = c(0, 0, 5, 4))
#' plot(vapour_read_raster_dbl(f, native = TRUE), pch = ".", ylim = c(273, 300))
vapour_read_raster_raw <- function(x, band = 1,
window,
resample = "nearestneighbour", ...,
sds = NULL, native = FALSE, set_na = TRUE, nara = FALSE) {
if (nara) message("nara incompatible with raw output, ignoring")
if (length(band) > 1) message("_raw output implies one band, using only the first")
vapour_read_raster(x, band = band, window = window, resample = resample, ..., sds = sds,
native = native, set_na = set_na, band_output_type = "Byte")[[1L]]
}
#' @name vapour_read_raster_raw
#' @export
vapour_read_raster_int <- function(x, band = 1,
window,
resample = "nearestneighbour", ...,
sds = NULL, native = FALSE, set_na = TRUE) {
if (length(band) > 1) message("_int output implies one band, using only the first")
vapour_read_raster(x, band = band, window = window, resample = resample, ..., sds = sds,
native = native, set_na = set_na, band_output_type = "Int32")[[1L]]
}
#' @name vapour_read_raster_raw
#' @export
vapour_read_raster_dbl <- function(x, band = 1,
window,
resample = "nearestneighbour", ...,
sds = NULL, native = FALSE, set_na = TRUE) {
if (length(band) > 1) message("_dbl output implies one band, using only the first")
vapour_read_raster(x, band = band, window = window, resample = resample, ..., sds = sds,
native = native, set_na = set_na, band_output_type = "Float64")[[1L]]
}
#' @name vapour_read_raster_raw
#' @export
vapour_read_raster_chr <- function(x, band = 1,
window,
resample = "nearestneighbour", ...,
sds = NULL, native = FALSE, set_na = TRUE) {
## band must be length 1, 3 or 4
if (length(band) == 2 || length(band) > 4) message("_chr output implies one, three or four bands ...")
if (length(band) == 2L) band <- band[1L]
if (length(band) > 4) band <- band[1:4]
bytes <- vapour_read_raster(x, band = band, window = window, resample = resample, ..., sds = sds,
native = native, set_na = set_na, band_output_type = "Byte")
## pack into character with as.raster ...
## note that we replicate out *3 if we only have one band ... (annoying of as.raster)
as.vector(grDevices::as.raster(array(unlist(bytes, use.names = FALSE), c(length(bytes[[1]]), 1, max(c(3, length(bytes)))))))
}
#' @name vapour_read_raster_raw
#' @export
vapour_read_raster_hex <- function(x, band = 1,
window,
resample = "nearestneighbour", ...,
sds = NULL, native = FALSE, set_na = TRUE) {
vapour_read_raster_chr(x, band = band, window = window, resample = resample, sds = sds,
native = native, set_na = set_na, ...)
}
#' Raster warper (reprojection)
#'
#' Read a window of data from a GDAL raster source through a warp specification.
#' The warp specification is provided by 'extent', 'dimension', and 'projection'
#' properties of the transformed output.
#'
#' Any bands may be read, including repeats.
#'
#' This function is not memory safe, the source is left on disk but the output
#' raster is all computed in memory so please be careful with very large values
#' for 'dimension'. `1000 * 1000 * 8` for 1000 columns, 1000 rows and floating
#' point double type will be 8Mb.
#'
#' There's control over the output type, and is auto-detected from the source
#' (raw/Byte, integer/Int32, numeric/Float64) or can be set with
#' 'band_output_type'.
#'
#' 'projection' refers to any projection string for a CRS understood by GDAL.
#' This includes the full Well-Known-Text specification of a coordinate
#' reference system, PROJ strings, "AUTH:CODE" types, and others. See
#' [vapour_srs_wkt()] for conversion from PROJ.4 string to WKT, and
#' [vapour_raster_info()] and [vapour_layer_info()] for various formats
#' available from a data source. Any string accepted by GDAL may be used for
#' 'projection' or 'source_projection', including EPSG strings, PROJ4 strings,
#' and file names. Note that this argument was named 'wkt' up until version
#' 0.8.0.
#'
#' 'extent' is the four-figure xmin,xmax,ymin,ymax outer corners of corner pixels
#'
#' 'dimension' is the pixel dimensions of the output, x (ncol) then y (nrow).
#'
#' Options for missing data are not yet handled, just returned as-is. Note that
#' there may be regions of "zero data" in a warped output, separate from
#' propagated missing "NODATA" values in the source.
#'
#' Argument 'source_projection' may be used to assign the projection of the
#' source, 'source_extent' to assign the extent of the source. Sometimes both
#' are required. Note, this is now better done by creating 'VRT', see [vapour_vrt()]
#' for assigning the source projection, extent, and some other options.
#'
#' If multiple sources are specified via 'x' and either 'source_projection' or
#' 'source_extent' are provided, these are applied to every source even if they
#' have valid values already. If this is not sensible please use VRT to wrap the
#' multiple sources first.
#'
#' Wild combinations of 'source_extent' and/or 'extent' may be used for
#' arbitrary flip orientations, scale and offset. For expert usage only. Old
#' versions allowed transform input for target and source but this is now
#' disabled (maybe we'll write a new wrapper for that).
#'
#' @section Options:
#'
#' The various options are convenience arguments for 'warp options -wo',
#' transformation options -to', 'open options -oo', and 'options' for any other
#' arguments in gdalwarp. There are no 'creation options -co' or 'dataset output
#' options -doo', because these are not supported by the MEM driver.
#'
#' All 'warp_options' are paired with a '-wo' declaration and similarly for '-to', and '-oo',
#' this is purely a convenience, since 'options' itself can be used for these as well but we recommend using
#' the individual arguments.
#' An example for warp options is `warp_options = c("SAMPLE_GRID=YES", "SAMPLE_STEPS=30")` and one for
#' general arguments might be
#' 'options = c("-ovr", "AUTO", "-nomd", "-cutline", "/path/to/cut.gpkg", "-crop_to_cutline")'. If they would
#' be separated by spaces on the command line then include as separate elements in the options character vector.
#'
#'
#' See [GDALWarpOptions](https://gdal.org/api/gdalwarp_cpp.html#_CPPv4N15GDALWarpOptions16papszWarpOptionsE) for '-wo'.
#'
#' See [GDAL transformation options](https://gdal.org/api/gdal_alg.html#_CPPv432GDALCreateGenImgProjTransformer212GDALDatasetH12GDALDatasetHPPc) for '-to'.
#'
#' See [GDALWARP command line app](https://gdal.org/programs/gdalwarp.html) for further details.
#'
#' Note we already apply the following gdalwarp arguments based on input R
#' arguments to this function.
#'
#' \describe{
#' \item{-of}{MEM is hardcoded, but may be extended in future}
#' \item{-t_srs}{set via 'projection'}
#' \item{-s_srs}{set via 'source_projection'}
#' \item{-te}{set via 'extent'}
#' \item{-ts}{set via 'dimension'}
#' \item{-r}{set via 'resample'}
#' \item{-ot}{set via 'band_output_type'}
#' \item{-te_srs}{ not supported}
#' \item{-a_ullr}{(not a gdalwarp argument, but we do analog) set via 'source_extent' use [vapour_vrt()] instead}
#' }
#'
#' In future all 'source_*' arguments may be deprecated in favour of
#' augmentation by 'vapour_vrt()'.
#'
#' Common inputs for `projection` are WKT variants, 'AUTH:CODE's e.g.
#' 'EPSG:3031', the 'OGC:CRS84' for lon,lat WGS84, 'ESRI:code' and other
#' authority variants, and datum names such as 'WGS84','NAD27' recognized by
#' PROJ itself.
#'
#' See help for 'SetFromUserInput' in 'OGRSpatialReference', and
#' 'proj_create_crs_to_crs'.
#'
#' [c.proj_create_crs_to_crs](https://proj.org/development/reference/functions.html#c.proj_create_crs_to_crs)
#'
#' [c.proj_create](https://proj.org/development/reference/functions.html#c.proj_create)
#'
#' [SetFromUserInput](https://gdal.org/doxygen/classOGRSpatialReference.html#aec3c6a49533fe457ddc763d699ff8796)
#'
#' @param x vector of data source names (file name or URL or database connection string)
#' @param bands index of band/s to read (1-based), may be new order or replicated, or NULL (all bands used, the default)
#' @param extent extent of the target warped raster 'c(xmin, xmax, ymin, ymax)'
#' @param source_extent extent of the source raster, used to override/augment incorrect source metadata
#' @param dimension dimensions in pixels of the warped raster (x, y)
#' @param projection projection of warped raster (in Well-Known-Text, or any projection string accepted by GDAL)
#' @param set_na NOT IMPLEMENTED logical, should 'NODATA' values be set to `NA`
#' @param resample resampling method used (see details in [vapour_read_raster])
#' @param source_projection optional, override or augment the projection of the source (in Well-Known-Text, or any projection string accepted by GDAL)
#' @param silent `TRUE` by default, set to `FALSE` to report messages
#' @param band_output_type numeric type of band to apply (else the native type if '') can be one of 'Byte', 'Int32', or 'Float64' but see details in [vapour_read_raster()]
#' @param ... unused
#' @param warp_options character vector of options, as in gdalwarp -wo - see Details
#' @param transformation_options character vector of options, as in gdalwarp -to see Details
#' @param open_options character vector of options, as in gdalwarp -oo - see Details
#' @param options character vectors of options as per the gdalwarp command line
#' @param nomd if `TRUE` the Metadata tag is removed from the resulting VRT (it can be quite substantial)
#' @param overview pick an integer overview from the source (0L is highest resolution, default -1L does nothing)
#' @param nara if 'TRUE' return in nativeRaster format
#' @export
#' @seealso vapour_read_raster vapour_read_raster_raw vapour_read_raster_int vapour_read_raster_dbl vapour_read_raster_chr vapour_read_raster_hex
#' @return list of vectors (only 1 for 'band') of numeric values, in raster order
#' @examples
#' b <- 4e5
#' f <- system.file("extdata", "sst.tif", package = "vapour")
#' prj <- "+proj=aeqd +lon_0=147 +lat_0=-42"
#' vals <- vapour_warp_raster(f, extent = c(-b, b, -b, b),
#' dimension = c(186, 298),
#' bands = 1,
#' projection = vapour_srs_wkt(prj),
#' warp_options = c("SAMPLE_GRID=YES"))
#'
#'
#' image(list(x = seq(-b, b, length.out = 187), y = seq(-b, b, length.out = 298),
#' z = matrix(unlist(vals, use.names = FALSE), 186)[,298:1]), asp = 1)
vapour_warp_raster <- function(x, bands = NULL,
extent = NULL,
dimension = NULL,
projection = "",
set_na = TRUE,
source_projection = NULL,
source_extent = 0.0,
resample = "near",
silent = TRUE, ...,
band_output_type = "",
warp_options = "",
transformation_options = "",
open_options = "",
options = "",
nomd = FALSE, overview = -1L, nara = FALSE) {
x <- .check_dsn_multiple(x)
if (!is.null(bands) && (anyNA(bands) || length(bands) < 1 || !is.numeric(bands))) {
stop("'bands' must be a valid set of band integers (1-based)")
}
if (is.null(projection) || is.na(projection[1]) || length(projection) != 1L || !is.character(projection)) {
stop("'projection' must be a valid character string for a GDAL/PROJ coordinate reference system (map projection)")
}
## deprecated arguments
if (!is.null(list(...)$source_wkt)) stop("'source_wkt' is defunct, please use 'source_projection'")
if(!is.null(list(...)$source_geotransform)) stop("'source_geotransform' is defunct and now ignored, please use 'source_extent'")
if (!is.null(list(...)$geotransform)) stop("'geotransform' is defunct and now ignored, used 'extent'")
if (band_output_type == "vrt") {
## special case, do nothing
} else {
band_output_type <- .r_to_gdal_datatype(band_output_type)
}
args <- list(...)
if (projection == "" && "wkt" %in% names(args)) {
projection <- args$wkt
message("please use 'projection = ' rather than 'wkt = ', use of 'wkt' is deprecated and will be removed in a future version")
}
## bands
if (is.numeric(bands) && any(bands < 1)) stop("all 'bands' index must be >= 1")
#if (is.null(bands)) bands <- 0
if (!is.null(bands)) {
if (length(bands) < 1 || anyNA(bands)) stop("'bands' must be numeric (integer), start at 1")
bands <- as.integer(bands)
}
##if ("projection" %in% names(list(...))) message("argument 'projection' input is ignored, warper functions use 'wkt = ' to specify target projection (any format is ok)")
# dud_extent <- FALSE
# if (is.null(extent)) {
# ## set a dummy and then pass in dud after following tests
# extent <- c(0, 1, 0, 1)
# dud_extent <- TRUE
# }
if(!is.numeric(extent)) {
if (isS4(extent)) {
extent <- c(extent@xmin, extent@xmax, extent@ymin, extent@ymax)
} else if (is.matrix(extent)) {
extent <- extent[c(1, 3, 2, 4)]
} else {
stop("'extent' must be numeric 'c(xmin, xmax, ymin, ymax)'")
}
}
if(!length(extent) == 4L) stop("'extent' must be of length 4")
if (any(diff(extent)[c(1, 3)] == 0)) stop("'extent' expected to be 'c(xmin, xmax, ymin, ymax)', zero x or y range not permitted")
if (length(source_extent) > 1 && any(diff(source_extent)[c(1, 3)] == 0)) stop("'extent' expected to be 'c(xmin, xmax, ymin, ymax)', zero x or y range not permitted")
if (!all(diff(extent)[c(1, 3)] > 0)) message("'extent' expected to be 'c(xmin, xmax, ymin, ymax)', negative values detected (ok for expert use)")
if (length(source_extent) > 1 && !all(diff(source_extent)[c(1, 3)] > 0)) message("'extent' expected to be 'c(xmin, xmax, ymin, ymax)', negative values detected (ok for expert use)")
## if (dud_extent) extent <- 0.0
## hmm, we can't rely on gdalwarp to give a sensibleish dimension if not specified, it goes for the native-res
#dud_dimension <- FALSE
## we dud it if no target projection is set, so you get native from the extent
if (!is.null(dimension)) {
## NO. We can't heuristic dimension or extent because we don't have a format to return those values with
## we make a simple raster, the image() thing and go with that
## FIXME: move this hardcode to C, and override with min(native_dimension, dimension)
#dimension <- c(512, 512)
## we could
## ## hardcode a default options(vapour.warp.default.dimension = c(512, 512))
## ## modify hardcode based on extent aspect ratio (not lat corrected)
## ## modify harcode to not exceed the native (I think I like this the best, because it reduces logic churn and delays when
## ## that has to be set in the C++, but we need to send down a message that the default is used (so do it all in C is the summ))
## set it to native with a max
## set it to native with a warn/override
#dud_dimension <- TRUE
#dimension <- c(2, 2)
if(!is.numeric(dimension)) stop("'dimension' must be numeric")
if(!length(dimension) == 2L) stop("'dimension must be of length 2'")
if(!all(dimension > 0)) stop("'dimension' values must be greater than 0")
if(!all(is.finite(dimension))) stop("'dimension' values must be finite and non-missing")
}
#if (dud_dimension) dimension <- 0
if (length(source_extent) > 1) {
if (!is.numeric(source_extent)) {
stop("'source_extent' must be numeric, of length 4 c(xmin, xmax, ymin, ymax)")
}
if (!all(is.finite(source_extent))) stop("'source_extent' values must be finite and non missing")
x <- vapour_vrt(x, source_extent = source_extent)
}
if(!is.null(source_projection)) {
if (!is.character(source_projection)) stop("source_projection must be character")
if(!silent) {
if(!nchar(source_projection) > 10) message("short 'source_projection', possibly invalid?")
}
x <- vapour_vrt(x, projection = source_projection)
}
if (!silent) {
if(!nchar(projection) > 0) message("target 'projection' not provided, read will occur from from source in native projection")
}
if (is.null(source_projection)) source_projection <- ""
resample <- tolower(resample[1L])
if (resample == "gauss") {
warning("Gauss resampling not available for warper, using NearestNeighbour")
resample <- "near"
}
rso <- c("near", "bilinear", "cubic", "cubicspline", "lanczos", "average",
"mode", "max", "min", "med", "q1", "q3", "sum") #, "rms")
if (!resample %in% rso) {
warning(sprintf("%s resampling not available for warper, using near", resample))
resample <- "near"
}
warp_options <- warp_options[!is.na(warp_options)]
#if (length(warp_options) < 1) warp_options <- ""
transformation_options <- transformation_options[!is.na(transformation_options)]
#if (length(transformation_options) < 1) transformation_options <- ""
open_options <- open_options[!is.na(open_options)]
# if (length(open_options) < 1) open_options <- ""
# dataset_output_options <- dataset_output_options[!is.na(dataset_output_options)]
## process all options into one big string list
if (any(grepl("-wo", options) |
grepl("-to", options) |
grepl("-oo", options))) {
##message("manually setting -wo, -to, -oo options detected, prefer use of 'warp_options', 'transformation_options', 'open_options'")
}
if (nchar(warp_options)[1L] > 0) options <- c(options, rbind("-wo", warp_options))
if (nchar(transformation_options)[1L] > 0) options <- c(rbind("-to", transformation_options))
if (nchar(open_options)[1L] > 0) options <- c(options, rbind("-oo", open_options))
#if (nchar(dataset_output_options)[1L] > 0) options <- c(rbind("-doo", dataset_output_options))
options <- options[!is.na(options)]
options <- options[nchar(options) > 0]
if (length(options) < 1) options <- ""
## no -r, -te, -t_srs, -ts, -of, -s_srs, -ot, -te_srs we set them manually
if (any(grepl("-r", options) |
grepl("-te", options) |
grepl("-t_srs", options) |
grepl("-ts", options) |
grepl("-of", options) |
grepl("-s_srs", options) |
grepl("-ot", options))) {
stop("manually setting -r, -te, -t_srs, -of, -s_srs, -ot options not allowed \n ( these controlled by arguments 'resample', 'target_extent', 'target_projection', '<MEM>', 'source_projection', 'band_output_type')")
}
if (any(grepl("-te_srs", options))) stop("setting '-te_srs' projection of target extent is not supported")
# vals <- warp_in_memory_gdal_cpp(x, source_WKT = source_projection,
# target_WKT = projection,
# target_extent = as.numeric(extent),
# target_dim = as.integer(dimension),
# bands = as.integer(bands),
# source_extent = as.numeric(source_extent),
# resample = resample,
# silent = silent,
# band_output_type = band_output_type,
# options = options, nomd = nomd, overview, nara = nara)
#
if (nara) {
vals <- gdal_raster_nara(x, target_crs = projection, target_ext = extent, target_dim = dimension,
bands = bands, resample = resample, band_output_type = band_output_type,
options = options, include_meta = !nomd)
} else {
## FIXME deal with nara and source projection source extent
if (is.null(bands)) bands <- 1L
vals <- gdal_raster_data(x, target_crs = projection, target_ext = extent, target_dim = dimension,
bands = bands, resample = resample, band_output_type = band_output_type,
options = options, include_meta = !nomd)
}
# ##// if we Dataset->RasterIO we don't have separated bands'
# nbands <- length(vals[[1L]]) / prod(as.integer(dimension))
# if (nbands > 1) vals <- split(vals[[1L]], rep(seq_len(nbands), each = prod(as.integer(dimension))))
#
# names(vals) <- make.names(sprintf("Band%i",seq_len(nbands)), unique = TRUE)
#
if (band_output_type == "vrt") return(vals[[1L]])
bands <- seq_along(vals)
names(vals) <- make.names(sprintf("Band%i",bands), unique = TRUE)
vals
}
#' type safe(r) raster warp
#'
#' These wrappers around [vapour_warp_raster()] guarantee single vector output of the nominated type.
#'
#' _hex and _chr are aliases of each other.
#' @inheritParams vapour_warp_raster
#' @aliases vapour_warp_raster_raw vapour_warp_raster_int vapour_warp_raster_dbl vapour_warp_raster_chr vapour_warp_raster_hex
#' @name vapour_warp_raster_raw
#' @export
#' @return atomic vector of the nominated type raw, int, dbl, or character (hex)
#' @examples
#' b <- 4e5
#' f <- system.file("extdata", "sst.tif", package = "vapour")
#' prj <- "+proj=aeqd +lon_0=147 +lat_0=-42"
#' bytes <- vapour_warp_raster_raw(f, extent = c(-b, b, -b, b),
#' dimension = c(18, 2),
#' bands = 1,
#' projection = prj)
#' # not useful given source type floating point, but works
#' str(bytes)
vapour_warp_raster_raw <- function(x, bands = NULL,
extent = NULL,
dimension = NULL,
projection = "",
set_na = TRUE,
source_projection = NULL,
source_extent = 0.0,
resample = "near",
silent = TRUE, ...,
warp_options = "",
transformation_options = "",
open_options = "",
options = "") {
if (length(bands) > 1 ) message("_raw output implies one band, ignoring all but the first")
vapour_warp_raster(x,
bands = bands,
extent = extent,
dimension = dimension,
projection = projection,
set_na = set_na,
source_projection = source_projection,
source_extent = source_extent,
resample = resample,
silent = silent,
band_output_type = "Byte",
warp_options = warp_options,
transformation_options = transformation_options,
open_options = open_options,
options = options,...)[[1L]]
}
vapour_warp_raster_vrt <- function(x, bands = NULL,
extent = NULL,
dimension = NULL,
projection = "",
set_na = TRUE,
source_projection = NULL,
source_extent = 0.0,
resample = "near",
silent = TRUE, ...,
warp_options = "",
transformation_options = "",
open_options = "",
options = "") {
vapour_warp_raster(x,
bands = bands,
extent = extent,
dimension = dimension,
projection = projection,
set_na = set_na,
source_projection = source_projection,
source_extent = source_extent,
resample = resample,
silent = silent,
band_output_type = "vrt",
warp_options = warp_options,
transformation_options = transformation_options,
open_options = open_options,
options = options,...)[[1L]]
}
#' @name vapour_warp_raster_raw
#' @export
vapour_warp_raster_int <- function(x, bands = NULL,
extent = NULL,
dimension = NULL,
projection = "",
set_na = TRUE,
source_projection = NULL,
source_extent = 0.0,
resample = "near",
silent = TRUE, ...,
warp_options = "",
transformation_options = "",
open_options = "",
options = "") {
if (length(bands) > 1 ) message("_int output implies one band, ignoring all but the first")
vapour_warp_raster(x,
bands = bands,
extent = extent,
dimension = dimension,
projection = projection,
set_na = set_na,
source_projection = source_projection,
source_extent = source_extent,
resample = resample,
silent = silent,
band_output_type = "Int32",
warp_options = warp_options,
transformation_options = transformation_options,
open_options = open_options,
options = options,...)[[1L]]
}
#' @name vapour_warp_raster_raw
#' @export
vapour_warp_raster_dbl <- function(x, bands = NULL,
extent = NULL,
dimension = NULL,
projection = "",
set_na = TRUE,
source_projection = NULL,
source_extent = 0.0,
resample = "near",
silent = TRUE, ...,
warp_options = "",
transformation_options = "",
open_options = "",
options = "") {
if (length(bands) > 1 ) message("_dbl output implies one band, ignoring all but the first")
vapour_warp_raster(x,
bands = bands,
extent = extent,
dimension = dimension,
projection = projection,
set_na = set_na,
source_projection = source_projection,
source_extent = source_extent,
resample = resample,
silent = silent,
band_output_type = "Float64",
warp_options = warp_options,
transformation_options = transformation_options,
open_options = open_options,
options = options,...)[[1L]]
}
#' @name vapour_warp_raster_raw
#' @export
vapour_warp_raster_chr <- function(x, bands = NULL,
extent = NULL,
dimension = NULL,
projection = "",
set_na = TRUE,
source_projection = NULL,
source_extent = 0.0,
resample = "near",
silent = TRUE, ...,
warp_options = "",
transformation_options = "",
open_options = "",
options = "") {
## band must be length 1, 3 or 4
if (length(bands) == 2 || length(bands) > 4) message("_chr output implies one, three or four bands ...")
if (length(bands) == 2L) bands <- bands[1L]
if (length(bands) > 4) bands <- bands[1:4]
bytes <- vapour_warp_raster(x,
bands = bands,
extent = extent,
dimension = dimension,
projection = projection,
set_na = set_na,
source_projection = source_projection,
source_extent = source_extent,
resample = resample,
silent = silent,
band_output_type = "Byte",
warp_options = warp_options,
transformation_options = transformation_options,
open_options = open_options,
options = options,...)
## note that we replicate out *3 if we only have one band ... (annoying of as.raster)
as.vector(grDevices::as.raster(array(unlist(bytes, use.names = FALSE), c(length(bytes[[1]]), 1, max(c(3, length(bytes)))))))
}
#' @name vapour_warp_raster_raw
#' @export
vapour_warp_raster_hex <- function(x, bands = NULL,
extent = NULL,
dimension = NULL,
projection = "",
set_na = TRUE,
source_projection = NULL,
source_extent = 0.0,
resample = "near",
silent = TRUE, ...,
warp_options = "",
transformation_options = "",
open_options = "",
options = "") {
vapour_warp_raster_chr(x,
bands = bands,
extent = extent,
dimension = dimension,
projection = projection,
set_na = set_na,
source_projection = source_projection,
source_extent = source_extent,
resample = resample,
silent = silent,
warp_options = warp_options,
transformation_options = transformation_options,
open_options = open_options,
options = options,...)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.