[R-sig-Geo] Extract with Large Rasters

Michael Treglia mtreglia at gmail.com
Wed Dec 17 00:02:53 CET 2014


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)

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]]



More information about the R-sig-Geo mailing list