Addressing just the outermost cells in a raster?

Question

I have a DEM cropped to the extends of a shapefile with different glaciers (see here for data: https://drive.google.com/open?id=1ERFdsqDGLH1a_FbxwawE_gPm0Au0Q9vT).

DEM <- raster("***Glacier_Clip1.tif")
glaciers <- readOGR("***/01_Clipped_Data", "Glacier_Clip_Polygon")

glaciers_sf <- st_as_sf(glaciers)
glaciers_comb <- st_union(glaciers_sf)
glaciers_comb <- st_cast(glaciers_comb, "POLYGON")
glaciers_large <- glaciers_comb[st_area(glaciers_comb) > units::set_units(50000,"m^2")]
glaciers_large_sp <- as(glaciers_large, Class =  "Spatial")
glaciers_ras <- mask(DEM,glaciers_large_sp)

For further calculations I need to differentiate between all the "inner cells" and the "inner boundary cells" inside the glacier directly at the border and the "outer boundary cells" directly outside (those are not part of the cropped raster, so I somehow need to get them from the original DEM).

This image shows the three cell types (light blue, grey and yellow) I need to address somehow. Is it possible to get them via extract? Can anyone help me please?

Answer 1

I think you can do this with sneaky masking tricks, although it will be slow on very large rasters. Using your demo data,

library(raster)
library(sf)
setwd('C:/data/glaciers')

DEM <- raster("Glacier_Clip1.tif")
glaciers <- read_sf("Glacier_Clip_Polygon.shp")

glaciers_comb <- st_union(glaciers)
glaciers_comb <- st_cast(glaciers_comb, "POLYGON")
glaciers_large <- glaciers_comb[st_area(glaciers_comb) > units::set_units(50000,"m^2")]
glaciers_large_sp <- as(glaciers_large, Class =  "Spatial")
glaciers_ras <- mask(DEM,glaciers_large_sp)

This is the workhorse function. On a raster with patchy areas of data/nodata, it finds the cells on the outer bounds of the data patches and calls them 0, and everything else NA. You can set diag = FALSE for a stricter choice of which cells are selected.

outer_cells <- function(x = NULL, diag = TRUE) {
  focal(x, w = matrix(1, ncol = 3, nrow = 3),
        fun = function(win) {
          if(!is.na(win[5])) {
            return(NA_integer_)
          }

          if(diag == TRUE) {
            if(any(!is.na(win[c(-5)]))) { 
              0L
            } else {
              NA_integer_
            }  
          } else {
            if(any(!is.na(win[c(2,4,6,8)]))) {  
              0L
            } else {
              NA_integer_
            }
          }

        }
  )
}

Now, make a data presence/absence layer,

g_fp <- glaciers_ras
g_fp[!is.na(getValues(g_fp))] <- 1L

To get a raster that only contains the 'yellow cells', make a masking layer with outer_cells() and add it to your original DEM. data + NA = NA, so only the cells you need are in the output.

y_msk <- outer_cells(g_fp)
yellow <- y_msk + DEM

For 'grey' cells, invert your data mask, take the outer_cells() of that, and again add it to your DEM.

rcl <- data.frame('f' = c(1L, NA_integer_), 't' = c(NA_integer_, 1L))
g_fp_inv <- subs(g_fp, rcl)
g_msk <- outer_cells(g_fp_inv)
grey <- g_msk + DEM

Finally, for the light blue zone: flip g_msk, combine it with g_fp, and add it to your DEM to overwrite the 'grey' cells with NA

lb_msk <- subs(g_msk, rcl)
lb_msk <- lb_msk + g_fp
light_blue <- DEM + lb_msk

This is the three *_msks in QGIS, coloured up and zoomed in: