Calculate structural connectivity — lsm_connectivity

Identifies patch and calculates area in hectare.

Usage

lsm_connectivity_structural(
  input,
  output = NULL,
  zero_as_null = FALSE,
  id_direction = 8,
  area_round_digit = 0,
  area_unit = "ha",
  map_connec_struct = TRUE,
  map_connec_struct_area = FALSE,
  nprocs = 1,
  memory = 300
)

Arguments

input: [character] Habitat map (binary classification: e.g., 1/0 or 1/NA) in GRASS.
output: [character] Output map base name in GRASS.
zero_as_null: [logical] If TRUE, non-habitat (0) cells are converted to NULL.
id_direction: [numeric] Neighborhood for clumping (4 or 8).
area_round_digit: [integer] Decimal digits for area rounding.
area_unit: [character] Area unit: "ha", "m2", or "km2".
map_id: [logical] Keep fragment ID raster?
map_ncell: [logical] Output raster with fragment cell counts?
map_area: [logical] Output raster with fragment area?

Examples

library(lsmetrics)
library(terra)

# read habitat data
r <- lsmetrics::lsm_toy_landscape(proj_type = "meters")

# plot
plot(r, legend = FALSE, axes = FALSE, main = "Binary habitat")
plot(as.polygons(r, dissolve = FALSE), lwd = .1, add = TRUE)
plot(as.polygons(r), add = TRUE)
text(r)


# find grass
path_grass <- system("grass --config path", inter = TRUE) # windows users need to find the grass gis path installation, e.g. "C:/Program Files/GRASS GIS 8.3"

# create grassdb
rgrass::initGRASS(gisBase = path_grass,
                  SG = r,
                  gisDbase = "grassdb",
                  location = "newLocation",
                  mapset = "PERMANENT",
                  override = TRUE)
#> gisdbase    grassdb 
#> location    newLocation 
#> mapset      PERMANENT 
#> rows        16 
#> columns     16 
#> north       7525600 
#> south       7524000 
#> west        234000 
#> east        235600 
#> nsres       100 
#> ewres       100 
#> projection:
#>  PROJCRS["WGS 84 / UTM zone 23S",
#>     BASEGEOGCRS["WGS 84",
#>         ENSEMBLE["World Geodetic System 1984 ensemble",
#>             MEMBER["World Geodetic System 1984 (Transit)"],
#>             MEMBER["World Geodetic System 1984 (G730)"],
#>             MEMBER["World Geodetic System 1984 (G873)"],
#>             MEMBER["World Geodetic System 1984 (G1150)"],
#>             MEMBER["World Geodetic System 1984 (G1674)"],
#>             MEMBER["World Geodetic System 1984 (G1762)"],
#>             MEMBER["World Geodetic System 1984 (G2139)"],
#>             ELLIPSOID["WGS 84",6378137,298.257223563,
#>                 LENGTHUNIT["metre",1]],
#>             ENSEMBLEACCURACY[2.0]],
#>         PRIMEM["Greenwich",0,
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>         ID["EPSG",4326]],
#>     CONVERSION["UTM zone 23S",
#>         METHOD["Transverse Mercator",
#>             ID["EPSG",9807]],
#>         PARAMETER["Latitude of natural origin",0,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8801]],
#>         PARAMETER["Longitude of natural origin",-45,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8802]],
#>         PARAMETER["Scale factor at natural origin",0.9996,
#>             SCALEUNIT["unity",1],
#>             ID["EPSG",8805]],
#>         PARAMETER["False easting",500000,
#>             LENGTHUNIT["metre",1],
#>             ID["EPSG",8806]],
#>         PARAMETER["False northing",10000000,
#>             LENGTHUNIT["metre",1],
#>             ID["EPSG",8807]]],
#>     CS[Cartesian,2],
#>         AXIS["(E)",east,
#>             ORDER[1],
#>             LENGTHUNIT["metre",1]],
#>         AXIS["(N)",north,
#>             ORDER[2],
#>             LENGTHUNIT["metre",1]],
#>     USAGE[
#>         SCOPE["Navigation and medium accuracy spatial referencing."],
#>         AREA["Between 48°W and 42°W, southern hemisphere between 80°S and equator, onshore and offshore. Brazil."],
#>         BBOX[-80,-48,0,-42]],
#>     ID["EPSG",32723]] 


# import raster from r to grass
rgrass::write_RAST(x = r, flags = c("o", "overwrite", "quiet"), vname = "r", verbose = FALSE)

# structural connectivity
lsmetrics::lsm_connectivity_structural(input = "r",
                                       zero_as_null = FALSE,
                                       map_connec_struct = TRUE,
                                       map_connec_struct_area = TRUE)
#> Converting zero as null
#> Identifying fragments
#> Identifying fragments
#> Calculation area
#> Assigning area
#>    0%   6%  12%  18%  25%  31%  37%  43%  50%  56%  62%  68%  75%  81%  87%  93% 100%
#> Cleaning data
#> Converting null to zero
#> Identifying the fragments
#> Converting zero as null
#> Identifying fragments
#> Creating matrix
#> Filling holes
#> Cleaning data
#> WARNING: Raster map
#>          <r_connec_struct_null_patch_binary_aux_fill_hole_matrix_id@PERMANENT>
#>          is a base map for
#>          <r_connec_struct_null_patch_binary_aux_fill_hole_matrix_id_fill@PERMANENT>.
#>          Remove forced.
#> Calculation area
#> Assigning area
#>    0%   6%  12%  18%  25%  31%  37%  43%  50%  56%  62%  68%  75%  81%  87%  93% 100%
#> Cleaning data
#>    0%   6%  12%  18%  25%  31%  37%  43%  50%  56%  62%  68%  75%  81%  87%  93% 100%
#> Cleaning data

# files
rgrass::execGRASS(cmd = "g.list", type = "raster")
#> r
#> r_struct_connec
#> r_struct_connec_area

# import from grass to r
r_struct_connec_area <- rgrass::read_RAST("r_struct_connec_area", flags = "quiet", return_format = "terra")

plot(r_struct_connec_area, legend = FALSE, axes = FALSE, main = "Structural connected area (ha)")
plot(as.polygons(r, dissolve = FALSE), lwd = .1, add = TRUE)
plot(as.polygons(r), add = TRUE)
text(r_struct_connec_area)


# import from grass to r
r_struct_connec <- rgrass::read_RAST("r_struct_connec", flags = "quiet", return_format = "terra")

plot(r_struct_connec, legend = FALSE, axes = FALSE, main = "Structural connectivity (ha)")
plot(as.polygons(r, dissolve = FALSE), lwd = .1, add = TRUE)
plot(as.polygons(r), add = TRUE)
text(r_struct_connec)


# delete grassdb
unlink("grassdb", recursive = TRUE)