[R-sig-Geo] raster::extract fails on brick but works on individual layers of brick
Robert J. Hijmans
r.hijmans at gmail.com
Thu Sep 11 18:04:57 CEST 2014
Frank,
Thanks for your feedback. I indeed added some code to get the weights
adding up to 1 in more (all) cases. The difference is proportional and
the results should otherwise be the same.
Robert
On Wed, Sep 10, 2014 at 3:22 PM, Frank Davenport
<frank.davenport at gmail.com> wrote:
> Hi Robert-
>
> I was unable to install the development version, however I did update to
> raster- 2.3-0 and the issue appears to be resolved. I can run
> extract(..,..,mean) with or without weights and get values each time,
> without having to use dissagregate first.
>
> However I noticed now that the 'weights' colunm sums to one. This what I
> would like from a weights column however this is not the how weights behaved
> in the past, as per this discussion:
>
> https://stackoverflow.com/questions/17766989/extract-data-from-raster-with-small-polygons-rounded-weights-too-small
>
>
> My new question: Has the weights function fundamentally changed or is it
> different only when the polygons are significantly smaller than the raster
> cells? Or put another way, if I re-run old code that uses weights=T, should
> I expect different answers?
>
> Thanks again for the great package!
>
> Frank
>
> Here is my code and session info:
>
>> library(raster)
> Loading required package: sp
> Warning message:
> package ‘raster’ was built under R version 3.1.1
>> #install.packages("raster", repos="http://R-Forge.R-project.org")
>>
>> load('~/Dropbox/Public/99_raster_bugreport.Rdata')
>> #https://dl.dropboxusercontent.com/u/9577903/99_raster_bugreport.Rdata
>> test1<-extract(g,dsd,fun=mean,na.rm=T,weights=T) #No Error
>> test2<-extract(g,dsd,fun=mean,na.rm=T,weights=F) #No Error
>>
>> test3<-extract(g,dsd,weights=T) #just get the weights
>>
>> test3[[1]][,'weight'] #the weights column sums to 1
> [1] 0.6 0.4
>
>>
>> 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.3-0 sp_1.0-15
>
> loaded via a namespace (and not attached):
> [1] grid_3.1.0 lattice_0.20-29 tools_3.1.0
> On 9/2/14 9:41 AM, Robert J. Hijmans wrote:
>>
>> 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.davenport at 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.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
>>>
>>> _______________________________________________
>>> 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