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

[Lyndon Estes]

Dear List,

I am currently trying to figure out the the proportions of two
different landcover classes (classified from Landsat imagery) that
fall within the extent of coarser MODIS pixels (I am trying to isolate
MODIS pixels associated with fairly pure cover types).  I want to do
this without changing the geometries of either raster, so the only
solution I can think of is to convert the MODIS pixels to polygons and
then use the extract function to do the work.  The code should look
something like this (adapting a solution I saw here
https://stat.ethz.ch/pipermail/r-sig-geo/2011-February/010999.html):

library(raster)
library(rgdal)

# Dummy of MODIS raster (but the actual extent I am working with),
each cell given a unique identifier
sincrs <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
+b=6371007.181 +units=m +no_defs"
r <- raster(xmn = 3230680, xmx = 3671058, ymn = -1712867, ymx =
-1343839, crs = sincrs)
res(r) <- 231.6564
r[] <- 1:ncell(r)

# Dummy 30m 7 class landcover raster (nonsense numbers, but is of the
extent I am working with)
lc <- raster(xmn = 3230756, xmx = 3671006, ymn = -1712837, ymx =
-1343777, crs = sincrs)
res(lc) <- 30
lc[] <- sample(0:7, ncell(lc), replace = TRUE)  # This is very slow

bs <- blockSize(r, minblocks = 200)  # Calculate size of rows to work on
rsp <- as(r,  "SpatialPixelsDataFrame")  # Convert MODIS dataset to
spatial pixels--allows ready conversion to polygons

# Output datasets
rcl3 <- r * 0
rcl4 <- r * 0
rtot <- r * 0

rcl3 <- writeStart(rcl3, "lc_class3.tif", overwrite = TRUE)
rcl4 <- writeStart(rcl4, "lc_class4.tif", overwrite = TRUE)
rtot <- writeStart(rtot, "lc_class4.tif", overwrite = TRUE)

# Start processing loop
for(i in 1:bs$n) {
  ids <- getValues(r, row = bs$row[i], nrows = bs$nrows[i]) # Get
pixel IDs from raster
  pix <- rsp[ids, ]  # Pull out spatialpixels from first block
  pol <- as(pix, "SpatialPolygonsDataFrame")  # Coerce them to polygons
  xv  <- extract(lc, pol[1:10, ])  #  Extract values
  writeValues(rcl3, sapply(xv, function(x) length(which(x == 3))),
bs$row[i])  # Write out sum class 3
  writeValues(rcl4, sapply(xv, function(x) length(which(x == 4))),
bs$row[i])  # Write out sum class 4
  writeValues(rtot, sapply(xv, function(x) length(x)), bs$row[i])  #
Write sum total pixels
}

rcl3 <- writeStop(rcl3)
rcl4 <- writeStop(rcl4)
rtot <- writeStop(rtot)

Note that I haven't actually run this all the way through to see if it
works properly (I stopped the run on the first iteration because it
hadn't finished after about 30-45 minutes), and my actual scenes have
plenty of NA values in them, but this is more or less the approach.

I would very much appreciate any pointers as to how I can do this more
efficiently.

Many thanks, Lyndon