Sampling

sample_srs

We have demonstrated a simple example of using the sample_srs() function in vignette("sgsR"). We will demonstrate additional examples below.

#--- 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: 431250 ymin: 5337710 xmax: 438550 ymax: 5343190
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#>                  geometry
#> 1  POINT (437870 5341770)
#> 2  POINT (436950 5341590)
#> 3  POINT (436690 5341230)
#> 4  POINT (431990 5338030)
#> 5  POINT (435830 5337870)
#> 6  POINT (437550 5337930)
#> 7  POINT (434250 5342030)
#> 8  POINT (435390 5341250)
#> 9  POINT (434850 5342630)
#> 10 POINT (433210 5339590)

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: 5337730 xmax: 438470 ymax: 5343230
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#>                  geometry
#> 1  POINT (435050 5342570)
#> 2  POINT (435530 5341970)
#> 3  POINT (436450 5341350)
#> 4  POINT (438010 5340410)
#> 5  POINT (437750 5338370)
#> 6  POINT (435270 5340230)
#> 7  POINT (435190 5338410)
#> 8  POINT (433590 5340550)
#> 9  POINT (435090 5342910)
#> 10 POINT (435670 5338630)

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 37 features and 0 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 431185.1 ymin: 5337734 xmax: 438460.2 ymax: 5343128
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#>                    geometry
#> 1  POINT (431515.3 5338355)
#> 2    POINT (432466 5338045)
#> 3  POINT (433416.6 5337734)
#> 4  POINT (432776.2 5338995)
#> 5  POINT (434677.5 5338375)
#> 6  POINT (435628.2 5338065)
#> 7  POINT (436578.8 5337754)
#> 8  POINT (431185.1 5340566)
#> 9  POINT (432135.8 5340256)
#> 10 POINT (433086.4 5339946)

#--- 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 165 features and 0 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 431221.3 ymin: 5337701 xmax: 438540.2 ymax: 5343215
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#>                    geometry
#> 1  POINT (431355.9 5343139)
#> 2  POINT (431732.1 5343099)
#> 3  POINT (431326.1 5342690)
#> 4  POINT (432236.3 5342823)
#> 5  POINT (431409.9 5342100)
#> 6  POINT (431972.2 5342482)
#> 7  POINT (432621.3 5343045)
#> 8  POINT (431761.4 5342106)
#> 9  POINT (432372.4 5342734)
#> 10 POINT (433334.1 5343150)

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 623 features and 0 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 431153.9 ymin: 5337701 xmax: 438423.7 ymax: 5343196
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#>                    geometry
#> 1  POINT (438423.7 5338002)
#> 2  POINT (437924.9 5337711)
#> 3  POINT (438423.7 5338002)
#> 4  POINT (437924.9 5337711)
#> 5  POINT (437923.7 5338000)
#> 6  POINT (438173.1 5338145)
#> 7  POINT (438423.7 5338002)
#> 8  POINT (438423.7 5338002)
#> 9  POINT (438173.1 5338145)
#> 10 POINT (438171.9 5338434)

sample_strat

The sample_strat() contains two methods 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.

#--- 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: 431150 ymin: 5337730 xmax: 438530 ymax: 5343210
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#>    strata type  rule               geometry
#> x       1  new rule1 POINT (435650 5339450)
#> x1      1  new rule1 POINT (433570 5340510)
#> x2      1  new rule1 POINT (433070 5343050)
#> x3      1  new rule1 POINT (434710 5341790)
#> x4      1  new rule1 POINT (437410 5338090)
#> x5      1  new rule1 POINT (436050 5342530)
#> x6      1  new rule1 POINT (436770 5338210)
#> x7      1  new rule1 POINT (437150 5338230)
#> x8      1  new rule1 POINT (438050 5339750)
#> x9      1  new rule1 POINT (435490 5342110)

#--- extract strata values to existing samples ---#
e.sr <- extract_strata(
  sraster = sraster, # input sraster
  existing = existing
) # existing samples to add strata value to

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 (434490 5341390)
#> 2       1 existing existing POINT (434970 5341730)
#> 3       1 existing existing POINT (438030 5341090)
#> 4       1 existing existing POINT (436870 5338030)
#> 5       1 existing existing POINT (437230 5337970)
#> 6       1 existing existing POINT (436410 5338390)
#> 7       1 existing existing POINT (434010 5340850)
#> 8       1 existing existing POINT (435590 5342730)
#> 9       1 existing existing POINT (433690 5340530)
#> 10      1 existing existing POINT (435530 5342490)

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: 431170 ymin: 5337730 xmax: 438530 ymax: 5343210
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#>    strata     type     rule               geometry
#> 1       1 existing existing POINT (434490 5341390)
#> 2       1 existing existing POINT (434970 5341730)
#> 3       1 existing existing POINT (438030 5341090)
#> 4       1 existing existing POINT (436870 5338030)
#> 5       1 existing existing POINT (437230 5337970)
#> 6       1 existing existing POINT (436410 5338390)
#> 7       1 existing existing POINT (434010 5340850)
#> 8       1 existing existing POINT (435590 5342730)
#> 9       1 existing existing POINT (433690 5340530)
#> 10      1 existing existing POINT (435530 5342490)

