rasterToVRT: 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
)

Value

Returns the VRT filename invisibly.

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.

Examples

Run this code

## 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()