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
GDALRasterfor 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
rasterwill 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
interpis 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 bysrs_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)
#> id xcoord ycoord
#> 1 1 324650.9 5103344
#> 2 2 324171.0 5103034
#> 3 3 323533.4 5103329
#> 4 4 325220.0 5103508
#> 5 5 325703.1 5102377
#> 6 6 326297.8 5103924
#> 7 7 325520.4 5104146
#> 8 8 326247.7 5102506
#> 9 9 327711.7 5104476
#> 10 10 324181.7 5103901
raster_file <- system.file("extdata/storml_elev.tif", package="gdalraster")
pixel_extract(raster_file, pts[-1])
#> extracting from band 1...
#> 0...10...20...30...40...50...60...70...80...90...100 - done.
#> b1
#> [1,] 2648
#> [2,] 2865
#> [3,] 2717
#> [4,] 2560
#> [5,] 2916
#> [6,] 2633
#> [7,] 2548
#> [8,] 2801
#> [9,] 2475
#> [10,] 2822
# or as GDALRaster object
ds <- new(GDALRaster, raster_file)
pixel_extract(ds, pts[-1])
#> extracting from band 1...
#> 0...10...20...30...40...50...60...70...80...90...100 - done.
#> b1
#> [1,] 2648
#> [2,] 2865
#> [3,] 2717
#> [4,] 2560
#> [5,] 2916
#> [6,] 2633
#> [7,] 2548
#> [8,] 2801
#> [9,] 2475
#> [10,] 2822
# interpolated values
pixel_extract(raster_file, pts[-1], interp = "bilinear")
#> extracting from band 1...
#> 0...10...20...30...40...50...60...70...80...90...100 - done.
#> b1
#> [1,] 2649.217
#> [2,] 2881.799
#> [3,] 2716.290
#> [4,] 2558.797
#> [5,] 2920.404
#> [6,] 2629.495
#> [7,] 2548.250
#> [8,] 2810.543
#> [9,] 2478.609
#> [10,] 2819.776
# individual pixel values within a kernel
pixel_extract(raster_file, pts[-1], krnl_dim = 3)
#> extracting from band 1...
#> 0...10...20...30...40...50...60...70...80...90...100 - done.
#> b1_p1 b1_p2 b1_p3 b1_p4 b1_p5 b1_p6 b1_p7 b1_p8 b1_p9
#> [1,] 2660 2654 2651 2652 2648 2646 2653 2647 2645
#> [2,] 2924 2895 2901 2893 2865 2869 2863 2838 2841
#> [3,] NA 2709 2704 NA 2717 2717 NA 2725 2728
#> [4,] 2554 2558 2562 2557 2560 2564 2559 2562 2566
#> [5,] 2932 2920 2908 2927 2916 2903 2923 2911 2899
#> [6,] 2633 2642 2650 2627 2633 2640 2619 2624 2630
#> [7,] 2548 2548 2550 2549 2548 2550 2550 2548 2550
#> [8,] 2785 2777 2776 2807 2801 2803 2833 2825 2827
#> [9,] 2489 2482 2475 2479 2475 2470 2474 2471 2468
#> [10,] 2832 2810 2789 2839 2822 2800 2833 2811 2788
# 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")
#> extracting from band 1...
#> 0...10...20...30...40...50...60...70...80...90...100 - done.
#> b1
#> [1,] 2648
#> [2,] 2865
#> [3,] 2717
#> [4,] 2560
#> [5,] 2916
#> [6,] 2633
#> [7,] 2548
#> [8,] 2801
#> [9,] 2475
#> [10,] 2822
ds$close()