Stratification
Fundamental to many structurally guided sampling approaches is the
use of stratification methods that allow for more effective and
representative sampling protocols. It is important to note that the data
sets being used as inputs are considered to be populations.
Currently, there are 5 functions associated with the
strat verb in the sgsR package:
strat_kmeans
strat_kmeans() uses kmeans clustering to produce an
sraster output.
#--- perform stratification using k-means ---#
strat_kmeans(mraster = mraster, # input
nStrata = 5) # algorithm will produce 4 strata
#> class : SpatRaster
#> dimensions : 277, 373, 1 (nrow, ncol, nlyr)
#> resolution : 20, 20 (x, y)
#> extent : 431100, 438560, 5337700, 5343240 (xmin, xmax, ymin, ymax)
#> coord. ref. : UTM Zone 17, Northern Hemisphere
#> source(s) : memory
#> name : strata
#> min value : 1
#> max value : 5
TIP!
plot = FALSE is the default for all functions.
plot = TRUE will visualize raster and vector ouputs.
strat_kmeans(mraster = mraster, # input
nStrata = 10, # algorithm will produce 10 strata
iter = 1000, # set minimum number of iterations to determine kmeans centers
algorithm = "MacQueen", # use MacQueen algorithm
plot = TRUE) # plot output
#> class : SpatRaster
#> dimensions : 277, 373, 1 (nrow, ncol, nlyr)
#> resolution : 20, 20 (x, y)
#> extent : 431100, 438560, 5337700, 5343240 (xmin, xmax, ymin, ymax)
#> coord. ref. : UTM Zone 17, Northern Hemisphere
#> source(s) : memory
#> name : strata
#> min value : 1
#> max value : 10
strat_quantiles
The strat_quantiles() algorithm divides data into
equally sized strata (nStrata). Similar to
strat_breaks() users can perform stratification on a single
mraster or also input a secondary mraster
(mraster2) and specify the desired number of strata
(nStrata2).
Note that the dual stratification output will result in a product of
\(nStrata * nStrata2\).
#--- perform quantiles stratification ---#
strat_quantiles(mraster = mraster$zq90,
nStrata = 6,
plot = TRUE)
#> class : SpatRaster
#> dimensions : 277, 373, 1 (nrow, ncol, nlyr)
#> resolution : 20, 20 (x, y)
#> extent : 431100, 438560, 5337700, 5343240 (xmin, xmax, ymin, ymax)
#> coord. ref. : UTM Zone 17, Northern Hemisphere
#> source(s) : memory
#> name : strata
#> min value : 1
#> max value : 6
#--- dual stratification - will produce 12 output strata ---#
strat_quantiles(mraster = mraster$zq90,
mraster2 = mraster$zsd,
nStrata = 3,
nStrata2 = 4)
#> class : SpatRaster
#> dimensions : 277, 373, 1 (nrow, ncol, nlyr)
#> resolution : 20, 20 (x, y)
#> extent : 431100, 438560, 5337700, 5343240 (xmin, xmax, ymin, ymax)
#> coord. ref. : UTM Zone 17, Northern Hemisphere
#> source(s) : memory
#> name : strata
#> min value : 1
#> max value : 12
strat_breaks
strat_breaks() stratifies data based on user-defined
breaks in mraster. A single mraster can be
used, or mraster2 can also be defined. breaks
and breaks2 coincide with the user defined breaks for
mraster and mraster2 respectively.
#--- perform stratification using user-defined breaks ---#
#--- define breaks for metric ---#
breaks <- c(seq(0,100,20))
breaks
#> [1] 0 20 40 60 80 100
#--- perform stratification using user-defined breaks ---#
values <- terra::values(mraster$zq90)
#--- define breaks for metric ---#
breaks2 <- c(5,10,15,20,25)
breaks2
#> [1] 5 10 15 20 25
Once the breaks are created, we can use them as input into the
strat_breaks function using the breaks and
breaks2 parameters.
#--- stratify on 1 metric only ---#
strat_breaks(mraster = mraster$pzabove2,
breaks = breaks,
plot = TRUE)
#> class : SpatRaster
#> dimensions : 277, 373, 1 (nrow, ncol, nlyr)
#> resolution : 20, 20 (x, y)
#> extent : 431100, 438560, 5337700, 5343240 (xmin, xmax, ymin, ymax)
#> coord. ref. : UTM Zone 17, Northern Hemisphere
#> source(s) : memory
#> name : strata
#> min value : 1
#> max value : 6
#--- stratify on 1 metric only ---#
strat_breaks(mraster = mraster$zq90,
breaks = breaks2,
plot = TRUE)
#> class : SpatRaster
#> dimensions : 277, 373, 1 (nrow, ncol, nlyr)
#> resolution : 20, 20 (x, y)
#> extent : 431100, 438560, 5337700, 5343240 (xmin, xmax, ymin, ymax)
#> coord. ref. : UTM Zone 17, Northern Hemisphere
#> source(s) : memory
#> name : strata
#> min value : 1
#> max value : 6
strat_poly
Forest inventories with polygon coverages summarizing forest
attributes such as species, management type, or photo-interpreted
estimates of volume can be stratified using
strat_poly().
TIP!
Users may wish to stratify based on categorical or empirical
variables that are not available through raster data (e.g. species from
forest inventory polygons).
Users define the input poly and its associated
attribute. A raster layer must be defined to
guide the spatial extent and resolution for the output stratification
polygon. Based on the vector or list of features,
stratification is applied and the polygon is rasterized into its
appropriate strata.
#--- load in polygon coverage ---#
poly <- system.file("extdata", "inventory_polygons.shp", package = "sgsR")
fri <- sf::st_read(poly)
#> Reading layer `inventory_polygons' from data source
#> `C:\Users\tgood\AppData\Local\Temp\RtmpCYyraS\Rinst5710142d19b\sgsR\extdata\inventory_polygons.shp'
#> using driver `ESRI Shapefile'
#> Simple feature collection with 632 features and 3 fields
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: 431100 ymin: 5337700 xmax: 438560 ymax: 5343240
#> Projected CRS: UTM_Zone_17_Northern_Hemisphere
#--- specify polygon attribute to stratify ---#
attribute <- "NUTRIENTS"
#--- specify features within attribute & how they should be grouped ---#
#--- as a single vector ---#
features <- c("poor", "rich", "medium")
In our example, attribute = "NUTRIENTS" and features
within, c("poor", "rich", "medium"), define the 3 desired
strata.
#--- stratify polygon coverage ---#
srasterpoly <- strat_poly(poly = fri, # input polygon
attribute = attribute, # attribute to stratify by
features = features, # features within attribute
raster = sraster, # raster to define extent and resolution for output
plot = TRUE) # plot output
#> Assigning a new crs. Use 'project' to transform a SpatRaster to a new crs
features can be grouped. In our example below,
rich and medium features are combined into a
single strata, while low is left in isolation. The 2
vectors are specified into a list, which will result in the output of 2
strata (low & rich/medium).
#--- or as multiple lists ---#
g1 <- "poor"
g2 <- c("rich", "medium")
features <- list(g1, g2)
strat_poly(poly = fri,
attribute = attribute,
features = features,
raster = sraster,
plot = TRUE,
details = TRUE)
#> Assigning a new crs. Use 'project' to transform a SpatRaster to a new crs
#> $raster
#> class : SpatRaster
#> dimensions : 277, 373, 1 (nrow, ncol, nlyr)
#> resolution : 20, 20 (x, y)
#> extent : 431100, 438560, 5337700, 5343240 (xmin, xmax, ymin, ymax)
#> coord. ref. : UTM Zone 17, Northern Hemisphere
#> source(s) : memory
#> name : strata
#> min value : 1
#> max value : 2
#>
#> $lookUp
#> strata features
#> 1 1 poor
#> 2 2 rich
#> 3 2 medium
#>
#> $poly
#> class : SpatVector
#> geometry : polygons
#> dimensions : 524, 2 (geometries, attributes)
#> extent : 431100, 438560, 5337700, 5343240 (xmin, xmax, ymin, ymax)
#> coord. ref. : UTM_Zone_17_Northern_Hemisphere
#> names : features strata
#> type : <chr> <int>
#> values : poor 1
#> poor 1
#> poor 1
details
details returns the output outRaster, the
stratification $lookUp table, and the polygon
($poly) used to drive the stratification based on
attributes and features specified by the users.
strat_map
Users may wish to pair stratifications. strat_map()
faciliates two stratifications of matching extent and resolution to be
mapped against one another to generate unique strata based on stratum
pairings.
This facilitates the user to generate stratifications detailing
quantitative and qualitative measures such as structure by species, or
multiple qualitative measures such as species by management type.
The total number of classes of the output sraster is
multiplicative of the number of input strata. For example, if the input
sraster has 3 strata and sraster2 has 4
strata, then the output of strat_map() will be 12
strata.
There may be occasions where stratum do not interact spatially, this
will result in fewer output strata.
#--- map srasters ---#
strat_map(sraster = srasterpoly, # strat_poly 3 class stratification
sraster2 = sraster, # strat_kmeans 4 class stratification
plot = TRUE)
#> class : SpatRaster
#> dimensions : 277, 373, 1 (nrow, ncol, nlyr)
#> resolution : 20, 20 (x, y)
#> extent : 431100, 438560, 5337700, 5343240 (xmin, xmax, ymin, ymax)
#> coord. ref. : UTM Zone 17, Northern Hemisphere
#> source(s) : memory
#> categories : label
#> name : strata
#> min value : 11
#> max value : 34
The convention for the numeric value of the output strata is the
concatenation (merging) of sraster strata and
sraster2 strata. Check $lookUP for a clear
depiction of this step.
strat_map(sraster = srasterpoly, # strat_poly 3 class stratification
sraster2 = sraster, # strat_poly 3 class stratification
stack = TRUE, # stack input and output strata into multi layer output raster
details = TRUE, # provide additional details
plot = TRUE) # plot output
#> Stacking sraster, sraster2, and their combination (stratamapped).
#> $raster
#> class : SpatRaster
#> dimensions : 277, 373, 3 (nrow, ncol, nlyr)
#> resolution : 20, 20 (x, y)
#> extent : 431100, 438560, 5337700, 5343240 (xmin, xmax, ymin, ymax)
#> coord. ref. : UTM Zone 17, Northern Hemisphere
#> source(s) : memory
#> names : strata, strata2, stratamapped
#> min values : 1, 1, 11
#> max values : 3, 4, 34
#>
#> $lookUp
#> strata strata2 stratamapped
#> 1 3 2 32
#> 2 3 1 31
#> 3 1 3 13
#> 4 1 4 14
#> 5 3 4 34
#> 6 3 3 33
#> 7 1 2 12
#> 8 1 1 11
#> 9 2 2 22
#> 10 2 3 23
#> 11 2 4 24
#> 12 2 1 21