[R-sig-Geo] Extract with Large Rasters
Michael Treglia
mtreglia at gmail.com
Thu Jan 15 16:03:40 CET 2015
Hi Agustin,
Nope - since the process that I needed worked after clipping the rasters a
bit, I just did that, though didn't come to an actual solution that enabled
me to use my original, larger rasters referenced in my earlier e-mails.
Cheers,
Mike
On Thu, Jan 15, 2015 at 1:43 AM, Agustin Lobo <alobolistas at gmail.com> wrote:
> Michael,
> was this issue ever fully clarifed?
> Thanks
> Agus
>
> On Thu, Dec 18, 2014 at 7:54 PM, Michael Treglia <mtreglia at gmail.com>
> wrote:
> > 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]]
> >
> > _______________________________________________
> > 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