Currently, there are 9 functions associated with the
sample
verb in the sgsR
package:
Algorithm | Description | Reference |
---|---|---|
sample_srs() |
Simple random | |
sample_systematic() |
Systematic | |
sample_strat() |
Stratified | Queinnec, White, & Coops (2021) |
sample_sys_strat() |
Systematic Stratified | |
sample_nc() |
Nearest centroid | Melville & Stone (2016) |
sample_clhs() |
Conditioned Latin hypercube | Minasny & McBratney (2006) |
sample_balanced() |
Balanced sampling | Grafström, A. Lisic, J (2018) |
sample_ahels() |
Adapted hypercube evaluation of a legacy sample | Malone, Minasny, & Brungard (2019) |
sample_existing() |
Sub-sampling an existing sample |
sample_srs
We have demonstrated a simple example of using the
sample_srs()
function in vignette("sgsR")
. We
will demonstrate additional examples below.
raster
The input required for sample_srs()
is a
raster
. This means that sraster
and
mraster
are supported for this function.
#--- perform simple random sampling ---#
sample_srs(
raster = sraster, # input sraster
nSamp = 200, # number of desired sample units
plot = TRUE
) # plot
#> Simple feature collection with 200 features and 0 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431130 ymin: 5337730 xmax: 438550 ymax: 5343230
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#> geometry
#> 1 POINT (431330 5341590)
#> 2 POINT (437250 5340990)
#> 3 POINT (437930 5341210)
#> 4 POINT (434750 5337990)
#> 5 POINT (432350 5339310)
#> 6 POINT (435730 5342650)
#> 7 POINT (434330 5340490)
#> 8 POINT (433490 5342290)
#> 9 POINT (434310 5339690)
#> 10 POINT (431170 5342990)
sample_srs(
raster = mraster, # input mraster
nSamp = 200, # number of desired sample units
access = access, # define access road network
mindist = 200, # minimum distance sample units must be apart from one another
buff_inner = 50, # inner buffer - no sample units within this distance from road
buff_outer = 200, # outer buffer - no sample units further than this distance from road
plot = TRUE
) # plot
#> Simple feature collection with 200 features and 0 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431130 ymin: 5337770 xmax: 438550 ymax: 5343230
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#> geometry
#> 1 POINT (438550 5338470)
#> 2 POINT (438010 5339650)
#> 3 POINT (432510 5342190)
#> 4 POINT (431850 5343050)
#> 5 POINT (434050 5343190)
#> 6 POINT (438130 5341530)
#> 7 POINT (437690 5339250)
#> 8 POINT (437910 5338990)
#> 9 POINT (437970 5339390)
#> 10 POINT (435930 5339050)
sample_systematic
The sample_systematic()
function applies systematic
sampling across an area with the cellsize
parameter
defining the resolution of the tessellation. The tessellation shape can
be modified using the square
parameter. Assigning
TRUE
(default) to the square
parameter results
in a regular grid and assigning FALSE
results in a
hexagonal grid.
The location of sample units can also be adjusted using the
locations
parameter, where centers
takes the
center, corners
takes all corners, and random
takes a random location within each tessellation. Random start points
and translations are applied when the function is called.
#--- perform grid sampling ---#
sample_systematic(
raster = sraster, # input sraster
cellsize = 1000, # grid distance
plot = TRUE
) # plot
#> Simple feature collection with 39 features and 0 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431130.3 ymin: 5337824 xmax: 438296.4 ymax: 5343198
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#> geometry
#> 1 POINT (437838.4 5342760)
#> 2 POINT (438289.7 5341420)
#> 3 POINT (437393.9 5341865)
#> 4 POINT (436498.1 5342309)
#> 5 POINT (435602.4 5342754)
#> 6 POINT (434706.6 5343198)
#> 7 POINT (437845.2 5340524)
#> 8 POINT (436053.6 5341413)
#> 9 POINT (435157.9 5341858)
#> 10 POINT (434262.1 5342302)
#--- perform grid sampling ---#
sample_systematic(
raster = sraster, # input sraster
cellsize = 500, # grid distance
square = FALSE, # hexagonal tessellation
location = "random", # randomly sample within tessellation
plot = TRUE
) # plot
#> Simple feature collection with 161 features and 0 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431248 ymin: 5337746 xmax: 438553.6 ymax: 5343173
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#> geometry
#> 1 POINT (432527.9 5337829)
#> 2 POINT (431687.4 5337982)
#> 3 POINT (435540.2 5337755)
#> 4 POINT (434756.5 5337792)
#> 5 POINT (432895.7 5338116)
#> 6 POINT (432013.1 5338033)
#> 7 POINT (431475.6 5338141)
#> 8 POINT (438553.6 5337814)
#> 9 POINT (437662.6 5337746)
#> 10 POINT (436998.8 5337992)
sample_systematic(
raster = sraster, # input sraster
cellsize = 500, # grid distance
access = access, # define access road network
buff_outer = 200, # outer buffer - no sample units further than this distance from road
square = FALSE, # hexagonal tessellation
location = "corners", # take corners instead of centers
plot = TRUE
)
#> Simple feature collection with 629 features and 0 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431302 ymin: 5337731 xmax: 438483.5 ymax: 5343198
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#> geometry
#> 1 POINT (431863.3 5337798)
#> 2 POINT (431863.3 5337798)
#> 3 POINT (431769.8 5338289)
#> 4 POINT (431958.3 5338070)
#> 5 POINT (431863.3 5337798)
#> 6 POINT (431863.3 5337798)
#> 7 POINT (431958.3 5338070)
#> 8 POINT (431769.8 5338289)
#> 9 POINT (431864.8 5338561)
#> 10 POINT (431769.8 5338289)
sample_strat
The sample_strat()
contains two method
s to
perform sampling:
"Queinnec"
- Hierarchical sampling using a focal
window to isolate contiguous groups of stratum pixels, which was
originally developed by Martin Queinnec.
"random"
- Traditional stratified random sampling.
This method
ignores much of the functionality of the
algorithm to allow users the capability to use standard stratified
random sampling approaches without the use of a focal window to locate
contiguous stratum cells.
method = "Queinnec"
Queinnec, M., White, J. C., & Coops, N. C. (2021). Comparing airborne and spaceborne photon-counting LiDAR canopy structural estimates across different boreal forest types. Remote Sensing of Environment, 262(August 2020), 112510.
This algorithm uses moving window (wrow
and
wcol
parameters) to filter the input sraster
to prioritize sample unit allocation to where stratum pixels are
spatially grouped, rather than dispersed individuals across the
landscape.
Sampling is performed using 2 rules:
Rule 1 - Sample within spatially grouped stratum
pixels. Moving window defined by wrow
and
wcol
.
Rule 2 - If no additional sample units exist to
satisfy desired sample size(nSamp
), individual stratum
pixels are sampled.
The rule applied to a select each sample unit is defined in the
rule
attribute of output samples. We give a few examples
below:
#--- perform stratified sampling random sampling ---#
sample_strat(
sraster = sraster, # input sraster
nSamp = 200
) # desired sample size # plot
#> Simple feature collection with 200 features and 3 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431390 ymin: 5337770 xmax: 438510 ymax: 5343070
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#> strata type rule geometry
#> x 1 new rule1 POINT (435650 5339310)
#> x1 1 new rule1 POINT (433930 5341230)
#> x2 1 new rule1 POINT (437910 5338810)
#> x3 1 new rule1 POINT (433750 5341130)
#> x4 1 new rule1 POINT (433550 5341310)
#> x5 1 new rule1 POINT (435970 5339410)
#> x6 1 new rule1 POINT (437730 5342550)
#> x7 1 new rule1 POINT (433910 5341330)
#> x8 1 new rule1 POINT (434710 5342370)
#> x9 1 new rule1 POINT (434370 5342530)
In some cases, users might want to include an existing
sample within the algorithm. In order to adjust the total number of
sample units needed per stratum to reflect those already present in
existing
, we can use the intermediate function
extract_strata()
.
This function uses the sraster
and existing
sample units and extracts the stratum for each. These sample units can
be included within sample_strat()
, which adjusts total
sample units required per class based on representation in
existing
.
#--- extract strata values to existing samples ---#
e.sr <- extract_strata(
sraster = sraster, # input sraster
existing = existing
) # existing samples to add strata value to
TIP!
sample_strat()
requires the sraster
input
to have an attribute named strata
and will give an error if
it doesn’t.
sample_strat(
sraster = sraster, # input sraster
nSamp = 200, # desired sample size
access = access, # define access road network
existing = e.sr, # existing sample with strata values
mindist = 200, # minimum distance sample units must be apart from one another
buff_inner = 50, # inner buffer - no sample units within this distance from road
buff_outer = 200, # outer buffer - no sample units further than this distance from road
plot = TRUE
) # plot
#> Simple feature collection with 400 features and 3 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431130 ymin: 5337730 xmax: 438530 ymax: 5343230
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#> strata type rule geometry
#> 1 1 existing existing POINT (433430 5340850)
#> 2 1 existing existing POINT (435610 5342270)
#> 3 1 existing existing POINT (438450 5339510)
#> 4 1 existing existing POINT (433150 5341670)
#> 5 1 existing existing POINT (433130 5341450)
#> 6 1 existing existing POINT (433610 5340730)
#> 7 1 existing existing POINT (433830 5340530)
#> 8 1 existing existing POINT (436030 5342910)
#> 9 1 existing existing POINT (438490 5338290)
#> 10 1 existing existing POINT (437790 5338830)
The code in the example above defined the mindist
parameter, which specifies the minimum euclidean distance that new
sample units must be apart from one another.
Notice that the sample units have type
and
rule
attributes which outline whether they are
existing
or new
, and whether
rule1
or rule2
were used to select them. If
type
is existing (a user provided
existing
sample), rule
will be
existing as well as seen above.
sample_strat(
sraster = sraster, # input
nSamp = 200, # desired sample size
access = access, # define access road network
existing = e.sr, # existing samples with strata values
include = TRUE, # include existing sample in nSamp total
buff_outer = 200, # outer buffer - no samples further than this distance from road
plot = TRUE
) # plot
#> Simple feature collection with 200 features and 3 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431130 ymin: 5337750 xmax: 438530 ymax: 5343210
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#> strata type rule geometry
#> 1 1 existing existing POINT (433430 5340850)
#> 2 1 existing existing POINT (435610 5342270)
#> 3 1 existing existing POINT (438450 5339510)
#> 4 1 existing existing POINT (433150 5341670)
#> 5 1 existing existing POINT (433130 5341450)
#> 6 1 existing existing POINT (433610 5340730)
#> 7 1 existing existing POINT (433830 5340530)
#> 8 1 existing existing POINT (436030 5342910)
#> 9 1 existing existing POINT (438490 5338290)
#> 10 1 existing existing POINT (437790 5338830)
The include
parameter determines whether
existing
sample units should be included in the total
sample size defined by nSamp
. By default, the
include
parameter is set as FALSE
.
method = "random
Stratified random sampling with equal probability for all cells
(using default algorithm values for mindist
and no use of
access
functionality). In essence this method perform the
sample_srs
algorithm for each stratum separately to meet
the specified sample size.
#--- perform stratified sampling random sampling ---#
sample_strat(
sraster = sraster, # input sraster
method = "random", # stratified random sampling
nSamp = 200, # desired sample size
plot = TRUE
) # plot
#> Simple feature collection with 200 features and 2 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431150 ymin: 5337710 xmax: 438550 ymax: 5343210
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#> strata type geometry
#> x 1 new POINT (438110 5338610)
#> x1 1 new POINT (438070 5343090)
#> x2 1 new POINT (438010 5343110)
#> x3 1 new POINT (435310 5342470)
#> x4 1 new POINT (438050 5341790)
#> x5 1 new POINT (437470 5343090)
#> x6 1 new POINT (438370 5338210)
#> x7 1 new POINT (437410 5337750)
#> x8 1 new POINT (434130 5341250)
#> x9 1 new POINT (434650 5341390)
sample_sys_strat
sample_sys_strat()
function implements systematic
stratified sampling on an sraster
. This function uses the
same functionality as sample_systematic()
but takes an
sraster
as input and performs sampling on each stratum
iteratively.
#--- perform grid sampling on each stratum separately ---#
sample_sys_strat(
sraster = sraster, # input sraster with 4 strata
cellsize = 1000, # grid size
plot = TRUE # plot output
)
#> Warning: [readStart] source already open for reading
#> Processing strata : 1
#> Warning: [extract] source already open for reading
#> Processing strata : 2
#> Warning: [extract] source already open for reading
#> Processing strata : 3
#> Warning: [extract] source already open for reading
#> Processing strata : 4
#> Warning: [extract] source already open for reading
#> Simple feature collection with 28 features and 1 field
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431394.4 ymin: 5337720 xmax: 438426 ymax: 5342671
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#> strata geometry
#> 1 1 POINT (433959.2 5338974)
#> 2 1 POINT (433850.1 5341208)
#> 3 1 POINT (436683 5338212)
#> 4 1 POINT (435593.2 5342188)
#> 5 1 POINT (438044.8 5337831)
#> 6 1 POINT (438426 5339192)
#> 7 1 POINT (437826.6 5342297)
#> 8 2 POINT (437623.6 5338249)
#> 9 2 POINT (437775.4 5341245)
#> 10 2 POINT (434678 5339400)
Just like with sample_systematic()
we can specify where
we want our samples to fall within our tessellations. We specify
location = "corners"
below. Note that the tesselations are
all saved to a list file when details = TRUE
should the
user want to save them.
sample_sys_strat(
sraster = sraster, # input sraster with 4 strata
cellsize = 500, # grid size
square = FALSE, # hexagon tessellation
location = "corners", # samples on tessellation corners
plot = TRUE # plot output
)
#> Processing strata : 1
#> Warning: [extract] source already open for reading
#> Processing strata : 2
#> Warning: [extract] source already open for reading
#> Processing strata : 3
#> Warning: [extract] source already open for reading
#> Processing strata : 4
#> Warning: [extract] source already open for reading
#> Simple feature collection with 1295 features and 1 field
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431123.3 ymin: 5337731 xmax: 438526.7 ymax: 5343234
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#> strata geometry
#> 1 1 POINT (438164.4 5337823)
#> 2 1 POINT (438164.4 5337823)
#> 3 1 POINT (437413.3 5337961)
#> 4 1 POINT (437698.7 5338005)
#> 5 1 POINT (436842.5 5337875)
#> 6 1 POINT (435986.3 5337745)
#> 7 1 POINT (438164.4 5337823)
#> 8 1 POINT (437698.7 5338005)
#> 9 1 POINT (437803.9 5338274)
#> 10 1 POINT (438089.3 5338317)
This sampling approach could be especially useful incombination with
strat_poly()
to ensure consistency of sampling accross
specific management units.
#--- read polygon coverage ---#
poly <- system.file("extdata", "inventory_polygons.shp", package = "sgsR")
fri <- sf::st_read(poly)
#> Reading layer `inventory_polygons' from data source
#> `/tmp/Rtmp4h4X3Y/Rinst152a1db126c5/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
#--- stratify polygon coverage ---#
#--- 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")
#--- get polygon stratification ---#
srasterpoly <- strat_poly(
poly = fri,
attribute = attribute,
features = features,
raster = sraster
)
#--- systematatic stratified sampling for each stratum ---#
sample_sys_strat(
sraster = srasterpoly, # input sraster from strat_poly() with 3 strata
cellsize = 500, # grid size
square = FALSE, # hexagon tessellation
location = "random", # randomize plot location
plot = TRUE # plot output
)
#> Processing strata : 1
#> Warning: [extract] source already open for reading
#> Processing strata : 2
#> Warning: [extract] source already open for reading
#> Processing strata : 3
#> Warning: [extract] source already open for reading
#> Simple feature collection with 167 features and 1 field
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431103.8 ymin: 5337754 xmax: 438550.7 ymax: 5343192
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#> strata geometry
#> 1 1 POINT (432444.2 5338991)
#> 2 1 POINT (432289.8 5339604)
#> 3 1 POINT (432029.4 5340588)
#> 4 1 POINT (431917.5 5341362)
#> 5 1 POINT (432356.1 5339195)
#> 6 1 POINT (432388.9 5340062)
#> 7 1 POINT (432042 5341746)
#> 8 1 POINT (432514.1 5339964)
#> 9 1 POINT (431999.2 5342180)
#> 10 1 POINT (433171.1 5338495)
sample_nc
sample_nc()
function implements the Nearest Centroid
sampling algorithm described in Melville &
Stone (2016). The algorithm uses kmeans clustering where the number
of clusters (centroids) is equal to the desired sample size
(nSamp
).
Cluster centers are located, which then prompts the nearest neighbour
mraster
pixel for each cluster to be selected (assuming
default k
parameter). These nearest neighbours are the
output sample units.
#--- perform simple random sampling ---#
sample_nc(
mraster = mraster, # input
nSamp = 25, # desired sample size
plot = TRUE
)
#> K-means being performed on 3 layers with 25 centers.
#> Simple feature collection with 25 features and 4 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431230 ymin: 5337890 xmax: 438510 ymax: 5343230
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#> zq90 pzabove2 zsd kcenter geometry
#> 88737 8.73 40.5 2.33 1 POINT (437810 5338490)
#> 83289 23.80 89.8 6.83 2 POINT (433290 5338770)
#> 44181 9.92 66.5 2.57 3 POINT (434430 5340870)
#> 57393 26.50 85.7 8.38 4 POINT (437570 5340170)
#> 59102 21.90 92.3 5.25 5 POINT (434450 5340070)
#> 51993 20.60 85.5 6.15 6 POINT (434010 5340450)
#> 6226 6.04 59.4 1.43 7 POINT (436250 5342910)
#> 87111 17.70 92.6 3.47 8 POINT (435130 5338570)
#> 36188 13.00 73.5 3.33 9 POINT (431230 5341290)
#> 94607 12.10 53.1 3.27 10 POINT (435850 5338170)
Altering the k
parameter leads to a multiplicative
increase in output sample units where total output samples = nSamp * k.
#--- perform simple random sampling ---#
samples <- sample_nc(
mraster = mraster, # input
k = 2, # number of nearest neighbours to take for each kmeans center
nSamp = 25, # desired sample size
plot = TRUE
)
#> K-means being performed on 3 layers with 25 centers.
Visualizing what the kmeans centers and sample units looks like is
possible when using details = TRUE
. The $kplot
output provides a quick visualization of where the centers are based on
a scatter plot of the first 2 layers in mraster
. Notice
that the centers are well distributed in covariate space and chosen
sample units are the closest pixels to each center (nearest
neighbours).
#--- perform simple random sampling with details ---#
details <- sample_nc(
mraster = mraster, # input
nSamp = 25, # desired sample number
details = TRUE
)
#> K-means being performed on 3 layers with 25 centers.
#--- plot ggplot output ---#
details$kplot
sample_clhs
sample_clhs()
function implements conditioned Latin
hypercube (clhs) sampling methodology from the clhs
package.
TIP!
A number of other functions in the sgsR
package help to
provide guidance on clhs sampling including calculate_pop()
and calculate_lhsOpt()
. Check out these functions to better
understand how sample numbers could be optimized.
The syntax for this function is similar to others shown above,
although parameters like iter
, which define the number of
iterations within the Metropolis-Hastings process are important to
consider. In these examples we use a low iter
value for
efficiency. Default values for iter
within the
clhs
package are 10,000.
sample_clhs(
mraster = mraster, # input
nSamp = 200, # desired sample size
plot = TRUE, # plot
iter = 100
) # number of iterations
The cost
parameter defines the mraster
covariate, which is used to constrain the clhs sampling. An example
could be the distance a pixel is from road access
(e.g. from calculate_distance()
see example below), terrain
slope, the output from calculate_coobs()
, or many
others.
#--- cost constrained examples ---#
#--- calculate distance to access layer for each pixel in mr ---#
mr.c <- calculate_distance(
raster = mraster, # input
access = access, # define access road network
plot = TRUE
) # plot
sample_clhs(
mraster = mr.c, # input
nSamp = 250, # desired sample size
iter = 100, # number of iterations
cost = "dist2access", # cost parameter - name defined in calculate_distance()
plot = TRUE
) # plot
sample_balanced
The sample_balanced()
algorithm performs a balanced
sampling methodology from the stratifyR / SamplingBigData
packages.
#> Simple feature collection with 200 features and 0 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431110 ymin: 5337710 xmax: 438550 ymax: 5343190
#> Projected CRS: +proj=utm +zone=17 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
#> First 10 features:
#> geometry
#> 1 POINT (431750 5343190)
#> 2 POINT (431770 5343190)
#> 3 POINT (433350 5343150)
#> 4 POINT (434350 5343150)
#> 5 POINT (435230 5343090)
#> 6 POINT (431570 5343030)
#> 7 POINT (432250 5342970)
#> 8 POINT (435150 5342950)
#> 9 POINT (433450 5342930)
#> 10 POINT (436210 5342890)
sample_balanced(
mraster = mraster, # input
nSamp = 100, # desired sample size
algorithm = "lcube", # algorithm type
access = access, # define access road network
buff_inner = 50, # inner buffer - no sample units within this distance from road
buff_outer = 200
) # outer buffer - no sample units further than this distance from road
#> Simple feature collection with 100 features and 0 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431350 ymin: 5337850 xmax: 438510 ymax: 5343190
#> Projected CRS: +proj=utm +zone=17 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
#> First 10 features:
#> geometry
#> 1 POINT (434170 5343190)
#> 2 POINT (433110 5343170)
#> 3 POINT (435510 5343170)
#> 4 POINT (432510 5343090)
#> 5 POINT (433590 5343050)
#> 6 POINT (435430 5343050)
#> 7 POINT (435350 5343030)
#> 8 POINT (434070 5343010)
#> 9 POINT (432610 5342990)
#> 10 POINT (433370 5342910)
sample_ahels
The sample_ahels()
function performs the adapted
Hypercube Evaluation of a Legacy Sample (ahels) algorithm
usingexisting
sample data and an mraster
. New
sample units are allocated based on quantile ratios between the
existing
sample and mraster
covariate
dataset.
This algorithm was adapted from that presented in the paper below, which we highly recommend.
Malone BP, Minansy B, Brungard C. 2019. Some methods to improve the utility of conditioned Latin hypercube sampling. PeerJ 7:e6451 DOI 10.7717/peerj.6451
This algorithm:
Determines the quantile distributions of existing
sample units and mraster
covariates.
Determines quantiles where there is a disparity between sample units and covariates.
Prioritizes sampling within those quantile to improve representation.
To use this function, user must first specify the number of quantiles
(nQuant
) followed by either the nSamp
(total
number of desired sample units to be added) or the
threshold
(sampling ratio vs. covariate coverage ratio for
quantiles - default is 0.9) parameters.
#--- remove `type` variable from existing - causes plotting issues ---#
existing <- existing %>% select(-type)
sample_ahels(
mraster = mraster,
existing = existing, # existing sample
plot = TRUE
) # plot
#> Simple feature collection with 346 features and 7 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431110 ymin: 5337750 xmax: 438530 ymax: 5343210
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#> type.x zq90 pzabove2 zsd strata type.y rule geometry
#> 1 existing 3.41 56.1 0.60 1 new rule1 POINT (433430 5340850)
#> 2 existing 9.53 76.3 2.34 1 new rule1 POINT (435610 5342270)
#> 3 existing 3.03 2.8 0.52 1 new rule1 POINT (438450 5339510)
#> 4 existing 4.12 2.2 0.89 1 new rule1 POINT (433150 5341670)
#> 5 existing 6.04 13.5 1.67 1 new rule1 POINT (433130 5341450)
#> 6 existing 4.08 32.2 0.80 1 new rule1 POINT (433610 5340730)
#> 7 existing 6.65 40.3 2.33 1 new rule1 POINT (433830 5340530)
#> 8 existing 4.55 24.9 1.00 1 new rule1 POINT (436030 5342910)
#> 9 existing 2.51 3.6 0.37 1 new rule1 POINT (438490 5338290)
#> 10 existing 2.86 24.8 0.47 1 new rule1 POINT (437790 5338830)
TIP!
Notice that no threshold
, nSamp
, or
nQuant
were defined. That is because the default setting
for threshold = 0.9
and nQuant = 10
.
The first matrix output shows the quantile ratios between the sample and the covariates. A value of 1.0 indicates that the sample is representative of quantile coverage. Values > 1.0 indicate over representation of sample units, while < 1.0 indicate under representation.
sample_ahels(
mraster = mraster,
existing = existing, # existing sample
nQuant = 20, # define 20 quantiles
nSamp = 300
) # desired sample size
#> Simple feature collection with 500 features and 7 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431110 ymin: 5337710 xmax: 438530 ymax: 5343230
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#> type.x zq90 pzabove2 zsd strata type.y rule geometry
#> 1 existing 3.41 56.1 0.60 1 new rule1 POINT (433430 5340850)
#> 2 existing 9.53 76.3 2.34 1 new rule1 POINT (435610 5342270)
#> 3 existing 3.03 2.8 0.52 1 new rule1 POINT (438450 5339510)
#> 4 existing 4.12 2.2 0.89 1 new rule1 POINT (433150 5341670)
#> 5 existing 6.04 13.5 1.67 1 new rule1 POINT (433130 5341450)
#> 6 existing 4.08 32.2 0.80 1 new rule1 POINT (433610 5340730)
#> 7 existing 6.65 40.3 2.33 1 new rule1 POINT (433830 5340530)
#> 8 existing 4.55 24.9 1.00 1 new rule1 POINT (436030 5342910)
#> 9 existing 2.51 3.6 0.37 1 new rule1 POINT (438490 5338290)
#> 10 existing 2.86 24.8 0.47 1 new rule1 POINT (437790 5338830)
Notice that the total number of samples is 500. This value is the sum
of existing units (200) and number of sample units defined by
nSamp = 300
.
sample_existing
Acknowledging that existing
sample networks are common
is important. There is significant investment into these samples, and in
order to keep inventories up-to-date, we often need to collect new data
for sample units. The sample_existing
algorithm provides
the user with methods for sub-sampling an existing
sample
network should the financial / logistical resources not be available to
collect data at all sample units. The functions allows users to choose
between algorithm types using (type = "clhs"
- default,
type = "balanced"
, type = "srs"
,
type = "strat"
). Differences in type result in calling
internal sample_existing_*()
functions
(sample_existing_clhs()
(default),
sample_existing_balanced()
,
sample_existing_srs()
,
sample_existing_strat()
). These functions are not exported
to be used stand-alone, however they employ the same functionality as
their sample_clhs()
etc counterparts.
While using sample_existing()
, should the user wish to
specify algorithm specific parameters
(e.g. algorithm = "lcube"
in sample_balanced()
or allocation = "equal"
in sample_strat()
),
they can specify within sample_existing()
as if calling the
function directly.
I give applied examples for all methods below that are based on the following scenario:
We have a systematic sample where sample units are 200m apart.
We know we only have resources to sample 300 of them.
We have some ALS data available (mraster
), which we
can use to improve knowledge of the metric populations.
See our existing
sample for the scenario below.
#--- generate existing samples and extract metrics ---#
existing <- sample_systematic(raster = mraster, cellsize = 200, plot = TRUE)
sample_existing(type = "clhs")
The algorithm is unique in that it has two fundamental approaches:
existing
and the attributes it
contains.#--- sub sample using ---#
sample_existing(existing = e, nSamp = 300, type = "clhs")
#> Simple feature collection with 300 features and 3 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431141.2 ymin: 5337704 xmax: 438554 ymax: 5343237
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#> zq90 pzabove2 zsd geometry
#> 770 21.60 96.1 5.62 POINT (433303 5338207)
#> 689 10.80 55.9 2.91 POINT (435671 5339352)
#> 474 14.30 96.4 2.12 POINT (432583.4 5339869)
#> 610 15.60 89.2 3.51 POINT (436509.6 5340028)
#> 678 8.07 72.3 1.98 POINT (437965.3 5340057)
#> 321 20.10 89.7 6.53 POINT (435598.8 5341632)
#> 4 16.90 85.1 4.02 POINT (431585.4 5343119)
#> 177 4.01 53.6 0.74 POINT (434599.2 5342162)
#> 712 22.70 94.7 4.55 POINT (438215.2 5339924)
#> 614 24.70 88.9 5.11 POINT (435744.8 5339794)
raster
distributionsOur systematic sample of ~900 plots is fairly comprehensive, however
we can generate a true population distribution through the inclusion of
the ALS metrics in the sampling process. The metrics will be included in
internal latin hypercube sampling to help guide sub-sampling of
existing
.
#--- sub sample using ---#
sample_existing(
existing = existing, # our existing sample
nSamp = 300, # desired sample size
raster = mraster, # include mraster metrics to guide sampling of existing
plot = TRUE
) # plot
#> Simple feature collection with 300 features and 3 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431127.7 ymin: 5337704 xmax: 438554 ymax: 5343235
#> Projected CRS: +proj=utm +zone=17 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
#> First 10 features:
#> zq90 pzabove2 zsd geometry
#> 91959 19.20 86.7 5.42 POINT (434641.3 5338618)
#> 91951 0.00 0.0 0.00 POINT (436362 5339146)
#> 91520 4.11 27.8 0.87 POINT (434642.9 5341338)
#> 91608 19.60 96.3 5.04 POINT (431127.7 5339840)
#> 92060 13.20 76.0 3.77 POINT (437229.1 5338366)
#> 91596 15.60 45.9 4.22 POINT (433613.2 5340604)
#> 91212 17.00 65.4 4.76 POINT (432467.6 5342972)
#> 91920 20.80 25.4 6.82 POINT (435347.4 5339044)
#> 91518 11.90 90.7 2.83 POINT (435025.3 5341456)
#> 91244 19.10 84.1 5.01 POINT (432585 5342589)
The sample distribution again mimics the population distribution quite well! Now lets try using a cost variable to constrain the sub-sample.
#--- create distance from roads metric ---#
dist <- calculate_distance(raster = mraster, access = access)
#--- sub sample using ---#
sample_existing(
existing = existing, # our existing sample
nSamp = 300, # desired sample size
raster = dist, # include mraster metrics to guide sampling of existing
cost = 4, # either provide the index (band number) or the name of the cost layer
plot = TRUE
) # plot
#> Simple feature collection with 300 features and 4 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431127.7 ymin: 5337778 xmax: 438525.4 ymax: 5343237
#> Projected CRS: +proj=utm +zone=17 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
#> First 10 features:
#> zq90 pzabove2 zsd dist2access geometry
#> 91228 11.20 67.3 2.63 361.44437 POINT (432335.1 5342722)
#> 91947 9.06 34.7 2.40 372.30853 POINT (437318 5339440)
#> 92101 1.96 0.0 0.20 7.46457 POINT (437978.8 5337969)
#> 91906 20.00 94.9 4.12 327.78388 POINT (438406.4 5339983)
#> 91197 20.30 75.7 6.20 38.31340 POINT (431967.8 5343237)
#> 91959 19.20 86.7 5.42 227.77954 POINT (434641.3 5338618)
#> 91561 6.40 21.5 1.67 91.03501 POINT (433745.6 5340853)
#> 91271 13.30 86.6 3.05 602.24629 POINT (431305.4 5341987)
#> 92068 18.40 88.7 3.52 309.33632 POINT (435699.6 5337897)
#> 91724 22.40 96.4 4.77 197.93306 POINT (434421.6 5340015)
Finally, should the user wish to further constrain the sample based
on access
like other sampling approaches in
sgsR
that is also possible.
#--- ensure access and existing are in the same CRS ---#
sf::st_crs(existing) <- sf::st_crs(access)
#--- sub sample using ---#
sample_existing(
existing = existing, # our existing sample
nSamp = 300, # desired sample size
raster = dist, # include mraster metrics to guide sampling of existing
cost = 4, # either provide the index (band number) or the name of the cost layer
access = access, # roads layer
buff_inner = 50, # inner buffer - no sample units within this distance from road
buff_outer = 300, # outer buffer - no sample units further than this distance from road
plot = TRUE
) # plot
#> Simple feature collection with 300 features and 4 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431142.8 ymin: 5337704 xmax: 438508.7 ymax: 5343235
#> Projected CRS: +proj=utm +zone=17 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
#> First 10 features:
#> zq90 pzabove2 zsd dist2access geometry
#> 91552 9.79 57.4 2.7400000 179.05833 POINT (434715.1 5339059)
#> 91558 22.70 94.7 4.5499997 141.02542 POINT (438215.2 5339924)
#> 91371 16.00 86.8 4.5299997 84.68565 POINT (431142.8 5340473)
#> 91404 20.40 87.2 5.8799996 152.89818 POINT (431833.7 5340266)
#> 91573 7.67 69.1 1.8099999 164.80159 POINT (437891.5 5339616)
#> 91333 13.70 77.0 3.1799998 237.86991 POINT (437966.9 5342777)
#> 91496 25.80 88.9 7.3999996 80.76568 POINT (438362.8 5340807)
#> 91196 20.40 86.9 4.1500001 98.29372 POINT (431776.6 5343178)
#> 91509 4.46 4.1 0.9899999 230.12394 POINT (438421.5 5340615)
#> 91531 9.33 85.0 2.4800000 103.35894 POINT (436377.1 5339778)
TIP!
The greater constraints we add to sampling, the less likely we will have strong correlations between the population and sample, so its always important to understand these limitations and plan accordingly.
sample_existing(type = "balanced")
When type = "balanced"
users can define all parameters
that are found within sample_balanced()
. This means that
one can change the algorithm
, p
etc.
sample_existing(existing = e, nSamp = 300, type = "balanced")
#> Simple feature collection with 300 features and 3 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431141.2 ymin: 5337706 xmax: 438538.9 ymax: 5343235
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#> zq90 pzabove2 zsd geometry
#> 7 16.40 80.5 4.36 POINT (432217.7 5343104)
#> 10 22.10 89.5 4.71 POINT (431644.1 5342928)
#> 11 19.80 81.2 6.57 POINT (431452.9 5342869)
#> 13 12.80 48.2 3.31 POINT (433232.4 5343206)
#> 14 17.70 75.2 4.74 POINT (433041.2 5343148)
#> 15 5.23 53.3 1.03 POINT (432850 5343089)
#> 29 17.30 89.0 4.60 POINT (433099.9 5342957)
#> 30 10.90 75.2 2.62 POINT (432908.7 5342898)
#> 31 9.62 70.3 2.17 POINT (432717.5 5342839)
#> 35 13.20 79.1 2.88 POINT (431952.7 5342604)
sample_existing(existing = e, nSamp = 300, type = "balanced", algorithm = "lcube")
#> Simple feature collection with 300 features and 3 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431112.6 ymin: 5337764 xmax: 438552.4 ymax: 5343205
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#> zq90 pzabove2 zsd geometry
#> 8 23.90 82.3 7.14 POINT (432026.5 5343045)
#> 9 18.30 91.5 3.32 POINT (431835.3 5342987)
#> 10 22.10 89.5 4.71 POINT (431644.1 5342928)
#> 12 19.30 87.1 5.65 POINT (431261.8 5342811)
#> 15 5.23 53.3 1.03 POINT (432850 5343089)
#> 18 13.80 65.5 3.46 POINT (432276.4 5342913)
#> 23 17.60 46.5 5.89 POINT (431320.5 5342619)
#> 28 18.80 88.2 3.69 POINT (433291.1 5343015)
#> 29 17.30 89.0 4.60 POINT (433099.9 5342957)
#> 32 10.10 55.8 2.37 POINT (432526.3 5342780)
sample_existing(type = "srs")
The simplest, type = srs
, randomly selects sample
units.
sample_existing(existing = e, nSamp = 300, type = "srs")
#> Simple feature collection with 300 features and 3 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431112.6 ymin: 5337707 xmax: 438540.5 ymax: 5343220
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#> zq90 pzabove2 zsd geometry
#> 1 1.42 0.000000 0.05 POINT (437717 5342910)
#> 2 14.60 61.799999 4.34 POINT (432393.8 5342530)
#> 3 16.20 94.800003 3.03 POINT (432920.6 5338089)
#> 4 25.80 88.900002 7.40 POINT (438362.8 5340807)
#> 5 2.51 9.400001 0.37 POINT (437584.5 5342660)
#> 6 17.50 70.900002 4.88 POINT (433452.1 5341809)
#> 7 14.70 72.700005 3.44 POINT (433187.1 5341310)
#> 8 6.93 30.800001 1.80 POINT (433628.2 5341236)
#> 9 7.42 41.600002 1.85 POINT (437611.5 5338484)
#> 10 18.10 94.099998 3.34 POINT (436511.1 5342749)
sample_existing(type = "strat")
When type = "strat"
, existing
must have an
attribute named strata
(just like how
sample_strat()
requires a strata
layer). If it
doesnt exist you will get an error. Lets define an sraster
so that we are compliant.
sraster <- strat_kmeans(mraster = mraster, nStrata = 4)
e_strata <- extract_strata(sraster = sraster, existing = e)
When we do have a strata attribute, the function works very much the
same as sample_strat()
in that is allows the user to define
the allocation
method ("prop"
- defaults,
"optim"
, "manual"
, "equal"
).
#--- proportional stratified sampling of existing ---#
sample_existing(existing = e_strata, nSamp = 300, type = "strat", allocation = "prop")
#> Simple feature collection with 301 features and 4 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431129.3 ymin: 5337706 xmax: 438552.4 ymax: 5343237
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#> strata zq90 pzabove2 zsd geometry
#> 1 4 17.8 96.0 3.78 POINT (434877.7 5340573)
#> 2 4 16.9 96.1 2.08 POINT (434025.7 5341985)
#> 3 4 15.2 95.8 3.07 POINT (437202.1 5342542)
#> 4 4 13.9 96.0 2.94 POINT (434289.1 5339765)
#> 5 4 17.3 89.0 4.60 POINT (433099.9 5342957)
#> 6 4 19.2 93.7 4.00 POINT (437333 5340072)
#> 7 4 14.6 93.2 2.84 POINT (432700.8 5339487)
#> 8 4 16.2 96.8 3.30 POINT (432011.4 5342413)
#> 9 4 15.0 91.1 2.82 POINT (433422 5340545)
#> 10 4 14.2 87.3 3.39 POINT (437819.3 5341895)
TIP!
Remember that when allocation = "equal"
, the
nSamp
value will be allocated for each strata.
We get 400 sample units in our output below because we have 4 strata
and nSamp = 100
.
#--- equal stratified sampling of existing ---#
sample_existing(existing = e_strata, nSamp = 100, type = "strat", allocation = "equal")
#> Simple feature collection with 400 features and 4 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431112.6 ymin: 5337704 xmax: 438552.4 ymax: 5343235
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#> strata zq90 pzabove2 zsd geometry
#> 1 4 16.5 86.4 4.09 POINT (431687.8 5342104)
#> 2 4 17.0 93.8 3.42 POINT (438275.5 5342454)
#> 3 4 14.1 84.2 3.20 POINT (433702 5341677)
#> 4 4 17.2 83.0 4.58 POINT (434496.9 5343176)
#> 5 4 14.1 92.6 3.38 POINT (431614 5341663)
#> 6 4 15.6 55.9 4.55 POINT (434773.8 5338868)
#> 7 4 15.8 89.0 3.75 POINT (431657.6 5340840)
#> 8 4 19.0 86.5 4.22 POINT (436687.3 5342175)
#> 9 4 16.1 80.5 4.37 POINT (435391 5338220)
#> 10 4 18.5 99.4 2.64 POINT (438067.7 5339042)
#--- manual stratified sampling of existing with user defined weights ---#
s <- sample_existing(existing = e_strata, nSamp = 100, type = "strat", allocation = "manual", weights = c(0.2, 0.6, 0.1, 0.1))
We can check the proportion of samples from each strata with:
Finally, type = "optim
allows for the user to define a
raster
metric to be used to optimize within strata
variances.
#--- manual stratified sampling of existing with user defined weights ---#
sample_existing(existing = e_strata, nSamp = 100, type = "strat", allocation = "optim", raster = mraster, metric = "zq90")
#> Simple feature collection with 100 features and 4 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 431156.2 ymin: 5337793 xmax: 438493.7 ymax: 5343163
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#> strata zq90 pzabove2 zsd geometry
#> 1 4 14.3 89.0 3.70 POINT (436452.4 5342940)
#> 2 4 16.4 91.9 3.34 POINT (433391.8 5339280)
#> 3 4 18.1 92.8 3.25 POINT (432318.5 5339369)
#> 4 4 17.1 96.3 2.35 POINT (433834.5 5341927)
#> 5 4 14.0 96.3 3.36 POINT (437405.2 5337793)
#> 6 4 13.0 83.4 3.55 POINT (437730.5 5340822)
#> 7 4 14.0 92.5 2.49 POINT (431156.2 5338385)
#> 8 4 17.6 71.0 4.50 POINT (431848.8 5340899)
#> 9 4 17.4 95.4 3.46 POINT (435053.8 5340000)
#> 10 4 19.3 92.8 2.59 POINT (431716.3 5340649)
We see from the output that we get 300 sample units that are a
sub-sample of existing
. The plotted output shows cumulative
frequency distributions of the population (all existing
samples) and the sub-sample (the 300 samples we requested).