pixel_extract: Extract pixel values at geospatial point locations

In gdalraster: Bindings to 'GDAL'

Extract pixel values at geospatial point locations

Description

pixel_extract() returns raster pixel values for a set of geospatial point locations. The coordinates are given as a two-column matrix of (x, y) values in the same spatial reference system as the input raster (unless xy_srs is specified). Values are extracted from all bands of the raster by default, or specific band numbers may be given. An optional interpolation method may be specified for bilinear (2 x 2 kernel), cubic convolution (4 x 4 kernel, GDAL >= 3.10), or cubic spline (4 x 4 kernel, GDAL >= 3.10). Alternatively, an optional kernel dimension may be given to extract values of the individual pixels within an N x N kernel centered on the pixel containing the point location. If xy_srs is given, the function will attempt to transform the input points to the projection of the raster with a call to transform_xy().

Usage

pixel_extract(
  raster,
  xy,
  bands = NULL,
  interp = NULL,
  krnl_dim = NULL,
  xy_srs = NULL,
  max_ram = 300
)

Arguments

raster	Either a character string giving the filename of a raster, or an object of class GDALRaster for the source dataset.
xy	A two-column numeric matrix or two-column data frame of geospatial (x, y) coordinates, or vector (x, y) for a single point, in the same spatial reference system as raster.
bands	Optional numeric vector of band numbers. All bands in raster will be processed by default if not specified.
interp	Optional character string specifying an interpolation method. Must be one of "bilinear", "cubic", "cubicspline", or "nearest" (the default if not specified, i.e., no interpolation). GDAL >= 3.10 is required for "cubic" and "cubicspline".
krnl_dim	Optional integer value specifying the dimension of an N x N kernel for which all individual pixel values will be returned. Only supported when extracting from a single raster band. Ignored if interp is specified as other than "nearest" (i.e., will always use the kernel implied by the interpolation method).
xy_srs	Optional character string specifying the spatial reference system for xy. May be in WKT format or any of the formats supported by srs_to_wkt().
max_ram	Numeric value giving the maximum amount of RAM (in MB) to use for potentially copying a remote raster into memory for processing (see Note). Defaults to 300 MB. Set to zero to disable potential copy of remote files into memory.

Value

A numeric matrix of pixel values with number of rows equal to the number of rows in xy, and number of columns equal to the number of bands, or if krnl_dim = N is used, number of columns equal to N * N. Named columns indicate the band number, e.g., "b1". If krnl_dim is used, named columns indicate band number and pixel, e.g., "b1_p1", "b1_p2", ..., "b1_p9" if krnl_dim = 3. Pixels are in left-to-right, top-to-bottom order in the kernel.

Note

Depending on the number of input points, extracting from a raster on a remote filesystem may require a large number of HTTP range requests which may be slow (i.e., URLs/remote VSI filesystems). In that case, it may be faster to copy the raster into memory first (either as MEM format or to a /vsimem filesystem). pixel_extract() will attempt to automate that process if the total size of file(s) that would be copied does not exceed the threshold given by max_ram, and length(xy) > 1 (requires GDAL >= 3.6).

For alternative workflows that involve copying to local storage, the data management functions (e.g., copyDatasetFiles()) and the VSI filesystem functions (e.g., vsi_is_local(), vsi_stat(), vsi_copy_file()) may be of interest.

Examples

pt_file <- system.file("extdata/storml_pts.csv", package="gdalraster")
# id, x, y in NAD83 / UTM zone 12N, same as the raster
pts <- read.csv(pt_file)
print(pts)

raster_file <- system.file("extdata/storml_elev.tif", package="gdalraster")

pixel_extract(raster_file, pts[-1])

# or as GDALRaster object
ds <- new(GDALRaster, raster_file)
pixel_extract(ds, pts[-1])

# interpolated values
pixel_extract(raster_file, pts[-1], interp = "bilinear")

# individual pixel values within a kernel
pixel_extract(raster_file, pts[-1], krnl_dim = 3)

# lont/lat xy
pts_wgs84 <- transform_xy(pts[-1], srs_from = ds$getProjection(),
                          srs_to = "WGS84")

# transform the input xy
pixel_extract(ds, xy = pts_wgs84, xy_srs = "WGS84")

ds$close()

gdalraster documentation built on June 8, 2025, 12:37 p.m.