R: Create a GDAL virtual raster derived from one source dataset

rasterToVRT {gdalraster}

R Documentation

Create a GDAL virtual raster derived from one source dataset

Description

rasterToVRT() creates a virtual raster dataset (VRT format) derived from one source dataset with options for virtual subsetting, virtually resampling the source data at a different pixel resolution, or applying a virtual kernel filter. (See buildVRT() for virtual mosaicing.)

Usage

rasterToVRT(
  srcfile,
  relativeToVRT = FALSE,
  vrtfile = tempfile("tmprast", fileext = ".vrt"),
  resolution = NULL,
  subwindow = NULL,
  src_align = TRUE,
  resampling = "nearest",
  krnl = NULL,
  normalized = TRUE
)

Arguments

`srcfile`	Source raster filename.
`relativeToVRT`	Logical. Indicates whether the source filename should be interpreted as relative to the .vrt file (`TRUE`) or not relative to the .vrt file (`FALSE`, the default). If `TRUE`, the .vrt file is assumed to be in the same directory as `srcfile` and `basename(srcfile)` is used in the .vrt file. Use `TRUE` if the .vrt file will always be stored in the same directory with `srcfile`.
`vrtfile`	Output VRT filename.
`resolution`	A numeric vector of length two (xres, yres). The pixel size must be expressed in georeferenced units. Both must be positive values. The source pixel size is used if `resolution` is not specified.
`subwindow`	A numeric vector of length four (xmin, ymin, xmax, ymax). Selects `subwindow` of the source raster with corners given in georeferenced coordinates (in the source CRS). If not given, the upper left corner of the VRT will be the same as source, and the VRT extent will be the same or larger than source depending on `resolution`.
`src_align`	Logical. `TRUE`: the upper left corner of the VRT extent will be set to the upper left corner of the source pixel that contains `subwindow` xmin, ymax. The VRT will be pixel-aligned with source if the VRT `resolution` is the same as the source pixel size, otherwise VRT extent will be the minimum rectangle that contains `subwindow` for the given pixel size. Often, `src_align=TRUE` when selecting a raster minimum bounding box for a vector polygon. `FALSE`: the VRT upper left corner will be exactly `subwindow` xmin, ymax, and the VRT extent will be the minimum rectangle that contains `subwindow` for the given pixel size. If `subwindow` is not given, the source raster extent is used in which case `src_align=FALSE` has no effect. Use `src_align=FALSE` to pixel-align two rasters of different sizes, i.e., when the intent is target alignment.
`resampling`	The resampling method to use if xsize, ysize of the VRT is different than the size of the underlying source rectangle (in number of pixels). The values allowed are nearest, bilinear, cubic, cubicspline, lanczos, average and mode (as character).
`krnl`	A filtering kernel specified as pixel coefficients. `krnl` is a array with dimensions (size, size), where size must be an odd number. `krnl` can also be given as a vector with length size x size. For example, a 3x3 average filter is given by: krnl <- c( 0.11111, 0.11111, 0.11111, 0.11111, 0.11111, 0.11111, 0.11111, 0.11111, 0.11111) A kernel cannot be applied to sub-sampled or over-sampled data.
`normalized`	Logical. Indicates whether the kernel is normalized. Defaults to `TRUE`.

Details

rasterToVRT() can be used to virtually clip and pixel-align various raster layers with each other or in relation to vector polygon boundaries. It also supports VRT kernel filtering.

A VRT dataset is saved as a plain-text file with extension .vrt. This file contains a description of the dataset in an XML format. The description includes the source raster filename which can be a full path (relativeToVRT = FALSE) or relative path (relativeToVRT = TRUE). For relative path, rasterToVRT() assumes that the .vrt file will be in the same directory as the source file and uses basename(srcfile). The elements of the XML schema describe how the source data will be read, along with algorithms potentially applied and so forth. Documentation of the XML format for .vrt is at: https://gdal.org/drivers/raster/vrt.html.

Since .vrt is a small plain-text file it is fast to write and requires little storage space. Read performance is not degraded for certain simple operations (e.g., virtual clip without resampling). Reading will be slower for virtual resampling to a different pixel resolution or virtual kernel filtering since the operations are performed on-the-fly (but .vrt does not require the up front writing of a resampled or kernel-filtered raster to a regular format). VRT is sometimes useful as an intermediate raster in a series of processing steps, e.g., as a tempfile (the default).

GDAL VRT format has several capabilities and uses beyond those covered by rasterToVRT(). See the URL above for a full discussion.

Value

Returns the VRT filename invisibly.

Note

Pixel alignment is specified in terms of the source raster pixels (i.e., srcfile of the virtual raster). The use case in mind is virtually clipping a raster to the bounding box of a vector polygon and keeping pixels aligned with srcfile (src_align = TRUE). src_align would be set to FALSE if the intent is "target alignment". For example, if subwindow is the bounding box of another raster with a different layout, then also setting resolution to the pixel resolution of the target raster and src_align = FALSE will result in a virtual raster pixel-aligned with the target (i.e., pixels in the virtual raster are no longer aligned with its srcfile). Resampling defaults to nearest if not specified. Examples for both cases of src_align are given below.

