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

[Lyndon Estes]

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