6

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?

enter image description here

4

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: enter image description here

  • 1
    I will try that. Thanks a lot for the very detailed answer! – Florian Mlehliv Jul 4 at 12:00
  • 1
    So I tried your solution, @obrl_soil, but I have a follow-up question, that I edited in into the original one. Or should I post it as a new one? Not sure whats the best way here as my post is getting quite long... – Florian Mlehliv Jul 4 at 13:37
  • I've rolled your edit back to the original question. – Spacedman Jul 4 at 14:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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