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})
4
  • 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. Nov 5, 2018 at 21:07
  • Thanks! I´ve edited my question.
    – user0405
    Nov 5, 2018 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). Nov 5, 2018 at 23:44
  • The raster package has many examples of how to do this. There are also 100s of examples of stackoverflow Nov 6, 2018 at 0:48

1 Answer 1

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 
4
  • Wow, thanks for this nice solution! Now I only have to figure out how to exclude all layers of the RasterBrick with more than 90% of pixels with the value 3 (90% cloud cover).
    – user0405
    Nov 6, 2018 at 17:21
  • I added some code to show how you can select layers Nov 6, 2018 at 18:34
  • csf shows me something different. It shows a count of 15.21429 ?
    – user0405
    Nov 6, 2018 at 19:25
  • I forgot a line! Sorry Nov 7, 2018 at 1:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.