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
_______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo