Extract with Large Rasters
Hello, see response inline.
On Wed Dec 17 2014 at 10:03:23 Michael Treglia <mtreglia at gmail.com> wrote:
Hi All, I'm working with some pretty huge rasters and trying to extract data to SpatialPointsDataFrames, and running into some issues. I'm using the raster package for this. (sessionInfo at the end of this e-mail) Here are the specs of my raster: class : RasterLayer dimensions : 59901, 49494, 2964740094 (nrow, ncol, ncell) resolution : 9.332575, 9.332575 (x, y) extent : 222855.6, 684762, 3762092, 4321122 (xmin, xmax, ymin, ymax) coord. ref. : +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 data source : D:\GIS\Projects\ETynerensis\GISData\n35w092_n39_w096_UTM_cubic.tif names : n35w092_n39_w096_UTM_cubic values : 0, 837.6777 (min, max) When I run the following line, with a vector points layer in the same CRS as the raster, I get the error below.
points at data <- data.frame(points at data, extract(elev, points)
There's a few problems with your code, you should never use @ for getting and setting slots, and besides you can assign to the data with $ in the usual way (the developers do use @ in internal code to provide high level methods that behave in defined ways). I would do it like this: library(raster) ## simulate a raster, not large r <- raster(volcano, crs = "+proj=laea") ## more fake data pts0 <- xyFromCell(r, sample(ncell(r), 10), sp = TRUE) pts <- SpatialPointsDataFrame(pts, data.frame(id = 1:nrow(coordinates(pts0)))) ## rather than points at data <- etc. pts$r <- extract(r, pts) That might help if the overall points data is large, but you didn't include that information. You could use print(points) (with raster loaded) for a succinct summary. (Also points is a commonly used function so best avoided as a name). Cheers, Mike.
Error in .readCellsGDAL(x, uniquecells, layers) : NAs are not allowed in subscripted assignments In addition: Warning message: In .readCells(x, cells, 1) : NAs introduced by coercion
traceback()
8: .readCellsGDAL(x, uniquecells, layers) 7: .readCells(x, cells, 1) 6: .cellValues(object, cells, layer = layer, nl = nl) 5: .xyValues(x, coordinates(y), ..., df = df) 4: .local(x, y, ...) 3: extract(elev, sdata) 2: extract(elev, sdata) 1: data.frame(sdata at data, extract(elev, sdata)) I had a tough time trouble-shooting from the error, but the raster are huge (~2x10^9 cells), so I tried to clip a small portion of it and do the extract, and it worked. Here's the code I used to Crop:
bbox <- as(extent(22285.6, 3762091.8, 23285.6, 3862091.8),
'SpatialPolygons')
y <- crop(elev,bbox)
The data are FLT4S, and as a TIF, take up ~11.5GB. Is this just out of the bounds for what is feasible? Or is there some trouble-shooting anybody recommends I do? (Or is there a good way to do this as tiles in R?) And if there's a way to do this with spatial.tools, I'm open to that, but wasn't able to figure it out on my own thus far. Thanks for any suggestions, Mike T
sessionInfo()
R version 3.1.2 (2014-10-31)
Platform: x86_64-w64-mingw32/x64 (64-bit)
locale:
[1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United
States.1252 LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C LC_TIME=English_United
States.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] raster_2.2-31 rgdal_0.8-16 sp_1.0-15
loaded via a namespace (and not attached):
[1] grid_3.1.2 lattice_0.20-29 tools_3.1.2
[[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