1

I have a RasterBrick consisting of daily snow cover data with the values 1, 2 and 3 (1= snow, 2= no snow, 3= cloud-obscured).

Example of snow cover of one day:

enter image description here

> snowcover
class       : Large RasterBrick 
dimensions  : 26, 26, 2938  (nrow, ncol, nlayers)
resolution  : 231, 232  (x, y)
extent      : 718990, 724996, 5154964, 5160996  (xmin, xmax, ymin, ymax)
crs         : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84       
              +towgs84=0,0,0

Now I wish to interpolate the cloud-obscured pixels (but only where there are less than 90 % cloud cover in a single RasterLayer, otherwise the original values should be retained for this Layers).

For spatial interpolation I want to use a digital elevation model (same study area and already in same resolution) to extract upper and lower snowline boundaries for each Layer of the RasterBrick respectively. The upper snow line represents the elevation where all cloud-free pixels are classified as snow. The lower snowline identifies the altitude below which all cloud-free pixels are also snow-free.

enter image description here

> dem
class       : RasterLayer 
resolution  : 231, 232  (x, y)
extent      : 718990.2, 724996.2, 5154964, 5160996 (xmin, xmax, ymin, ymax)
crs         : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84       
              +towgs84=0,0,0 
values      : 1503, 2135  (min, max)

For the upper snowlines I need the minimum elevation of the snow-covered pixels (value = 1). Now all pixels of value 3 in a RasterLayer of the RasterBrick above this minimum elevation, should be reclassified as value 1 (assumed to be snow-covered).

For the lower snowline on the other hand I need to identify the maximum elevation of the no-snow-pixels (value = 2). Now all pixels of value 3 in a RasterLayer of the RasterBrick above this maximum elevation should be reclassified as value 2 (assumed to be snow-free).

Is this possible using R?

I tried to make use of the overlay function, but I got stuck there.

# For the upper snowline:
overlay <- overlay(snowcover, dem, fun=function(x,y){ x[y>=minValue(y[x == 1])] <- 1; x})
user0405
  • 15
  • 6
  • Yes, it is possible. To allow us to help you efficiently, please provide some example data generated by code (or from an example in the manual), and some code to show what you have tried. – Robert Hijmans Nov 05 '18 at 21:07
  • Thanks! I´ve edited my question. – user0405 Nov 05 '18 at 21:59
  • 1
    I may be nitpicking but Google Drive is not the best place to store the data that makes your question reproducible (link may break, unlikely to be archived by archive.org). The ideal would be to include the data in the text of the question (but in this case the XML is too big so you would have to use another format). Or else put the XML code on something more html based like pastebin (I just checked, the license there is the same as on SO). – Gabriel Devillers Nov 05 '18 at 23:44
  • The `raster` package has many examples of how to do this. There are also 100s of examples of stackoverflow – Robert Hijmans Nov 06 '18 at 00:48

1 Answers1

3

Here is some example data

library(raster)
dem <- raster(ncol=8, nrow=7, xmn=720145, xmx=721993, ymn=5158211, ymx=5159835, crs='+proj=utm +zone=32 +datum=WGS84')
values(dem) <- ncell(dem):1
snow <- setValues(dem, c(1, 1, rep(1:3, each=18)))
snow[,c(2,5)] <- NA
snow[3] <- 3


plot(snow)
lines(as(dem, 'SpatialPolygons'))
text(dem)

The plot shows the snow classes (1, 2, 3) with the elevation values on top. snow map

We can use mask, but need to deal with the missing values. 

msnow <- reclassify(snow, cbind(NA, 0))
# mask to get only the snow elevations
x <- mask(dem, msnow, maskvalue=1, inverse=TRUE)

# minimum elevation of the snow-covered cells
minsnow <- minValue(x)
minsnow 
#[1] 37

# snow elevation = 1
snowy <- reclassify(dem, rbind(c(-Inf, minsnow, NA), c(minsnow, Inf, 1)))
newsnow <- cover(snow, snowy)

s <- stack(dem, snow, newsnow)
names(s) <- c("elevation", "old_snow", "new_snow")

enter image description here

You were very close, as you can do

 r <- overlay(dem, snow, fun=function(e, s){ s[e >= minsnow] <- 1; s})

But note that that also overwrites high cells with no snow.

enter image description here

Which could be fixed like this:

r <- overlay(dem, snow, fun=function(e, s){ s[e >= minsnow & is.na(s)] <- 1; s})

To select layers with more than x% cells with value 3 (here I use a threshold of 34%):

threshold = .34
s <- stack(snow, snow+1, snow+2)
f <- freq(snow)
f 
#     value count
#[1,]     1    14
#[2,]     2    13
#[3,]     3    15
#[4,]    NA    14

nas <- f[is.na(f[,1]), 2]

ss <- subs(s, data.frame(from=3, to=1, subsWithNA=TRUE))
cs <- cellStats(ss, sum)
csf <- cs / (ncell(snow) - nas)
csf
#  layer.1   layer.2   layer.3 
#0.3571429 0.3095238 0.3333333 

i <- which(csf < threshold)
use <- s[[i]]
#use
class       : RasterStack 
dimensions  : 7, 8, 56, 2  (nrow, ncol, ncell, nlayers)
resolution  : 231, 232  (x, y)
extent      : 720145, 721993, 5158211, 5159835  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=32 +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
names       : layer.2, layer.3 
min values  :       2,       3 
max values  :       4,       5
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63