Identifies fragments and calculates area in a selected unit.
Usage
lsm_area_fragment(
input,
output = NULL,
zero_as_null = FALSE,
id_direction = 8,
area_round_digit = 0,
area_unit = "ha",
map_fragment_id = FALSE,
map_fragment_ncell = FALSE,
map_fragment_area = TRUE,
table_fragment_area = FALSE
)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_fragment_id
[logical]Keep fragment ID raster?- map_fragment_ncell
[logical]Output raster with fragment cell counts?- map_fragment_area
[logical]Output raster with fragment area?- table_fragment_area
[logical]Output CSV with fragment area and summary?
Examples
library(lsmetrics)
library(terra)
#> terra 1.8.60
# 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)
# area
lsmetrics::lsm_area_fragment(input = "r", map_fragment_id = TRUE, table_fragment_area = TRUE)
#> Converting zeros to null
#> Identifying fragments
#> Calculation area
#> Assigning area
#> 0% 6% 12% 18% 25% 31% 37% 43% 50% 56% 62% 68% 75% 81% 87% 93% 100%
#> Exporting table
#> Cleaning data
# files
rgrass::execGRASS(cmd = "g.list", type = "raster")
#> r
#> r_fragment_area
#> r_fragment_id
# import from grass to r
r_fragment_id <- rgrass::read_RAST("r_fragment_id", flags = "quiet", return_format = "terra")
# plot
plot(r_fragment_id, legend = FALSE, axes = FALSE, main = "Fragment id")
plot(as.polygons(r, dissolve = FALSE), lwd = .1, add = TRUE)
plot(as.polygons(r), add = TRUE)
text(r_fragment_id)
# import from grass to r
r_fragment_area <- rgrass::read_RAST("r_fragment_area", flags = "quiet", return_format = "terra")
plot(r_fragment_area, legend = FALSE, axes = FALSE, main = "Fragment area (ha)")
plot(as.polygons(r, dissolve = FALSE), lwd = .1, add = TRUE)
plot(as.polygons(r), add = TRUE)
text(r_fragment_area)
# table
da_fragment_area <- readr::read_csv("r_fragment_area.csv", show_col_types = FALSE)
da_fragment_area
#> # A tibble: 6 × 3
#> fid area ncell
#> <dbl> <dbl> <dbl>
#> 1 1 9 9
#> 2 2 25 25
#> 3 3 39 39
#> 4 4 2 2
#> 5 5 34 34
#> 6 6 1 1
da_fragment_area_summary <- readr::read_csv("r_fragment_area_summary.csv",
show_col_types = FALSE)
da_fragment_area_summary
#> # A tibble: 1 × 6
#> n_frag area_mean area_sd area_cv area_max area_min
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 6 18.3 16.6 90.3 39 1
# delete grassdb
unlink("grassdb", recursive = TRUE)
unlink("r_fragment_area.csv")
unlink("r_fragment_area_summary.csv")