How to speed up calculating proportion of landcover classes falling within larger pixels?

Wed, Aug 29, 2012 7:10 PM

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

How to speed up calculating proportion of landcover classes falling within larger pixels?

Thread (7 messages)