Error when extracting values from a rasterBrick or rasterStack

1.3k views Asked by At

I am having trouble extracting values or a point from a multi band raster of class rasterStack or rasterBrick. 'extract' works well with the individual rasters but posts an error when applied to the rasterStack or brick.

> all.var
class       : RasterBrick 
dimensions  : 89, 180, 16020, 34  (nrow, ncol, ncell, nlayers)
resolution  : 2, 2  (x, y)
extent      : -179, 181, -89, 89  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       :    period_max, pct_signif_periodmax, pct_coi_periodmax,    pct_ispos_signif, events_pos_periodmax, events_neg_periodmax, events_pos_all, events_neg_all,  maxpower_pos,  maxpower_neg, maxpower_events_pos, maxpower_events_neg, maxpower_pos_norm, maxpower_neg_norm, maxpower_events_pos_norm, ... 

> point
Lon      Lat
1  166.2790 -10.2690
2   26.9000 -33.6000
3  153.6209 -28.7001
4  113.8333 -28.6833
5  153.6335 -28.6591
6  153.5836 -28.4643
7   32.6833 -27.5333
8   32.6880 -27.5260
9   32.6880 -27.5260
10  32.6880 -27.5260


> point.extract<-extract(all.var, point, buffer=50000,na.rm=TRUE,fun=mean)
Error in apply(x, 2, fun2) : dim(X) must have a positive length

This works with individual rasters but fails with stack/brick and elicits an error only when I use a buffer argument.

1

There are 1 answers

15
Robert Hijmans On

Here is a working R example that illustrates the error:

library(raster)
b <- brick(nrow=89, ncol=180, nl=34, xmn=-179, xmx=181, ymn=-89, ymx=89, crs="+proj=longlat +datum=WGS84")
b[] <- 1

p <- matrix(c(166.2790,-10.2690,26.9000,-33.6000,153.6209,-28.7001,113.8333,-28.6833,153.6335,-28.6591,153.5836,-28.4643,32.6833,-27.5333,32.6880,-27.5260,32.6880,-27.5260,32.6880,-27.5260), ncol=2, byrow=TRUE)

v <- extract(b, p, buffer=50000, na.rm=TRUE, fun=mean)

That indeed gives the error you reported, probably due a bug in the raster package. Here is a work-around:

v <- extract(b, p, buffer=15000000)
# get the mean for each point (buffer) by layer
vv <- lapply(v, function(x) ifelse(is.matrix(x), colMeans(x, na.rm=TRUE), x))
# combine
do.call(rbind, vv)