raster::extract fails on brick but works on individual layers of brick
Hi Lydon, Thanks for the work around! I can use that for now, but it is still not clear to me why the extract/mean works on each individual layer but not on the brick. When I've used extract() on bricks in the past it will return a mean (or whatever function) value for each layer of the brick (but much faster than doing it individually). In this case, since I can use extract(..,..,mean) on each individual layer without error, which led me to believe this may be a bug. Regardless, I appreciate the help. Thanks again, F
On 7/23/14 10:07 AM, Lyndon Estes wrote:
Hi Frank, I think the mean function is getting in the way, since if you want to extra the values for all cells each polygon overlaps, the outputs will first end up in a list. test1 <- extract(g, dsd, weights = TRUE, small = TRUE) Will get the cell values for each polygon from each layer in the brick, along with their weight (while allowing very small polygons to get their underlying cell values). I would then process the mean function on the list. v <- t(sapply(test1, function(x) apply(t(sapply(1:nrow(x), function(y) x[y, 1:10] * x[y, 11])), 2, sum))) # matrix of annual values per polygon Hope this helps. Cheers, Lyndon On Wed, Jul 23, 2014 at 12:17 PM, Frank Davenport <frank.davenport at gmail.com> wrote:
Hello,
I am using extract() to aggregate values from a raster to a polygon. The
raster is a brick object. Extract fails on the brick object but succeeds
when applied on each individual layer of the brick (via a for() loop). I
would use this as a work around but in my actual scenario the brick is much
bigger and the individual approach is too slow.
The specific error message is: "Error in t(sapply(res, meanfunc)) : error
in evaluating the argument 'x' in selecting a method for function 't':
Error in apply(x, 1, function(X) { : dim(X) must have a positive length"
I believe this has something to do with the relative sizes of the polygons
(small) and the grid cells. I have successfully extracted this brick to
another shapefile that had much larger polygons. Likewise I can also do the
extraction if I resample the cells to a smaller resolution (not shown
below).
Finally the extract will work if I set 'na.rm=F' but then produces mostly
NA's, even though there are no NAs in the brick. I realize this might be
something to do with the dataType(). However that does not explain why
extract works on individual layers, but not the whole brick.
Reproducible code is below and prettier example can be found here:
http://rpubs.com/frank_davenport/22698
The data necessary to run the example can be found here:
https://dl.dropboxusercontent.com/u/9577903/99_raster_bugreport.Rdata
Thank you for your help and apologies if a solution is already posted.
Frank
rm(list=ls()) library(raster)
Loading required package: sp
##-Download data from here:
##contains 'dsd' a spatial polygons data.frame and 'g' a raster brick
with 10 layers
load('~/Dropbox/Public/99_raster_bugreport.Rdata')
dsd
class : SpatialPolygonsDataFrame features : 36 extent : 34.04665, 39.70319, -4.67742, 1.19921 (xmin, xmax, ymin, ymax) coord. ref. : +proj=longlat +a=6378249.145 +b=6356514.96582849 +no_defs variables : 4 names : Province, District, Division, uidu min values : CENTRAL, BUNGOMA, ABOTHUGUCHI WEST, 1 max values : WESTERN, VIHIGA, WINAM, 44
g
class : RasterBrick dimensions : 24, 20, 480, 10 (nrow, ncol, ncell, nlayers) resolution : 0.5, 0.5 (x, y) extent : 33, 43, -6, 6 (xmin, xmax, ymin, ymax) coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 data source : in memory names : X1999.03.01, X1999.03.02, X1999.03.03, X1999.03.04, X1999.03.05, X1999.03.06, X1999.03.07, X1999.03.08, X1999.03.09, X1999.03.10 min values : 22.47562, 22.44415, 22.25507, 22.23166, 22.12387, 22.42477, 22.27802, 22.15134, 22.36218, 22.33447 max values : 41.47818, 40.53116, 41.54944, 42.33093, 41.67810, 40.79260, 41.83319, 40.83359, 41.12604, 42.00555
#---Show the data plot(g[[1]]) plot(dsd,add=T) #--Try to extract, with weights yeilds an error test1<-extract(g,dsd,fun=mean,na.rm=T,weights=T)
Error in t(sapply(res, meanfunc)) :
error in evaluating the argument 'x' in selecting a method for function
't': Error in apply(x, 1, function(X) { : dim(X) must have a positive length
#--Extract without weights--produces NA's for most polygons test2<-extract(g,dsd,fun=mean,na.m=T) head(test2)
X1999.03.01 X1999.03.02 X1999.03.03 X1999.03.04 X1999.03.05 X1999.03.06 X1999.03.07 X1999.03.08 X1999.03.09 X1999.03.10 [1,] NA NA NA NA NA NA NA NA NA NA [2,] NA NA NA NA NA NA NA NA NA NA [3,] NA NA NA NA NA NA NA NA NA NA [4,] NA NA NA NA NA NA NA NA NA NA [5,] NA NA NA NA NA NA NA NA NA NA [6,] NA NA NA NA NA NA NA NA NA NA
summary(test2)
X1999.03.01 X1999.03.02 X1999.03.03 X1999.03.04 X1999.03.05 X1999.03.06 X1999.03.07 X1999.03.08 Min. :24.47 Min. :25.21 Min. :24.67 Min. :25.47 Min. :27.42 Min. :25.38 Min. :25.93 Min. :24.11 1st Qu.:29.29 1st Qu.:30.99 1st Qu.:29.55 1st Qu.:30.83 1st Qu.:31.58 1st Qu.:29.78 1st Qu.:29.51 1st Qu.:28.35 Median :30.31 Median :33.63 Median :31.13 Median :31.88 Median :33.32 Median :30.31 Median :30.52 Median :29.71 Mean :29.87 Mean :32.75 Mean :30.41 Mean :31.78 Mean :32.61 Mean :29.71 Mean :29.90 Mean :29.19 3rd Qu.:31.86 3rd Qu.:36.08 3rd Qu.:32.59 3rd Qu.:34.37 3rd Qu.:34.73 3rd Qu.:30.60 3rd Qu.:31.32 3rd Qu.:30.82 Max. :32.81 Max. :37.00 Max. :33.45 Max. :35.77 Max. :35.42 Max. :31.98 Max. :31.69 Max. :32.53 NA's :30 NA's :30 NA's :30 NA's :30 NA's :30 NA's :30 NA's :30 NA's :30 X1999.03.09 X1999.03.10 Min. :25.98 Min. :25.53 1st Qu.:29.53 1st Qu.:30.73 Median :30.57 Median :31.06 Mean :30.12 Mean :31.24 3rd Qu.:31.41 3rd Qu.:33.52 Max. :32.75 Max. :34.83 NA's :30 NA's :30
#--Extract each layer seperatley--works but is SLOW
glist<-vector(mode='list',length=nlayers(g))
for(i in 1:length(glist)){
+ glist[[i]]<-extract(g[[i]],dsd,fun=mean,na.rm=T,weigths=T) + }
sessionInfo()
R version 3.1.0 (2014-04-10)
Platform: x86_64-apple-darwin10.8.0 (64-bit)
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] raster_2.2-31 sp_1.0-15
loaded via a namespace (and not attached):
[1] grid_3.1.0 lattice_0.20-29 tools_3.1.0
[[alternative HTML version deleted]]
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Frank Davenport, Ph.D Climate Hazards Group UCSB Geography