Hi Lyndon,

I have added this approach to 'raster' under the name of 'layerize' (is
there a better term?). There are two variations.

1) In the simple case the function transforms a RasterLayer with values
representing classes (e.g. land cover type) into a RasterBrick with a
presence/absence layer for each class. This is a simple wrapper around a
call to 'calc', but I might change that to increase speed.

2) In the other (your) case you need to also provide a second RasterLayer
that overlaps with the first, but has larger cells. The layers in the
resulting RasterBrick will have counts for each class.

Note that if the larger cell RasterLayer can be obtained by aggregation of
the smaller cell RasterLayer it should be more efficient to use the simple
case, followed by aggregate(x, , sum)   (( in your case, you could also
choose to aggregate the output of the simple case to the approximate
resolution/origin desired and then use resample to get the values where you
want them; but that is a more approximate solution ))

Robert

On Tue, Aug 28, 2012 at 1:06 PM, Lyndon Estes <lyndon.estes@gmail.com>wrote:

> Hi Robert,
>
> Just wanted to report back on the solution you gave me, which worked
> perfectly as you wrote it.
>
> >> library(raster)
> >>
> >> library(rgdal)
> >>
> >> r <- raster(ncol=10, nrow=10)  # EP MODIS pixels
> >> eplc <- raster()
> >> set.seed(0)
> >> eplc[] <- round(  runif(ncell(eplc)) + 3)
> >> eplc[10000:11000] <- NA
> >>
> >> bs <- blockSize(r, minblocks = 200)  # Calculate size of rows to work on
> >>
> >> # Output datasets
> >> rcl3 <- writeStart(r, "lc_class3.tif", overwrite = TRUE)
> >> rcl4 <- writeStart(r, "lc_class4.tif", overwrite = TRUE)
> >> rtot <- writeStart(r, "lc_aall.tif", overwrite = TRUE)
> >>
> >> # Start processing loop
> >> for(i in 1:bs$n) {
> >>    print(paste("Iteration", i))
> >>    # extent of this (modis) block
> >>    e <- extent(xmin(r), xmax(r), yFromRow(r, bs$row[i]+bs$nrows[i]-1) -
> 0.5 * yres(r), yFromRow(r, bs$row[i])+0.5 * yres(r))
> >>
> >>    vo1 <- vo2 <- vo3 <- rep(NA, bs$nrows[i] * ncol(r))
> >>
> >>    int <- intersect(e, extent(eplc))
> >>    if (!is.null(int)) {
> >>    b <- crop(eplc, int)
> >>    xy <- xyFromCell(b, 1:ncell(b))
> >>    mc <- cellFromXY(r, xy)
> >>
> >>            # here I use table, there could be other ways (e.g. tapply)
> to achieve the below
> >> v <- table(mc, getValues(b))
> >>    cells <- as.integer(rownames(v))
> >> modcells <- cellFromRowCol(r, bs$row[i], 1) : cellFromRowCol(r,
> bs$row[i]+bs$nrows[i]-1, ncol(r))
> >> m <- match(cells, modcells)
> >>
> >> cn <- as.integer(colnames(v))
> >> if (3 %in% cn) {
> >> vo1[m] <- v[, '3']
> >> }
> >> if (4 %in% cn) {
> >> vo2[m] <- v[, '4']
> >> }
> >> vo3[m] <- apply(v, 1, sum)
> >> rcl3 <- writeValues(rcl3, vo1, bs$row[i])  # Write out sum class 3
> >> rcl4 <- writeValues(rcl4, vo2, bs$row[i])  # Write out sum class 4
> >> rtot <- writeValues(rtot, vo3, bs$row[i])  # Write sum total pixels
> >>   }
> >> }
> >>
> >> rcl3 <- writeStop(rcl3)
> >> rcl4 <- writeStop(rcl4)
> >> rtot <- writeStop(rtot)
> >>
>
> I tweaked it to run with my actual grids, and it took only 6 minutes
> to run, producing 4 outputs each with the following characteristics:
>
> class       : RasterLayer
> dimensions  : 1593, 1901, 3028293  (nrow, ncol, ncell)
> resolution  : 231.6564, 231.6564  (x, y)
> extent      : 3230680, 3671058, -1712867, -1343839  (xmin, xmax, ymin,
> ymax)
> coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
> +b=6371007.181 +units=m +no_defs
> values      :
> layer name  : eplc_ag
> min value   : 0
> max value   : 64
>
> Needless to say, that is an enormous improvement on the 3 day run I
> had using the extract approach.
>
> Thanks again for the advice and for the raster package.
>
> Cheers, Lyndon
>

	[[alternative HTML version deleted]]

