[R-sig-Geo] raster::extract fails on brick but works on individual layers of brick

Frank Davenport frank.davenport at gmail.com
Fri Jul 25 21:16:06 CEST 2014


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:
>> https://dl.dropboxusercontent.com/u/9577903/99_raster_bugreport.Rdata
>>> ##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



More information about the R-sig-Geo mailing list