GISFrag metric

gisfrag.R

## GISFrag metric
## https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19950017676.pdf
## 1) Produce a proximity (distance) map between 'patches'
## 2) GISFrag == mean of all values on the proximity map
## 3) Large GISFrag values reflect low forest fragmentation, low values == high fragmentation

gisFrag <- function(x, ...) {
  ## x needs to be a raster where cells with suitable habitat are coded as 1
  ## unsuitable cells coded with 0
  ## extract cell numbers for suitable cells
  require(raster)
  cell_ext <- Which(x == 1, cells = TRUE)
  if (length(cell_ext) == 0) {
    ## if no cells are suitable return NA
    return(NA)
  } else if (length(cell_ext)) {
    x[x == 1] <- NA
    if (all(is.na(values(x)))) {
      return(0)
    } else {
      d <- raster::distance(x, doEdge = TRUE)
      ## convert to km and calculate mean distance
      d <- mean(extract(d, cell_ext), na.rm = TRUE) / 1000
    }
  } else {
    ## if all cells are suitable return 0
    d <- 0
  }
  ## invert so larger values == more fragmentation
  ## in the original pub. above, larger values == less fragmentation
  return(1/d)
}

## TEST ##
test <- FALSE
if (test) {
  r <- brick(nl = 5, ncol = 36, nrow = 18)
  ## high frag == low connectivity
  set.seed(64); r[[1]][] <- 0; r[[1]][c(1, 648)] <- 1
  ## mid frag
  r[[2]][] <- rbinom(n = ncell(r), size = 1, prob = 0.5)
  ## low frag == high connectivity
  r[[3]][] <- rbinom(n = ncell(r), size = 1, prob = 0.90)
  ## no fragmentation
  r[[4]][] <- 1
  ## no suitable cells
  r[[5]][] <- 0
  ## green = suitable, red = unsuitable
  spplot(r, col.regions = c("red", "darkgreen"), colorkey = FALSE)
  sapply(X = 1:nlayers(r), FUN = function(i) gisFrag(r[[i]]))
}

Thanks for sharing this. Any suggestions for speeding up on larger rasters?

What size raster are you working with?

Actually, found the issue on my end - all raster values were equal to 1 so the resulting raster was all NAs. I modified your function to deal with that scenario:

gisFrag <- function(x, ...) {
cell_ext <- Which(x == 1, cells = TRUE)
if (length(cell_ext)){
x[x == 1] <- NA
if(all(is.na(values(x)))){
return(0)
} else{
d <- raster::distance(x, doEdge = TRUE, )
d <- mean(extract(d, cell_ext), na.rm = TRUE)/1000
}
} else{
d <- 0
}
return(1/d)
}