Calculating landscape metrics using R?

Question

I would like to calculate area-weighted mean shape index (AWMSI) and mean nearest neighbor distance (MMN) and other landscape metrics for a classified map (tif) using R.

I've found a github package called "LandscapeMetrics" https://rdrr.io/github/LandascapeMetrics/LandscapeMetrics/

My code is

install.packages("devtools")
library(devtools)
install_github("LandascapeMetrics/LandscapeMetrics")
library(LandscapeMetrics)
library(raster)
imagg <- raster("R/01.tif")
AWMSI(imagg)

I am getting

Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘is.projected’ for signature ‘"NULL"’
In addition: Warning message:
In rasterToPolygons(fragment_forest, fun = NULL, na.rm = TRUE, dissolve = TRUE) :
  no values in selection

How can I specify the class that I would like to calculate AWMSI for and other metrics?

Answer 1

Function AWMSI:

AWMSI <- function(x){
  class_forest=x==1
  fragment_forest=clump(class_forest)
  polygons=rasterToPolygons(fragment_forest, fun=NULL, na.rm=TRUE, dissolve=TRUE)
  areas_m2=gArea(polygons, byid = T)
  area_pixel=res(x)[1] * res(x)[2]
  areas_pixels=areas_m2/area_pixel
  perims_m=gLength(polygons, byid = T)
  perims_edges=perims_m/res(x)[1]
  d=data.frame(id = names(areas_m2), perims_edges, perims_m, areas_pixels, areas_m2)
  n=length(polygons)
  {
    LSI = function(a,TE) {
      b=trunc(sqrt(a))
      c=a - b^2
      minTE=rep(0,length(c))
      for (ii in 1:length(c)){
        if (c[ii]==0) minTE[ii]=4*b[ii]
        if (b[ii]^2<a[ii] & a[ii]<=b[ii]*(1+b[ii])) minTE[ii]=4 * b[ii] + 2
        if (a[ii] > b[ii]*(1+b[ii])) minTE[ii]=4 * b[ii] + 4
      }
      return(TE/minTE)
    }
  }
  d$SHAPE=LSI(d$areas_pixels, d$perims_edges)
  d$area_weighted=d$areas_m2/sum(areas_m2)
  d$AWMSI=(d$SHAPE)*(d$area_weighted)
  AWMSI=sum(d$AWMSI)
  return(AWMSI)
}

The function takes the raster and converts it to polygons. Trying to reproduce the error (example data):

library(LandscapeMetrics)
library(raster)

r <- raster(nrows=180, ncols=360, xmn=-10, xmx=10, ymn=-10, ymx=10, 
            crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0", 
            resolution = 1, vals = rep(x = 1:5, each = ncell(r)/5))

Is a geographic projection error?

AWMSI(r)
## [1] 1.333333
## Warning messages:
##   1: In RGEOSMiscFunc(spgeom, byid, "rgeos_area") :
##   Spatial object is not projected; GEOS expects planar coordinates
## 2: In RGEOSMiscFunc(spgeom, byid, "rgeos_length") :
##   Spatial object is not projected; GEOS expects planar coordinates

Warning message, no error... The same thing with meter projection:

# web mercator (in meters)
r_2 <- projectRaster(from = r, crs = '+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs')

AWMSI(r_2)
## [1] 1.5

If error came from rasterToPolygons(), is not projection, is raster's values. The same example with an empty raster:

r_3 <- raster(nrows=180, ncols=360, xmn=-10, xmx=10, ymn=-10, ymx=10, 
              crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0", 
              resolution = 1, vals = NULL)

AWMSI(r_3)
## Error in (function (classes, fdef, mtable)  : 
##             unable to find an inherited method for function ‘is.projected’ for signature ‘"NULL"’
##           In addition: Warning message:
##             In rasterToPolygons(fragment_forest, fun = NULL, na.rm = TRUE, dissolve = TRUE) :
##             no values in selection

Check your raster, there is an error with raster values

Answer 2

The warnings that you are receiving when adding the LandscapeMetrics package are not important, the package will function just fine. It is looking like your issue is that you either do not have a projection defined or the projection is geographic and not projected. It looks like there is a lot of hard-coding in the package that assumes projection units in meters. Also, not sure if you have realized that this package only returns "landscape level" (single value for entire landscape) metrics and nothing local?

If you want metrics assigned to a raster, sample locations or patches you will need to look elsewhere or create your own code. The spatialEco package has two functions for creating local landscape metrics. The land.metrics function calculates metrics for points (using radius) and polygons whereas, focal.lmetrics uses a specified sized moving window, much like Fragstats, and results in a raster of the specified metric.