#--- 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: 431130 ymin: 5337770 xmax: 438530 ymax: 5343210
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#>    strata type               geometry
#> x       1  new POINT (437230 5338430)
#> x1      1  new POINT (437870 5342290)
#> x2      1  new POINT (438090 5338750)
#> x3      1  new POINT (434210 5340790)
#> x4      1  new POINT (434470 5342850)
#> x5      1  new POINT (434030 5338810)
#> x6      1  new POINT (437810 5338390)
#> x7      1  new POINT (433590 5340590)
#> x8      1  new POINT (435250 5340430)
#> x9      1  new POINT (435470 5341750)

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 51 features and 1 field
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 431108.5 ymin: 5337716 xmax: 438392.5 ymax: 5343223
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#>    strata                 geometry
#> 1       1 POINT (432313.7 5343223)
#> 2       1 POINT (433712.4 5343014)
#> 3       1 POINT (435111.1 5342806)
#> 4       1 POINT (434307.3 5342211)
#> 5       1 POINT (433503.5 5341616)
#> 6       1   POINT (431896 5340426)
#> 7       1 POINT (434098.5 5340812)
#> 8       1 POINT (437908.5 5342388)
#> 9       1 POINT (434693.4 5340008)
#> 10      1 POINT (437699.7 5340989)

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 1213 features and 1 field
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 431101 ymin: 5337707 xmax: 438557.2 ymax: 5343214
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#>    strata                 geometry
#> 1       1 POINT (431180.5 5340689)
#> 2       1 POINT (431180.5 5340689)
#> 3       1 POINT (431180.5 5340689)
#> 4       1   POINT (431274 5341944)
#> 5       1   POINT (431274 5341944)
#> 6       1   POINT (431274 5341944)
#> 7       1   POINT (431274 5341944)
#> 8       1 POINT (433639.3 5338181)
#> 9       1 POINT (433639.3 5338181)
#> 10      1 POINT (433733.4 5338672)

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/RtmpJ0rxyo/Rinst12f07010d104/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 165 features and 1 field
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 431117.8 ymin: 5337705 xmax: 438553.2 ymax: 5343226
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#>    strata                 geometry
#> 1       1   POINT (431531 5338870)
#> 2       1 POINT (431148.3 5340112)
#> 3       1 POINT (432078.5 5339670)
#> 4       1 POINT (432383.5 5339203)
#> 5       1 POINT (432093.3 5340065)
#> 6       1   POINT (432530 5339923)
#> 7       1   POINT (432118 5340342)
#> 8       1 POINT (433888.6 5338107)
#> 9       1 POINT (431782.3 5341669)
#> 10      1 POINT (433799.5 5338482)

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: 432090 ymin: 5337710 xmax: 438530 ymax: 5343130
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#>        zq90 pzabove2  zsd kcenter               geometry
#> 16935 15.50     94.8 2.71       1 POINT (434090 5342330)
#> 77033 10.70     17.2 3.07       2 POINT (434990 5339110)
#> 46704 19.40     91.8 4.26       3 POINT (432670 5340730)
#> 59233 20.50     83.9 6.02       4 POINT (437070 5340070)
#> 47939 17.80     92.8 3.54       5 POINT (434990 5340670)
#> 98002  8.41     87.0 1.81       6 POINT (436610 5337990)
#> 25200  4.69     33.0 1.03       7 POINT (435270 5341890)
#> 38710 16.60     83.3 4.26       8 POINT (436910 5341170)
#> 63222  6.56     14.0 1.74       9 POINT (434790 5339850)
#> 78491  3.00      7.6 0.58      10 POINT (434310 5339030)

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.

#--- total samples = nSamp * k (25 * 2) = 50 ---#
nrow(samples)
#> [1] 50

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.

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.

sample_balanced(
  mraster = mraster, # input
  nSamp = 200, # desired sample size
  plot = TRUE
) # plot

