extract: Extract values from Raster objects

Description

Extract data from a Raster object using cell numbers or the locations of other spatial data (i.e., a spatial query). You can use coordinates (points), lines, polygons or an Extent (rectangle) object. This function has replaced five other functions in this package (cellValues, xyValues, lineValues, polygonValues, and focalValues)

Usage

extract(x, y, ...)

Arguments

Raster* object

A vector (representing cell numbers); or a SpatialPoints object or a two-column data.frame or matrix (representing points); or a SpatialPolygons, S

...

Additional arguments, see under Details

Value

A vector for RasterLayer objects, and a matrix for a RasterStack or RasterBrick object. A list if y is a SpatialPolygons* or SpatialLines* object or if a buffer argument is used (but not a fun argument). focal values are returned as a n column matrix. The first column has the column numbers, the remaining has the focal values for each layer.

Details

Below are additional arguments that can be used, depending on the type of objects used. = all objects = 1) fun. Function to summarize the values (e.g. mean). The function should take a single numeric vector as argument and return a single values (e.g. mean, min or max), and have a na.rm argument. If y represents points supplying a fun argument is only useful when non-zero buffers are used. 2) na.rm. Only useful when fun is supplied. If na.rm=TRUE (the default value), NA values are removed before fun is applied. This argument may be ignored if the function used has a ... argument and ignres an additional na.rm argument. = points = If y represents points, extract returns the values of a Raster* object for the cells in which a set of points fall. Additional arguments that can be suplied to this function if y represents points: 1) method. If method='simple' (the default), values for the cell a point falls in are returned. The alternative is method='bilinear', in which case the returned values are interpolated from the values of the four nearest raster cells. 2) buffer. The radius of a buffer around each point for which to extract cell values. It can be a single value, or a vector of the length of the number of points. If the data are not projected (latitude/longitude), the unit should be meters. Otherwise it should be in map-units (typically also meters). 3) cellnumbers. Logical. If cellnumbers=TRUE and buffer is not NULL, cellnumbers will also be returned. cellnumbers is set to FALSE when an argument fun is used. = polygons = If y represents polygons, the extract method returns the values of the cells of a Raster* object that is covered by a polygon. A cell is covered if its center is inside the polygon (but see the weights option for considering partly covered cells). Relatively small polygons (e.g. those smaller than a single cell of the RasterLayer), or polygons with an odd shape, may not return any values if they do not cover any RasterLayer cell centers. In some cases, you could use the centroids of such polygons, for example using extract(x, coordinates(y)) or extract(x, coordinates(y), method='bilinear'). Additional arguments that can be suplied to this function if y represents polygons are: If no fun argument is supplied: 1) weights. If TRUE, the function returns, for each polygon, a matrix with the cell values and the approximate fraction of each cell that is covered by the polygon(rounded to 1/100). The weights can be used for averaging; see examples. This option can be useful if the polygons are small relative to the cells size of the RasterLayer 2) cellnumbers. If cellnumbers=TRUE and weights=TRUE, cellnumbers will also be returned. = Focal values = If y is missing, you can supply two arguments to get 'focal' values (i.e. the values in a neighborhood around each cell) 1) row. row number (this is a required argument). 2) ngb. The neighborhood to consider. Either a single integer or a vector of two integers. Default is 3. See focal = All multi-layer objects = If x is a RasterStack or RasterBrick object, extract accepts, in addition to the arguments listed above, these two additional arguments: 1)layer. Integer. First layer for which you want values 2) nl. Integer. Number of layers for which you want values

Examples

Run this code

r <- raster(ncol=36, nrow=18)
r[] <- 1:ncell(r)

###############################
# extract values by cell number
###############################
extract(r, c(1:2, 10, 100))
s <- stack(r, sqrt(r), r/r)
extract(s, c(1, 10, 100), layer=2, n=2)

###############################
# extract values with points
###############################
xy <- cbind(-50, seq(-80, 80, by=20))
extract(r, xy)

sp <- SpatialPoints(xy)
extract(r, sp, method='bilinear')

# examples with a buffer
extract(r, xy[1:3,], buffer=1000000)
extract(r, xy[1:3,], buffer=1000000, fun=mean)

## illustrating the varying size of a buffer (expressed in meters) on a longitude/latitude raster
 z <- extract(r, xy, buffer=1000000)
 s <- raster(r)
 for (i in 1:length(z)) { s[z[[i]]] <- i }
## compare with raster that is not longitude/latitude
 projection(r) <- "+proj=UTM +zone=17" 
 xy[,1] <- 50
 z <- extract(r, xy, buffer=8)
 for (i in 1:length(z)) { s[z[[i]]] <- i }
 plot(s)
# library(maptools)
# data(wrld_simpl)
# plot(wrld_simpl, add=TRUE)

###############################
# extract values with lines
###############################

cds1 <- rbind(c(-50,0), c(0,60), c(40,5), c(15,-45), c(-10,-25))
cds2 <- rbind(c(80,20), c(140,60), c(160,0), c(140,-55))
lines <- SpatialLines(list(Lines(list(Line(cds1)), "1"), Lines(list(Line(cds2)), "2") ))

extract(r, lines)

###############################
# extract values with polygons
###############################
cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20))
cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))
polys <- SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1), Polygons(list(Polygon(cds2)), 2)))

#plot(r)
#plot(polys, add=TRUE)
v <- extract(r, polys)
v
# mean for each polygon
unlist(lapply(v, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA ))

# v <- extract(r, polys, weights=TRUE)
# v
## weighted mean for each polygon
# unlist(lapply(v, function(x) if (!is.null(x)) {sum(apply(x, 1, prod)) / sum(x[,2])} else NA  ))


###############################
# extract values with an extent
###############################
e <- extent(150,170,-60,-40)
extract(r, e)
#plot(r)
#plot(e, add=T)


###############################
# focal values
###############################
f <- extract(r, row=5, ngb=3)
f[1:15, ]