Extract with Large Rasters

From the traceback

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

        [[alternative HTML version deleted]]