rasterToVRT() assumes srcfile is a north-up raster.

Examples

## resample

evt_file <- system.file("extdata/storml_evt.tif", package="gdalraster")
ds <- new(GDALRaster, evt_file)
ds$res()
ds$bbox()
ds$close()

# table of the unique pixel values and their counts
tbl <- buildRAT(evt_file)
print(tbl)
sum(tbl$COUNT)

# resample at 90-m resolution
# EVT is thematic vegetation type so use a majority value
vrt_file <- rasterToVRT(evt_file,
                        resolution=c(90,90),
                        resampling="mode")

# .vrt is a small xml file pointing to the source raster
file.size(vrt_file)

tbl90m <- buildRAT(vrt_file)
print(tbl90m)
sum(tbl90m$COUNT)

ds <- new(GDALRaster, vrt_file)
ds$res()
ds$bbox()
ds$close()
vsi_unlink(vrt_file)


## clip

evt_file <- system.file("extdata/storml_evt.tif", package="gdalraster")
ds_evt <- new(GDALRaster, evt_file)
ds_evt$bbox()

# WKT string for a boundary within the EVT extent
bnd = "POLYGON ((324467.3 5104814.2, 323909.4 5104365.4, 323794.2
5103455.8, 324970.7 5102885.8, 326420.0 5103595.3, 326389.6 5104747.5,
325298.1 5104929.4, 325298.1 5104929.4, 324467.3 5104814.2))"

# src_align = TRUE
vrt_file <- rasterToVRT(evt_file,
                        subwindow = bbox_from_wkt(bnd),
                        src_align=TRUE)
ds_vrt <- new(GDALRaster, vrt_file)

# VRT is a virtual clip, pixel-aligned with the EVT raster
bbox_from_wkt(bnd)
ds_vrt$bbox()
ds_vrt$res()

ds_vrt$close()
vsi_unlink(vrt_file)

# src_align = FALSE
vrt_file <- rasterToVRT(evt_file,
                        subwindow = bbox_from_wkt(bnd),
                        src_align=FALSE)
ds_vrt_noalign <- new(GDALRaster, vrt_file)

# VRT upper left corner (xmin, ymax) is exactly bnd xmin, ymax
ds_vrt_noalign$bbox()
ds_vrt_noalign$res()

ds_vrt_noalign$close()
vsi_unlink(vrt_file)
ds_evt$close()


## subset and pixel align two rasters

# FARSITE landscape file for the Storm Lake area
lcp_file <- system.file("extdata/storm_lake.lcp", package="gdalraster")
ds_lcp <- new(GDALRaster, lcp_file)

# Landsat band 5 file covering the Storm Lake area
b5_file <- system.file("extdata/sr_b5_20200829.tif", package="gdalraster")
ds_b5 <- new(GDALRaster, b5_file)

ds_lcp$bbox()  # 323476.1 5101872.0  327766.1 5105082.0
ds_lcp$res()   # 30 30

ds_b5$bbox()   # 323400.9 5101815.8  327870.9 5105175.8
ds_b5$res()    # 30 30

# src_align = FALSE because we need target alignment in this case:
vrt_file <- rasterToVRT(b5_file,
                        resolution = ds_lcp$res(),
                        subwindow = ds_lcp$bbox(),
                        src_align = FALSE)
ds_b5vrt <- new(GDALRaster, vrt_file)

ds_b5vrt$bbox() # 323476.1 5101872.0  327766.1 5105082.0
ds_b5vrt$res()  # 30 30

# read the the Landsat file pixel-aligned with the LCP file
# summarize band 5 reflectance where FBFM = 165
# LCP band 4 contains FBFM (a classification of fuel beds):
ds_lcp$getMetadata(band=4, domain="")

# verify Landsat nodata (0):
ds_b5vrt$getNoDataValue(band=1)
# will be read as NA and omitted from stats
rs <- new(RunningStats, na_rm=TRUE)

ncols <- ds_lcp$getRasterXSize()
nrows <- ds_lcp$getRasterYSize()
for (row in 0:(nrows-1)) {
    row_fbfm <- ds_lcp$read(band=4, xoff=0, yoff=row,
                            xsize=ncols, ysize=1,
                            out_xsize=ncols, out_ysize=1)
    row_b5 <- ds_b5vrt$read(band=1, xoff=0, yoff=row,
                            xsize=ncols, ysize=1,
                            out_xsize=ncols, out_ysize=1)
	   rs$update(row_b5[row_fbfm == 165])
}
rs$get_count()
rs$get_mean()
rs$get_min()
rs$get_max()
rs$get_sum()
rs$get_var()
rs$get_sd()

ds_b5vrt$close()
vsi_unlink(vrt_file)
ds_lcp$close()
ds_b5$close()

[Package gdalraster version 1.11.1 Index]