#> Simple feature collection with 200 features and 0 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 431110 ymin: 5337730 xmax: 438510 ymax: 5343230
#> 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 (433070 5343230)
#> 2  POINT (436250 5343210)
#> 3  POINT (431810 5343150)
#> 4  POINT (435510 5343150)
#> 5  POINT (431130 5343110)
#> 6  POINT (431350 5343110)
#> 7  POINT (434970 5343110)
#> 8  POINT (434230 5343090)
#> 9  POINT (435430 5343070)
#> 10 POINT (432710 5343050)

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: 431150 ymin: 5337750 xmax: 438390 ymax: 5343230
#> 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 (435790 5343230)
#> 2  POINT (434610 5343210)
#> 3  POINT (433530 5343150)
#> 4  POINT (431690 5343090)
#> 5  POINT (435630 5343090)
#> 6  POINT (432570 5343030)
#> 7  POINT (433910 5343010)
#> 8  POINT (435430 5342930)
#> 9  POINT (432870 5342910)
#> 10 POINT (436070 5342870)

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:

1. Determines the quantile distributions of existing sample units and mraster covariates.

2. Determines quantiles where there is a disparity between sample units and covariates.

3. 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 337 features and 7 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 431170 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 5.11 13.800000 1.21      1    new rule1 POINT (434490 5341390)
#> 2  existing 3.45 48.100002 0.62      1    new rule1 POINT (434970 5341730)
#> 3  existing 7.44 67.200005 1.71      1    new rule1 POINT (438030 5341090)
#> 4  existing 6.96 73.099998 1.57      1    new rule1 POINT (436870 5338030)
#> 5  existing 8.29 95.800003 1.50      1    new rule1 POINT (437230 5337970)
#> 6  existing 9.23 53.799999 2.53      1    new rule1 POINT (436410 5338390)
#> 7  existing 3.06 30.500000 0.54      1    new rule1 POINT (434010 5340850)
#> 8  existing 4.48 17.400000 1.03      1    new rule1 POINT (435590 5342730)
#> 9  existing 8.96  8.400001 2.32      1    new rule1 POINT (433690 5340530)
#> 10 existing 3.65 50.700001 0.65      1    new rule1 POINT (435530 5342490)

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: 431130 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 5.11 13.800000 1.21      1    new rule1 POINT (434490 5341390)
#> 2  existing 3.45 48.100002 0.62      1    new rule1 POINT (434970 5341730)
#> 3  existing 7.44 67.200005 1.71      1    new rule1 POINT (438030 5341090)
#> 4  existing 6.96 73.099998 1.57      1    new rule1 POINT (436870 5338030)
#> 5  existing 8.29 95.800003 1.50      1    new rule1 POINT (437230 5337970)
#> 6  existing 9.23 53.799999 2.53      1    new rule1 POINT (436410 5338390)
#> 7  existing 3.06 30.500000 0.54      1    new rule1 POINT (434010 5340850)
#> 8  existing 4.48 17.400000 1.03      1    new rule1 POINT (435590 5342730)
#> 9  existing 8.96  8.400001 2.32      1    new rule1 POINT (433690 5340530)
#> 10 existing 3.65 50.700001 0.65      1    new rule1 POINT (435530 5342490)

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)

#--- sub sample using ---#
e <- existing %>%
  extract_metrics(mraster = mraster, existing = .)

#--- 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: 431101.5 ymin: 5337703 xmax: 438488.6 ymax: 5343230
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#>      zq90 pzabove2  zsd                 geometry
#> 806 15.20     67.3 3.59 POINT (433227.7 5342252)
#> 60  10.50     89.9 2.54 POINT (436951.6 5338337)
#> 273 21.60     83.3 6.42 POINT (435380.9 5339591)
#> 890 12.90     84.2 3.11 POINT (431810.2 5342518)
#> 10   2.27      2.4 0.31   POINT (438226 5338210)
#> 556  7.98     28.1 2.17 POINT (434655.3 5341138)
#> 327 14.40     90.3 2.96 POINT (437613.2 5342163)
#> 564  8.38     40.9 2.25 POINT (435773.1 5342283)
#> 674  6.55     17.7 1.60 POINT (433384.3 5340981)
#> 874  7.42     52.0 1.82 POINT (432934.8 5343097)

