Hi,

To follow up on my last with an update.  I modified the code a bit to run
with my data and set it on its way on our linux server.  It seemed to be
running, and after the 3 days or so it took to go through all iterations,
it failed on me.

Here's the contents of the .Rout file:

R version 2.15.1 (2012-06-22) -- "Roasted Marshmallows"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-unknown-linux-gnu (64-bit)

<snip>

During startup - Warning message:
Setting LC_CTYPE failed, using "C"

> library(raster)
Loading required package: sp
raster 2.0-08 (27-June-2012)
> library(rgdal)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.9.1, released 2012/05/15
Path to GDAL shared files: /usr/local/share/gdal
Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
Path to PROJ.4 shared files: (autodetected)
>
> # Paths
> p.wk <- "/home/lestes/Zambia/LC/"
>
> # Read in base rasters
> setwd(p.wk)
> r <- raster("epmodf.tif")  # EP MODIS pixels
> # spm <- raster("spmodf.tif")  # SP MODIS pixels
> eplc <- raster("zam_eplc123_sinus.tif")  # EP LC
> # splc <- raster("zam_splc_sinus.tif")  # SP LC
>
> bs <- blockSize(r, minblocks = 200)  # Calculate size of rows to work on
> rsp <- as(r,  "SpatialPixelsDataFrame")  # Convert MODIS dataset to
spatial pixels
>
> # 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_aall.tif", overwrite = TRUE)
>
> # Start processing loop
> for(i in 1:bs$n) {
+   print(paste("Iteration", i))
+   v <- getValues(r, row = bs$row[i], nrows = bs$nrows[i])
+   vo1 <- v
+   vo2 <- v
+   vo3 <- v
+   ids <- v[!is.na(v)]
+   pix <- rsp[ids, ]  # Pull out pixels from first block
+   pol <- as(pix, "SpatialPolygonsDataFrame")  # Coerce them to polygons
+   xv  <- extract(eplc, pol)  #  Extract values
+   vo1[!is.na(v)] <- sapply(xv, function(x) length(which(x == 3)))
+   vo2[!is.na(v)] <- sapply(xv, function(x) length(which(x == 4)))
+   vo3[!is.na(v)] <- sapply(xv, function(x) length(x))
+   writeValues(rcl3, vo1, bs$row[i])  # Write out sum class 3
+   writeValues(rcl4, vo2, bs$row[i])  # Write out sum class 4
+   writeValues(rtot, vo3, bs$row[i])  # Write sum total pixels
+ }
[1] "Iteration 1"
<snip>
[1] "Iteration 200"

Error in validObject(.Object) :
  invalid class "SpatialPixels" object: bbox should never contain infinite
values
Calls: [ ... SpatialPixels -> new -> initialize -> initialize -> validObject
In addition: Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf
3: In min(x) : no non-missing arguments to min; returning Inf
4: In max(x) : no non-missing arguments to max; returning -Inf

I haven't fully verified because I am not sure where exactly it happened in
the run, but I believe the failure was caused by this:

ids <- v[!is.na(v)]

With a case where length(ids) was 0. So a conditional statement in the code
should fix that:

 if(length(ids) > 0) {
    pix <- rsp[ids, ]  # Pull out pixels from first block
    pol <- as(pix, "SpatialPolygonsDataFrame")  # Coerce them to polygons
    xv  <- extract(eplc, pol)  #  Extract values
    vo1[!is.na(v)] <- sapply(xv, function(x) length(which(x == 3)))
    vo2[!is.na(v)] <- sapply(xv, function(x) length(which(x == 4)))
    vo3[!is.na(v)] <- sapply(xv, function(x) length(x))
    writeValues(rcl3, vo1, bs$row[i])  # Write out sum class 3
    writeValues(rcl4, vo2, bs$row[i])  # Write out sum class 4
    writeValues(rtot, vo3, bs$row[i])  # Write sum total pixels
  }

But again, not yet verified.

However, the main thing that I can't figure is why my results were never
written into the raster files as each block completed?  As written, I
expected the output *.tif files to appears in my output directory, but they
never showed up, Instead, I have a whole of temporary rasters in /tmp
(which later caused a failure because /tmp hit its limit).

So, if anyone could point me in the right direction for getting things
written properly to disk, I would appreciate it. I would also very much
appreciate any hints for how to write a speedier code for addressing this
particular application.

Thanks, Lyndon


error is due to the fact that I ended up with a line for


On Sat, Aug 18, 2012 at 6:38 PM, Lyndon Estes <lyndon.estes@gmail.com>wrote:

> 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
>

	[[alternative HTML version deleted]]

