County-to-Watershed Artifact Spatial Apportioning

Reading in the Data

library(terra)

terra 1.9.34

library(stringr)

# directory structure
main_dir <- "S:/ursa/campbell/Other/tune/data/for_analysis"

# inputs
counties <- file.path(main_dir, "counties.shp") |> vect()
huc04 <- file.path(main_dir, "huc04.shp") |> vect()
huc06 <- file.path(main_dir, "huc06.shp") |> vect()
huc08 <- file.path(main_dir, "huc08.shp") |> vect()
streamdist <- file.path(main_dir, "dist_to_steams_250cfs_clip.tif") |> rast()
artifacts <- file.path(main_dir, "artifacts.csv") |> read.csv()

# join artifacts data to counties
artifacts$FIPS <- str_pad(artifacts$FPS, 5, "left", "0")
counties <- merge(counties, artifacts, by = "FIPS")

# select the watersheds within the three main basins
keepers <- c("0601", "0602", "0603", "0604")
huc04 <- huc04[huc04$huc4 %in% keepers,]
huc06 <- huc06[str_sub(huc06$huc6, 1, 4) %in% keepers,]
huc08 <- huc08[str_sub(huc08$huc8, 1, 4) %in% keepers,]

# map artifacts by county and overlay watersheds
plot(counties, "Total.Sample", breaks = 10, bg = "white")
lines(huc08, col = "yellow", lwd = 0.5)
lines(huc06, col = "orange", lwd = 2)
lines(huc04, col = "red", lwd = 4)

Approach #1: Random Apportioning

# create sample points
counties_gt0 <- counties[counties$Total.Sample > 0,]
pts_rand <- spatSample(counties_gt0, counties_gt0$Total.Sample)
plot(pts_rand, cex = 0.1)
lines(counties, col = "red")

Approach #2: Stream Distance-Weighted Apportioning

# map distance to streams
plot(streamdist)
lines(counties, col = "red")

# get inverse distance weighted probability
decay_fact <- 5000
prob_rast <- exp(-streamdist/decay_fact)
plot(prob_rast)
lines(counties, col = "red")

# define function to generate weighted random points for a single county
sample_county_sites <- function(i, counties, prob_raster) {
  
  # select county
  county <- counties[i, ]
  
  # get sample size from county attribute
  n_points <- county$Total.Sample
  
  # return NULL or empty vector if count is zero or missing
  if (is.na(n_points) || n_points <= 0) return(NULL)
  
  # crop and mask the probability raster to the individual county boundary
  c_prob <- crop(prob_raster, county, snap = "out")
  c_prob <- mask(c_prob, county)
  
  # sample cell coordinates weighted by probability
  pts <- spatSample(
    x = c_prob,
    size = n_points,
    method = "weights",
    replace = TRUE,
    as.points = TRUE,
    values = FALSE,
    na.rm = TRUE
  )
  
  # append FIPS to points for tracking
  if (!is.null(pts) && nrow(pts) > 0) {
    pts$county_id <- county$FIPS
  }
  
  # return
  return(pts)
  
}

# run the function over all counties
sampled_list <- lapply(seq_len(nrow(counties_gt0)), function(i) {
  sample_county_sites(i, counties_gt0, prob_rast)
})

# combine all county points into a single SpatVector
sampled_list <- sampled_list[!sapply(sampled_list, is.null)]
pts_weight <- do.call(rbind, sampled_list)

# plot them out
plot(pts_weight, cex = 0.1)
lines(counties, col = "red")

Summarize Counts by Watershed

# count points by watershed -- random
huc04_rand <- extract(huc04, pts_rand)
huc04_rand <- table(huc04_rand$huc4) |> as.data.frame()
colnames(huc04_rand) <- c("huc4", "n_rand")

huc06_rand <- extract(huc06, pts_rand)
huc06_rand <- table(huc06_rand$huc6) |> as.data.frame()
colnames(huc06_rand) <- c("huc6", "n_rand")

huc08_rand <- extract(huc08, pts_rand)
huc08_rand <- table(huc08_rand$huc8) |> as.data.frame()
colnames(huc08_rand) <- c("huc8", "n_rand")

# count points by watershed -- weighted
huc04_weight <- extract(huc04, pts_weight)
huc04_weight <- table(huc04_weight$huc4) |> as.data.frame()
colnames(huc04_weight) <- c("huc4", "n_weight")

huc06_weight <- extract(huc06, pts_weight)
huc06_weight <- table(huc06_weight$huc6) |> as.data.frame()
colnames(huc06_weight) <- c("huc6", "n_weight")

huc08_weight <- extract(huc08, pts_weight)
huc08_weight <- table(huc08_weight$huc8) |> as.data.frame()
colnames(huc08_weight) <- c("huc8", "n_weight")

# join them back to the watersheds
huc04 <- merge(huc04, huc04_rand) |>
  merge(huc04_weight)
huc06 <- merge(huc06, huc06_rand) |>
  merge(huc06_weight)
huc08 <- merge(huc08, huc08_rand) |>
  merge(huc08_weight)

# plot side-by-side comparisons
par(mfrow = c(3,2))
plot(huc04, "n_rand", main = "HUC4 Random")
plot(huc04, "n_weight", main = "HUC4 Weighted")
plot(huc06, "n_rand", main = "HUC6 Random")
plot(huc06, "n_weight", main = "HUC6 Weighted")
plot(huc08, "n_rand", main = "HUC8 Random")
plot(huc08, "n_weight", main = "HUC8 Weighted")

Get Differences

# calculate differences
huc04$n_diff <- huc04$n_rand - huc04$n_weight
huc06$n_diff <- huc06$n_rand - huc06$n_weight
huc08$n_diff <- huc08$n_rand - huc08$n_weight

# plot differences
par(mfrow = c(3,1))
plot(huc04, "n_diff")
plot(huc06, "n_diff")
plot(huc08, "n_diff")