#--- 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: 431101.5 ymin: 5337703 xmax: 438508.8 ymax: 5343230
#> 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
#> 91703 12.60     77.6 3.49 POINT (432423.1 5338565)
#> 91272  2.38      3.7 0.34 POINT (437646.9 5339335)
#> 91276 21.40     96.1 4.17 POINT (438205.8 5339907)
#> 91399 17.30     91.8 3.43 POINT (434409.6 5338023)
#> 91573  2.85     17.4 1.96   POINT (434253 5339294)
#> 91428 15.80     93.4 2.74 POINT (434126.7 5338020)
#> 92085 12.90     84.2 3.11 POINT (431810.2 5342518)
#> 91564 17.30     72.3 5.12 POINT (438308.4 5343161)
#> 91342 17.80     79.7 5.68 POINT (437912.8 5340752)
#> 91949 20.60     62.6 6.04   POINT (432955 5341400)

#--- 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: 431101.5 ymin: 5337710 xmax: 438505.4 ymax: 5343234
#> 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
#> 91325  9.50     82.7 2.26    2.766268 POINT (435118.3 5337890)
#> 91872  3.89     46.5 0.75   18.496631 POINT (433803.5 5341410)
#> 91637 25.70     94.9 4.34  623.097176 POINT (432569.5 5338143)
#> 92018 16.90     84.1 3.79   75.052924 POINT (432525.7 5341820)
#> 91886 17.70     94.1 3.48  602.745800   POINT (431285 5339118)
#> 91291 22.90     92.5 6.91  205.595256 POINT (438342.1 5340333)
#> 91862 16.50     80.8 5.05  214.735502 POINT (432406.2 5339980)
#> 91313 15.30     92.6 3.20  664.483757 POINT (437077.8 5339611)
#> 91589 11.20     83.7 3.06  149.500836 POINT (436628.3 5341727)
#> 91960 17.80     87.0 4.70   63.599623   POINT (434492 5342974)

#--- 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: 431101.5 ymin: 5337710 xmax: 438505.4 ymax: 5343230
#> 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
#> 91207 10.8     58.4 2.87   109.42745 POINT (438082.9 5338350)
#> 91412 11.5     29.9 2.51    89.53668   POINT (435088 5340436)
#> 91317 12.0     88.0 3.11   187.69366 POINT (434406.2 5338306)
#> 91414 19.5     66.7 4.95   282.75431 POINT (436066.1 5341437)
#> 91341 16.4     94.8 3.67   193.23233 POINT (436079.5 5340306)
#> 91220 13.6     85.0 3.52   149.17740 POINT (437377.5 5338200)
#> 91401 11.2     83.7 3.06   149.50084 POINT (436628.3 5341727)
#> 91247 17.5     87.4 4.82   285.19361 POINT (437786.6 5339478)
#> 91404 21.2     87.6 7.21   135.76876   POINT (437327 5342443)
#> 91455 24.0     91.5 5.90   120.42013 POINT (436754.6 5343001)

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: 431101.5 ymin: 5337703 xmax: 438512.2 ymax: 5343230
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#>     zq90 pzabove2  zsd                 geometry
#> 2   2.27      0.5 0.30 POINT (438512.2 5337931)
#> 4  10.30     88.8 2.68 POINT (438229.4 5337927)
#> 8  10.90     76.4 2.84 POINT (437946.5 5337924)
#> 9   5.87     15.4 1.37 POINT (438086.3 5338067)
#> 21  7.33     56.6 1.81 POINT (437241.2 5337774)
#> 23  7.73     77.8 1.73 POINT (437520.6 5338060)
#> 29 17.00     93.1 3.54   POINT (438359 5338919)
#> 30  2.56      6.7 0.40 POINT (438498.7 5339062)
#> 32 10.10     93.9 2.31 POINT (437098.1 5337914)
#> 36  9.41     48.4 2.51   POINT (437657 5338486)

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: 431108.2 ymin: 5337703 xmax: 438512.2 ymax: 5343220
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#>     zq90 pzabove2  zsd                 geometry
#> 2   2.27      0.5 0.30 POINT (438512.2 5337931)
#> 3  14.40     12.8 4.13 POINT (438089.6 5337784)
#> 8  10.90     76.4 2.84 POINT (437946.5 5337924)
#> 9   5.87     15.4 1.37 POINT (438086.3 5338067)
#> 10  2.27      2.4 0.31   POINT (438226 5338210)
#> 11  2.95     12.3 0.50 POINT (438365.7 5338353)
#> 16  8.07      2.4 2.20 POINT (437943.2 5338207)
#> 17 10.80     58.4 2.87 POINT (438082.9 5338350)
#> 18 12.60     62.1 3.65 POINT (438222.6 5338493)
#> 21  7.33     56.6 1.81 POINT (437241.2 5337774)

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: 431138.6 ymin: 5337710 xmax: 438498.7 ymax: 5343237
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#>     zq90 pzabove2  zsd                 geometry
#> 1  11.90     88.2 3.09 POINT (437234.4 5338340)
#> 2  21.10     93.6 4.62 POINT (431288.4 5338835)
#> 3   6.64      9.7 1.75 POINT (434958.4 5339444)
#> 4   4.07     18.8 0.82 POINT (434355.6 5342549)
#> 5   2.51      7.7 0.37 POINT (437742.8 5343155)
#> 6  12.90     93.2 2.82 POINT (431823.7 5341387)
#> 7  16.00     94.3 3.22   POINT (434689 5338310)
#> 8  15.40     91.2 3.30 POINT (432832.1 5339843)
#> 9   7.69     40.4 1.87 POINT (434695.8 5337744)
#> 10 22.80     84.1 4.90 POINT (438461.6 5342173)

