[R-sig-Geo] Extract with Large Rasters

[Michael Sumner]

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
>

	[[alternative HTML version deleted]]