2

I have a raster cropped to the extent of some glaciers and divided it into three masks to address inner glacial cells, border cells and non-glacier cells adjacent to the glacier. Data can be found here: https://drive.google.com/open?id=1ERFdsqDGLH1a_FbxwawE_gPm0Au0Q9vT

EDIT: Data preparation was as follows

glaciers <- readOGR("***", "Glacier_Clip_Polygon")
DEM <- raster("***Glacier_Clip1.tif")
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)

To divide the masked DEM into my three classes I used the function proposed by @obrl_soil:

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_
            }
          }

        }
  )
}

## presence/absence layer with no DEM-information
g_fp <- glaciers_ras
g_fp[!is.na(getValues(g_fp))] <- 1L

## glacier-adjacent cells 
y_msk <- outer_cells(g_fp)
glacier_adj_cells <- y_msk + DEM
plot(glacier_adj_cells)

## glacier marginal cells 
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)
glacier_marg_cells <- g_msk + DEM
plot(glacier_marg_cells)

## inner glacier cells
lb_msk <- subs(g_msk, rcl)
lb_msk <- lb_msk + g_fp
glacier_inn_cells <- DEM + lb_msk
plot(glacier_inn_cells)

Now I have three raster layers with my three classes for various glaciers. From the inner cells mask of each glacier I now want to draw random cells (in order to lay a buffer and calculate the slope and do some more stuff). The ratio between the random cells and the number of glacier marginal cells should be approx. 0.3.

I want to draw the same ratio from each of my glaciers, so I would have to address them separately. Is there a possibility of doing this in a raster object, maybe by using the polygons as markers of some kind?

I hope, this edit makes my problem a bit clearer.

  • I suspect you're not getting a response here because its a bit unclear. Could you describe what your data files are (with their names) and maybe read them in and plot them so we can see what you want? – Spacedman Jul 6 '19 at 7:33
  • Thanks for the advice! I hope, my question is clearer now? – Florian Mlehliv Jul 6 '19 at 9:27
1

Here is a suggestion for a framework. I assume it could be done more effectively and I might have misunderstood some aspects of your question. A couple of things to note is the case where the 0.3 ratio of marginal cells exceeds the number of inner cells, where I have simply defined all inner raster cells as the sample. There are also a couple of polygons that hold glacier rasters, but do not hold any inner raster cells (as can be seen in the figure below).

I originally provided the answer drawing random raster cells, but seeing as you want to create buffers around the samples for further analysis it seemed better to create random points.

library(raster)
library(rgdal)

# Reads in the data, provided it is located in the working directory
glaciersPoly <- readOGR(".", "Glacier_Clip_Polygon")
adjacent <- raster('adjacent_cells.tif')
marginal <- raster('border_cells.tif')
inner <- raster('inner_cells.tif')

# Combine glacier rasters into one
glaciersRast <- merge(marginal, inner)

# Plot for visualisation
plot(glaciersPoly)
plot(glaciersRast, add = TRUE, alpha = 0.5)


# Loop over each polygon to draw samples
for(i in 1:length(glaciersPoly)){
  # Ignore polygons not overlapping any raster cells
  if(all(is.na(extract(glaciersRast, glaciersPoly[i,])[[1]]))){
    next
  } else {
    # Find the number of cells of the marginal raster that fall within the polygon
    margin <- extract(marginal, glaciersPoly[i,])
    # Find the number of cells that the desired 0.3 ratio equates to (rounded)
    marginRatio <- round(length(na.omit(margin[[1]])) * 0.3)
    # Clip the inner raster by the polygon
    innerClipped <- mask(inner, glaciersPoly[i,])
    # If the number of cells in innerClipped is less than what is required to meet 
    # the ratio I simply create a point for each raster cell and mark it red, as I'm
    # not sure how you would like to handle this case
    if(ncell(na.omit(innerClipped[])) < marginRatio){
      samplePts <- rasterToPoints(innerClipped)
      points(samplePts, pch = 3, cex = 0.4, col = "red") 
    } else {
      # Draw the sample of random points from the inner raster cells
      samplePts <- sampleRandom(innerClipped, size = marginRatio, sp = TRUE)
      points(samplePts, pch = 3, cex = 0.4)
    }
  }
} 

Here is the output, where the sample points are in black and red. enter image description here

  • 1
    Thanks, that brought me on the right track – Florian Mlehliv Jul 8 '19 at 14:20

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.