[R-sig-Geo] Extract with Large Rasters
Michael Treglia
mtreglia at gmail.com
Thu Dec 18 19:54:08 CET 2014
Thanks Mike and Robert,
Avoiding using @ didn't seem to help this problem. Also, I can do the
extract using the same exact points with different rasters (e.g., Bioclim
variables or a coarser DEM), so it doesn't seem to be a problem with the
points [Robert, if you still think they might pose an issue, let me know
and I'll send them off-list; I've included the specs of the points below
too though].
I realized I could clip the raster a bit and did so in with gdalwarp
(outside of R), and the extract worked without issue. I'm still curious to
figure out why it won't work with the original raster file though and I'm
all ears for suggestions, in case it can save me or others some trouble in
the future
Also, the reason I was using the @ was because I am ultimately doing this
for a stack of rasters, and would like to have all the raster fields added
to the SpatialPointsDataFrame. Is there a more preferred way to do that? I
see that if you do a hypothetical
>pts$field <- extract(stack, pts)
you get the values from the whole stack, but the columns have the prefix of
'field'. (I realize you can just rename columns after doing the extract)
Here's are the specs of the clipped Raster (for which the extract() works):
class : RasterLayer
dimensions : 35410, 39925, 1413744250 (nrow, ncol, ncell)
resolution : 9.332523, 9.332448 (x, y)
extent : 273952, 646553, 3885751, 4216213 (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:\GISData\Clipped_10mDEM.tif
names : Clipped_10mDEM
Here are the specs of the original (for which the extract() doesn't work):
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:\GISData\n35w092_n39_w096_UTM_cubic.tif
names : n35w092_n39_w096_UTM_cubic
Here are the Specs of the Points [excluding field names and max/min values]:
class : SpatialPointsDataFrame
features : 159
extent : 300852.6, 611979.9, 3910111, 4196814 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs
+ellps=WGS84 +towgs84=0,0,0
variables : 16
Thanks again!
Mike T
On Wed, Dec 17, 2014 at 2:06 PM, Robert J. Hijmans <r.hijmans at gmail.com>
wrote:
>
> 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
>
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list