Skip to contents

Create buffers and calculate statistics

Usage

lsm_buffer_statistic(
  input,
  output = NULL,
  landscape_metric,
  landscape_metric_has_null = FALSE,
  point,
  distance,
  column_prefix,
  method = "average",
  percentile = NULL
)

Arguments

input

[character=""]

output

[character=""]
Map name output inside GRASS Data Base.

landscape_metric

[character=""]

landscape_metric_has_null

[character=""]

distance

[character=""]

column_prefix

[character=""]

method

[character=""]
Univariate statistics: number, null_cells, minimum ,maximum, range, average, stddev, variance, coeff_var, sum, first_quartile ,median, third_quartile, percentile

percentile

[character=""]

points

[character=""]

Examples

library(lsmetrics)
library(terra)
#> terra 1.7.55

# read habitat data
f <- system.file("raster/toy_landscape_habitat.tif", package = "lsmetrics")
r <- terra::rast(f)

# 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["Engineering survey, topographic mapping."],
#>         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)

# points
p <- terra::vect(data.frame(x = c(235150, 234450, 235150),
                            y = c(7525250, 7524850, 7524350)),
                 geom = c("x", "y"),
                 crs = "+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs +type=crs")

# plot
plot(r, legend = FALSE, axes = FALSE, main = "Points")
plot(as.polygons(r, dissolve = FALSE), lwd = .1, add = TRUE)
plot(as.polygons(r), add = TRUE)
plot(p, cex = 3, add = TRUE)


# import vector from r to grass
rgrass::write_VECT(x = p, flags = c("o", "overwrite", "quiet"), vname = "p")

# area
lsmetrics::lsm_fragment_area(input = "r")
#> Converting zero as null
#> Identifying the fragmentes
#> Counting the cell number of fragmentes
#> First pass
#>    0%   6%  12%  18%  25%  31%  37%  43%  50%  56%  62%  68%  75%  81%  87%  93% 100%
#> Writing output map
#>    0%   6%  12%  18%  25%  31%  37%  43%  50%  56%  62%  68%  75%  81%  87%  93% 100%
#> Calculating the area of fragmentes

# buffer statistics
lsmetrics::lsm_buffer_statistic(input = "r",
                                landscape_metric = "r_fragment_area_ha",
                                landscape_metric_has_null = TRUE,
                                point = "p",
                                distance = 200,
                                column_prefix = "area",
                                method = "average")
#> Buffering features...
#>   33%  66% 100%
#> Cleaning buffers...
#> Building parts of topology...
#> Building topology for vector map <r_p_buffer200@PERMANENT>...
#> Registering primitives...
#> 
#> Snapping boundaries...
#> Reading features...
#> Snap vertices Pass 1: select points
#>    0%  33%  66% 100%
#> Snap vertices Pass 2: assign anchor vertices
#>    4%  10%  15%  20%  26%  31%  36%  42%  47%  52%  57%  62%  68%  73%  78%  84%  89%  94% 100%
#> Snap vertices Pass 3: snap to assigned points
#>    0%  33%  66% 100%
#> Breaking polygons...
#> Breaking polygons (pass 1: select break points)...
#>   33%  66% 100%
#> Breaking polygons (pass 2: break at selected points)...
#>   33%  66% 100%
#> Removing duplicates...
#>   33%  66% 100%
#> Breaking boundaries...
#>    0%  33%  66% 100%
#> Removing duplicates...
#>   33%  66% 100%
#> Cleaning boundaries at nodes
#>   33%  66% 100%
#> Building topology for vector map <r_p_buffer200@PERMANENT>...
#> Building areas...
#>    0%  33%  66% 100%
#> Removing dangles...
#>   33%  66% 100%
#> Removing bridges...
#>   33%  66% 100%
#> Attaching islands...
#> Building topology for vector map <r_p_buffer200@PERMANENT>...
#> Attaching islands...
#>    0%  33%  66% 100%
#> Calculating centroids for areas...
#>   33%  66% 100%
#> Building topology for vector map <r_p_buffer200@PERMANENT>...
#> Registering primitives...
#> 
#> Building areas...
#>    0%  33%  66% 100%
#> Attaching islands...
#>    0%  33%  66% 100%
#> Attaching centroids...
#>    0%  33%  66% 100%
#> Processing input data (3 categories)...
#> Adding columns '['area_average DOUBLE PRECISION']'
#> Adding column <area_average> to the table
#>    6%  12%  18%  25%  31%  37%  43%  50%  56%  62%  68%  75%  81%  87%  93% 100%
#>  100%
#> Updating the database ...
#> Statistics calculated from raster map <r_fragment_area_ha_zero> and
#> uploaded to attribute table of vector map <r_p_buffer200>.

# files
# rgrass::execGRASS(cmd = "g.list", type = "vector")

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

# import buffer
v_buffer <- rgrass::read_VECT("r_p_buffer200", flags = "quiet")

# plot
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)
plot(v_buffer, "area_average", legend = FALSE, alpha = .5, add = TRUE)
plot(p, col = "gray30", cex = 3, add = TRUE)
text(v_buffer, col = "white", labels = "area_average", cex = .7)


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