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). Coordinates can also be given in a data frame with
an optional column of point IDs.
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 N
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,
as_data_frame = NULL
)
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 coordinates (x, y), or a vector for a single point (x, y), in the same spatial reference system as
raster
. Can also be given as a data frame with three or more columns, in which the first column must be a vector of point IDs (character
ornumeric
), and the second and third columns must contain the geospatial coordinates (x, y).- bands
Optional numeric vector of band numbers. All bands in
raster
will be processed by default if not specified, or if0
is given.- 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 numeric value specifying the dimension
N
pixels of anN x N
kernel for which all individual pixel values will be returned. Should be a positive whole number (will be coerced to integer by truncation). Currently only supported when extracting from a single raster band. Ignored ifinterp
is specified as other than"nearest"
(i.e., the kernel implied by the interpolation method will always be used).- 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.
- as_data_frame
Logical value,
TRUE
to return output as a data frame. The default is to return a numeric matrix unless point IDs are present in the first column ofxy
given as a data frame. In that latter case, the output will always be a data frame with the point IDs in the first column.
Value
A numeric matrix or data frame of pixel values with number of rows
equal to the number of rows in xy
. The number of columns is equal to the
number of bands
(plus optional point ID column), or if krnl_dim = N
is used, number of columns is equal to N * N
(plus optional point ID
column). Output is always a data frame if xy
is given as a data frame with
a column of point IDs (as_data_frame
will be ignored in that case).
Named columns indicate the band, e.g., "b1"
. If krnl_dim
is used, named
columns indicate band 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 nrow(xy) > 10
(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)
#> extracting from band 1...
#> 0...10...20...30...40...50...60...70...80...90...100 - done.
#> id storml_elev
#> 1 1 2648
#> 2 2 2865
#> 3 3 2717
#> 4 4 2560
#> 5 5 2916
#> 6 6 2633
#> 7 7 2548
#> 8 8 2801
#> 9 9 2475
#> 10 10 2822
# or as GDALRaster object
ds <- new(GDALRaster, raster_file)
pixel_extract(ds, pts)
#> extracting from band 1...
#> 0...10...20...30...40...50...60...70...80...90...100 - done.
#> id storml_elev
#> 1 1 2648
#> 2 2 2865
#> 3 3 2717
#> 4 4 2560
#> 5 5 2916
#> 6 6 2633
#> 7 7 2548
#> 8 8 2801
#> 9 9 2475
#> 10 10 2822
# interpolated values
pixel_extract(raster_file, pts, interp = "bilinear")
#> extracting from band 1...
#> 0...10...20...30...40...50...60...70...80...90...100 - done.
#> id storml_elev
#> 1 1 2649.217
#> 2 2 2881.799
#> 3 3 2716.290
#> 4 4 2558.797
#> 5 5 2920.404
#> 6 6 2629.495
#> 7 7 2548.250
#> 8 8 2810.543
#> 9 9 2478.609
#> 10 10 2819.776
# individual pixel values within a kernel
pixel_extract(raster_file, pts, krnl_dim = 3)
#> extracting from band 1...
#> 0...10...20...30...40...50...60...70...80...90...100 - done.
#> id b1_p1 b1_p2 b1_p3 b1_p4 b1_p5 b1_p6 b1_p7 b1_p8 b1_p9
#> 1 1 2660 2654 2651 2652 2648 2646 2653 2647 2645
#> 2 2 2924 2895 2901 2893 2865 2869 2863 2838 2841
#> 3 3 NA 2709 2704 NA 2717 2717 NA 2725 2728
#> 4 4 2554 2558 2562 2557 2560 2564 2559 2562 2566
#> 5 5 2932 2920 2908 2927 2916 2903 2923 2911 2899
#> 6 6 2633 2642 2650 2627 2633 2640 2619 2624 2630
#> 7 7 2548 2548 2550 2549 2548 2550 2550 2548 2550
#> 8 8 2785 2777 2776 2807 2801 2803 2833 2825 2827
#> 9 9 2489 2482 2475 2479 2475 2470 2474 2471 2468
#> 10 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.
#> storml_elev
#> [1,] 2648
#> [2,] 2865
#> [3,] 2717
#> [4,] 2560
#> [5,] 2916
#> [6,] 2633
#> [7,] 2548
#> [8,] 2801
#> [9,] 2475
#> [10,] 2822
ds$close()