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.es...@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]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
