Frank, I hope this issue has been solved in the development version. install.packages("raster", repos="http://R-Forge.R-project.org") I would appreciate feedback. Best, Robert
On Fri, Jul 25, 2014 at 12:40 PM, Frank Davenport <frank.davenp...@gmail.com> wrote: > 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.davenp...@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@r-project.org >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo > > > -- > Frank Davenport, Ph.D > Climate Hazards Group > UCSB Geography > > _______________________________________________ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo