[R-sig-Geo] How to speed up calculating proportion of landcover classes falling within larger pixels?

[Lyndon Estes]

Hi Robert,

Excellent, many thanks for adding that.  I will give that a try when
it is released.

For a name, something like "distill" or "fractionate" might provide a
reasonable description of what the function does (or at least in terms
of how I was using the approach)?

One thing to point out, in case it is of interest to others who will
use this approach, is that if the second case is used, as I did, the
resulting raster may have a quilt-like appearance if the edges of the
cells of the coarser raster do not line up with the edges of the outer
pixels of the finer resolution raster underneath it.  In my case, I
had results ranging from the low 40s to 65 for the total counts of 30
m cells whose centers intersected each of the 232 m MODIS pixels.

Thanks again,

Lyndon


On Wed, Aug 29, 2012 at 1:48 PM, Robert J. Hijmans <r.hijmans at gmail.com> wrote:
> 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 at 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
>
>