Re: [R-sig-Geo] raster::extract fails on brick but works on individual layers of brick
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
Re: [R-sig-Geo] raster::extract fails on brick but works on individual layers of brick
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
