[R-sig-Geo] "subset" a large raster object

[Robert J. Hijmans]

Don,

I am guessing that you did not provide a function to rasterToPoints to
subset the values. That should work if the number of cells within the
range specified is not extremely large.

r = raster(filename, band=1)
p = rasterToPoints(r, fun=function(x){x > 21 & x < 25})


If you need to combine various bands, you could first do some raster
processing as in:

r = raster(filename, band=1)
s = r > 21 & r < 25
s[!s] = NA

#or
s = reclass(r, c(-Inf, 11, NA, 22, 24, 1, 25, Inf, NA))
# 1/NA result

r = raster(filename, band=2)
s2 = reclass(r, c(-Inf, 11, NA, 22, 24, 1, 25, Inf, NA))

etc

s = s * s2 * ...

all of that should worn on a file of (practically) any size.
then

p = pointsToRaster(s)

Robert


On Wed, Jan 27, 2010 at 2:03 PM, Don MacQueen <macq at llnl.gov> wrote:
> I have a SpatialGridDataFrame named tmp.lc [see below for part of the str()
> information on it] obtained by using readGDAL() on an image file.
>
> I would like to subset it in the sense of finding the portions of the raster
> image that meet a defined criterion, so that I can use, for example, chull()
> to find a polygon that surrounds that portion of the raster.
>
> For a relatively small object (150x150) I can do the job by converting to
> SpatialPointsDataFrame. Then subsetting is easly. See example below.
>
> But when I go to the full 5100x5400 image the method fails due to
> insufficient memory. Similarly if I use rasterToPoints() from the raster
> package from R-forge.
>
> I would appreciate suggestions for an alternative.
> (perhaps GRASS, to which I have just a very small amount of exposure?)
>
> Thanks
> -Don
>
>
> ## start with object tmp.lc
> ick <- SpatialPointsDataFrame(coordinates(tmp.lc), tmp.lc at data,
> proj4string=CRS(proj4string(tmp.lc)))
> bah <- ick[ ick at data$band1 %in% 22:24 ,]
>
>
>>  str(tmp.lc)
>
> Formal class 'SpatialGridDataFrame' [package "sp"] with 6 slots
>  ..@ data       :'data.frame': 22500 obs. of  1 variable:
>  .. ..$ band1: int [1:22500] 11 11 11 11 11 11 11 11 11 11 ...
>  ..@ grid       :Formal class 'GridTopology' [package "sp"] with 3 slots
>
>
>>  str(bah)
>
> Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
>  ..@ data       :'data.frame': 4032 obs. of  1 variable:
>  .. ..$ band1: int [1:4032] 22 22 22 22 22 22 22 22 22 23 ...
>
>
>>  sessionInfo()
>
> R version 2.10.1 (2009-12-14)
> i386-apple-darwin8.11.1
>
> locale:
> [1] C
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
>  [1] raster_0.9.9-4  spatstat_1.17-5 deldir_0.0-11   mgcv_1.6-1 rmacq_1.0-2
>     rgdal_0.6-12    maptools_0.7-29
>  [8] lattice_0.17-26 sp_0.9-56       foreign_0.8-39
>
> loaded via a namespace (and not attached):
> [1] Matrix_0.999375-33 grid_2.10.1        nlme_3.1-96 tcltk_2.10.1
> tools_2.10.1
>
>
> --
> --------------------------------------
> Don MacQueen
> Environmental Protection Department
> Lawrence Livermore National Laboratory
> Livermore, CA, USA
> 925-423-1062
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>