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

Frank Davenport Fri, 25 Jul 2014 12:17:19 -0700

Hi Lydon,

Thanks for the work around! I can use that for now, but it is still notclear to me why the extract/mean works on each individual layer but noton the brick. When I've used extract() on bricks in the past it willreturn 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 individuallayer 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.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

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

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

Reply via email to