sraster <- strat_kmeans(mraster = mraster, nStrata = 4)

e_strata <- extract_strata(sraster = sraster, existing = e)

#--- proportional stratified sampling of existing ---#
sample_existing(existing = e_strata, nSamp = 300, type = "strat", allocation = "prop")
#> Simple feature collection with 300 features and 4 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 431115 ymin: 5337724 xmax: 438505.4 ymax: 5343237
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#>    strata  zq90  pzabove2  zsd                 geometry
#> 1       1 11.50 21.800001 3.69 POINT (437752.9 5342306)
#> 2       1  8.87  8.900001 2.35 POINT (437803.4 5338064)
#> 3       1  3.43  8.000000 0.74 POINT (433530.8 5340559)
#> 4       1  4.47 18.700001 0.89 POINT (437606.4 5342729)
#> 5       1  4.77 11.900001 1.12 POINT (438372.4 5337788)
#> 6       1  1.64  0.000000 0.11 POINT (435517.3 5340017)
#> 7       1  4.16  9.600000 0.90 POINT (438338.7 5340616)
#> 8       1  4.20 15.600000 0.95   POINT (437790 5339195)
#> 9       1  4.67 20.000000 1.06 POINT (434788.3 5341847)
#> 10      1  1.98  0.100000 0.22 POINT (438072.8 5339198)

#--- 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: 431101.5 ymin: 5337703 xmax: 438512.2 ymax: 5343234
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#>    strata  zq90 pzabove2  zsd                 geometry
#> 1       1 11.90      5.0 3.91 POINT (435244.5 5339165)
#> 2       1  5.35     45.1 1.26 POINT (436052.6 5342569)
#> 3       1  1.93      0.0 0.19 POINT (435094.7 5339870)
#> 4       1  7.11     43.7 1.32 POINT (432508.9 5343234)
#> 5       1  2.95     12.3 0.50 POINT (438365.7 5338353)
#> 6       1 14.40     12.8 4.13 POINT (438089.6 5337784)
#> 7       1  6.28     26.8 1.59 POINT (434924.7 5342272)
#> 8       1  9.32      0.5 2.53 POINT (434392.7 5339437)
#> 9       1  5.30     26.1 1.21 POINT (435064.4 5342416)
#> 10      1  3.30     28.3 0.86 POINT (434931.4 5341707)

#--- 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))

#--- check proportions match weights ---#
table(s$strata) / 100
#> 
#>   1   2   3   4 
#> 0.2 0.6 0.1 0.1

#--- 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 99 features and 4 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 431391.1 ymin: 5337717 xmax: 438468.4 ymax: 5343128
#> Projected CRS: UTM Zone 17, Northern Hemisphere
#> First 10 features:
#>    strata  zq90 pzabove2  zsd                 geometry
#> 1       1  7.29     24.1 1.99 POINT (437037.4 5343005)
#> 2       1  3.89     46.5 0.75 POINT (433803.5 5341410)
#> 3       1  4.07     18.8 0.82 POINT (434355.6 5342549)
#> 4       1  1.64      0.0 0.11 POINT (435517.3 5340017)
#> 5       1  4.20     15.6 0.95   POINT (437790 5339195)
#> 6       1 11.50     21.8 3.69 POINT (437752.9 5342306)
#> 7       1  8.83     28.2 2.29 POINT (437796.7 5338629)
#> 8       1  6.55     17.7 1.60 POINT (433384.3 5340981)
#> 9       1  2.22      0.3 0.29 POINT (434256.4 5339012)
#> 10      1  6.80     47.0 1.60 POINT (434648.6 5341703)