Spencer,

That is an interesting challenge. Your approaches are fine in principle,
but you have a raster of 30 meter resolution for the conterminous US ( ~17
billion cells ) and, if I understand you well, you want to tabulate the
values by polygon, using irregular polygons (relative to the raster) of
0.5x0.5 degrees.  That is a good case study for me to improve the speed of
rasterize (which is used in both of your methods as extract will call it
for each polygon); and I might work on that this summer.

For now, I have not much help to offer, except that you can probably make
things more efficient by removing the polygons that are outside the US:

#load appropriate libraries
library(rgdal)
library(raster)
library(rgeos)

r <- raster(xmn=-124.5, xmx=-67, ymn=25, ymx=49, nrow=48, ncol=115)
r[] <- 1:ncell(r)
p <- as(r, 'SpatialPolygonsDataFrame')

# remove the cells that are outside of the USA
# that is 35% and could saves a lot of processing time

usa <- getData('GADM', country='USA', level=0)
i <- gIntersects(p, usa, byid=TRUE)  #some patience required here
pp <- p[as.vector(i), ]

#  projection (changed to the exact spec of the landcover data)
prj <- CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0
+y_0=0 +ellps=GRS80 +units=m")
us <- spTransform(pp, prj)

#load land cover data raster (16 GB)
env <- raster('nlcd2006_landcover_4-20-11_se5.img', RAT=F)

plot(env) #remakably fast!
plot(us, add=TRUE, border='red')

NAvalue(env) <- 0

rasterize/crosstab is probably faster (according to you numbers), but if
you are going the 'extract' route, you need to do some looping, otherwise
you will attempt to create a list with 10 billion values, and that will
likely not work. Something like this, which is indeed rather slow:

res <- list()
for (i in 1:nrow(us)) {
print( system.time( ex <- extract(env, us[i,]) ) )
res[i] <- table(ex[[1]])
cat(i, '\n')
flush.console()
}

# and now process res, but that should be easy

Best, Robert

On Mon, Jun 18, 2012 at 4:07 PM, Spencer Scheidt <scheidt14@gmail.com>wrote:

> Hello,
>
> I am attempting to extract land cover values from the NLCD 2006 dataset
> (http://www.mrlc.gov/nlcd2006.php) using a grid that I have converted
> into a
> SpatialPolygons* object.  However, the NLCD raster is approximately 16 GB
> in
> size (30 m resolution), and the grid that I would like to work with
> contains
> 5520 polygons (0.5 degree resolution).
>
> I have written a script based on previous postings that have asked a
> similar
> question and Applied Spatial Data Analysis in R by Bivand et. al, and used
> this approach successfully for smaller SpatialPolygons* objects with the
> same NLCD dataset.
>
> What I am wondering is if there is a quicker way to extract my polygon
> values, either by breaking the raster and polygons into smaller chunks, or
> using an entirely different approach.  I realize that the extreme sizes of
> my 2 objects may be the ultimate constraint, and if necessary I can reduce
> the resolution of my grid.
>
> Below is my working code, which uses 2 different approaches, via rasterize
> (estimated to take 10+ days to run and cross tabulate) and extract
> (estimated to take over 30 days to run):
>
> #load appropriate libraries
> library(maps)
> library(rgdal)
> library(raster)
> library(sp)
> library(maptools)
>
> #create bounding box
> us.box<-matrix(c(-124.5,-67,25,49), nrow = 2, ncol = 2, byrow = TRUE)
>
> #create grid
> getClass("GridTopology")
>
> cs<-c(0.5, 0.5)
> cc<-us.box[,1] + (cs/2)
> cd<-ceiling(diff(t(us.box))/cs)
> us.grid<-GridTopology(cellcentre.offset = cc, cellsize = cs, cells.dim =
> cd)
> us.grid
>
> #create spatial polygons
> #projection is albers equal area conic
> prj.string <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96
> +units=m"
> us.sp <- as.SpatialPolygons.GridTopology(us.grid, proj4string =
> CRS("+proj=longlat +ellps=WGS84"))
> us.SP<-spTransform(us.sp, CRS(prj.string))
>
> summary(us.SP)
>
> #check to make sure ID values are appropriate - should be listed by row
> from
> top left
> IDvaluesGridTopology(us.grid)
> slot(us.SP, "polygons")
>
> #load land cover data raster (16 GB)
> env.data = raster('nlcd2006_landcover_4-20-11_se5.img')
> projection(env.data) = CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23
> +lon_0=-96")
>
> #create raster of polygons
> poly.raster <- rasterize(us.SP, env.data, progress = 'window')
> df <- crosstab(poly.raster, env.data, progress = 'window')
>
> #maybe another way?
> ex.test <- extract(env.data, us.SP, df=TRUE, progress = 'window')
>
> Any help or advice on working with large datasets is much appreciated.
> Thank you for your time.
>
> Spencer
>
>
>
> --
> View this message in context:
> http://r-sig-geo.2731867.n2.nabble.com/Extract-polygon-values-from-large-raster-dataset-tp7580250.html
> Sent from the R-sig-geo mailing list archive at Nabble.com.
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

	[[alternative HTML version deleted]]

