Skip to content

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

7 messages · Lyndon Estes, Robert J. Hijmans

Original
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
5 days later
3 days later
Lyndon Estes
# 
Hi Robert,

Just wanted to report back on the solution you gave me, which worked
perfectly as you wrote it.
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
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: