[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:40:17 CEST 2014
Sorry for the mutliple postings but I found another solution using
raster::disaggregate(). Essentially the same as resample but prerserves
all the cell values. The fundamental issue was the small size of the
polygons compared to the cells (which can be problematic when using
weights()).
See the discussion here:
https://stackoverflow.com/questions/17766989/extract-data-from-raster-with-small-polygons-rounded-weights-too-small
g2<-raster::disaggregate(g,10)
test1<-extract(g2,dsd,fun=mean,na.rm=T,weights=T,small=T)
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