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 <alobolis...@gmail.com> wrote: > Michael, > was this issue ever fully clarifed? > Thanks > Agus > > On Thu, Dec 18, 2014 at 7:54 PM, Michael Treglia <mtreg...@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.hijm...@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 <mdsum...@gmail.com> > >> wrote: > >> > Hello, see response inline. > >> > > >> > On Wed Dec 17 2014 at 10:03:23 Michael Treglia <mtreg...@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@data <- data.frame(points@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@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@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@r-project.org > >> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo > >> >> > >> > > >> > [[alternative HTML version deleted]] > >> > > >> > _______________________________________________ > >> > R-sig-Geo mailing list > >> > R-sig-Geo@r-project.org > >> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > >> > > > > [[alternative HTML version deleted]] > > > > _______________________________________________ > > R-sig-Geo mailing list > > R-sig-Geo@r-project.org > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo