1

My input raster consists of multiple layers, each with an image area surrounded by no data values. These layers don't overlap completely, and I am trying to output a file which consists of only the intersection of all bands (the zone which has no NoData values on any layer).

The following works for a few layers, but not for the 50+ that I have in my real files (at least 3000x3000 pixels):

library(raster)

fin = "D:\\temp\\all_modes.pix"
fout = "D:\\temp\\test.pix"

inbands = stack(fin, bands = c(3:20))
NAvalue(inbands) = 0

# Not great:
#out = all(is.na(inbands) == FALSE) * inbands
#writeRaster(out, filename=fout, format="PCIDSK", dtype="INT2U", overwrite=TRUE, NAflag=0)

# A little better:
#mymask = all(as.logical(inbands))
#mask(inbands, mymask, filename=fout, format="PCIDSK", dtype="INT2U", overwrite=TRUE, NAflag=0)

# Even better, don't need to keep everything (but still not efficient):
#trim(all(as.logical(inbands)) * inbands, filename=fout, format="PCIDSK", dtype="INT2U", overwrite=TRUE, NAflag=0)

# Even better, calculations get smaller as we progress (is it possible to do even better?)
for(i in 1:nlayers(inbands)){
 band_i = subset(inbands, i)
 inbands = trim(as.logical(band_i) * inbands)
}
writeRaster(inbands, filename=fout, format="PCIDSK", dtype="INT2U", overwrite=TRUE, NAflag=0)

Any ideas on how to do this more efficiently / get it working with a large number of layers?

Benjamin
  • 11,560
  • 13
  • 70
  • 119

3 Answers3

2

I first proposed this:

x <- trim(inband, filename='fout')

But I now realize that would not get you what you want, as it would return the area (rows/columns) where at least one layer has a value that is not NA; rather than the area where all have a value that is not NA.

The below could be efficient. All cells with at least one NA value will become NA with sum (by default, na.rm=FALSE).

x <- sum(inband)
x <- trim(x)
r <- crop(inband, x)

perhaps followed by

r <- mask(r, x)

to set all cells in r to NA if they are NA in x

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • That is the union, not the intersection... and its equal to my original file. – Benjamin Apr 09 '11 at 04:00
  • A simple and concise solution, but I run out of memory with my large rasterStacks... :( – Benjamin Apr 15 '11 at 03:35
  • At what step do you run out of memory? If your RasterStack points to files on disk this should not happen. But if it happens anyway, you can use setOptions(chunksize= ) to circumvent this problem. – Robert Hijmans Apr 18 '11 at 00:18
1

Thanks for the answers, they gave me good ideas. I came up with this, which is much faster:

myraster = stack(fn, bands) # You get the idea
NAvalue(myraster) = 0

# Tranform to 1 where there is data    
logical_raster = as.logical(myraster)

# Make a raster with 1 in the zone of intersection
a = subset(myraster, 1)
values(a) = TRUE
for(i in 1:nlayers(myraster)) {
  a = a & logical_raster[[i]]
}

# Apply the "mask" and trim to intersection extent
myraster = myraster * a
intersect_only = trim(myraster)
Benjamin
  • 11,560
  • 13
  • 70
  • 119
0

I discovered that all/any are the ways to do logical AND/OR on a long list. Here are demonstrated two identical plots that accomplish your goal on small rasters.

#dummy data
m1 = cbind(matrix(NA, nrow = 4, ncol = 2), matrix(1, nrow = 4, ncol = 2))
m2 = t(m1)
m3 = matrix(rep(c(1, NA), 8), nrow = 4)
inbands = stack(lapply(list(m1, m2, m3), raster))

# first method, using & operator 
plot(inbands[[1]] & inbands[[2]] & inbands[[3]])

# second method, using `all`
plot(all(inbands))

On my system there are warnings for coercing numeric to logical. The following two methods avoid the warnings but might be slower? They are logically equivalent but you might time them against each other and against the second method above.

plot(all(!is.na(inbands)))
plot(!any(is.na(inbands)))
J. Win.
  • 6,662
  • 7
  • 34
  • 52
  • Thanks for the answer, will try on Monday. The other idea I had was to use extent() to figure out the minimum extent I need to keep, then crop() at the end, and then apply something similar to your solution. – Benjamin Apr 09 '11 at 13:47
  • @Benjamin. I've posted a better way. But, here is why do.call didnt work. http://stackoverflow.com/questions/2714159/how-to-wrap-a-function-that-only-takes-individual-elements-to-make-it-take-a-list – J. Win. Apr 12 '11 at 07:34
  • Thanks. There is also is.finite() which prevents the extra step. I think with raster I can simply use all(raster) which would simplify things. – Benjamin Apr 12 '11 at 11:53