[R-sig-Geo] Extract with Large Rasters

[Robert J. Hijmans]

>From the traceback

8: .readCellsGDAL(x, uniquecells, layers)

it seems to me that the error occurs here:

extract(elev, points)

and that it has nothing to do with the size of the raster.  I am not
sure what causes it.
Could you perhaps email me the points?

Thanks, Robert


On Wed, Dec 17, 2014 at 2:27 AM, Michael Sumner <mdsumner at gmail.com> wrote:
> 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]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo