[R-sig-Geo] Extract with Large Rasters
Agustin Lobo
alobolistas at gmail.com
Thu Jan 15 08:43:37 CET 2015
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
More information about the R-sig-Geo
mailing list