## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse  = TRUE,
  comment   = "#>",
  eval      = FALSE,
  fig.width = 9,
  fig.height = 6
)

## ----install_optional---------------------------------------------------------
# install.packages(c("terra", "exactextractr", "ggrepel"))

## ----run_ms-------------------------------------------------------------------
# library(manureshed)
# 
# ms <- run_builtin_analysis(
#   scale        = "county",
#   year         = 2016,
#   nutrients    = c("nitrogen", "phosphorus"),
#   include_wwtp = TRUE,
#   verbose      = TRUE
# )

## ----download_cafo------------------------------------------------------------
# # First call downloads (~5 MB RDS); subsequent calls use the cache
# cafo_path <- download_cafo_data()

## ----cache_option-------------------------------------------------------------
# options(manureshed.cache_dir = "~/manureshed_cache")

## ----score_basic--------------------------------------------------------------
# hub <- score_hub_sites(
#   ms_output       = ms,
#   catchment_miles = 50      # user-controlled; default 50
# )
# 
# print(hub)

## ----score_subset-------------------------------------------------------------
# hub9 <- score_hub_sites(
#   ms_output       = ms,
#   scores          = "Score9_S3_NP",
#   catchment_miles = 50
# )

## ----score_cdl----------------------------------------------------------------
# hub <- score_hub_sites(
#   ms_output        = ms,
#   cdl_path         = "path/to/cdl_2016_conus.tif",
#   catchment_miles  = 50,
#   n_threads        = 8        # increase on HPC
# )

## ----inspect------------------------------------------------------------------
# # Top-10 counties for the flagship score
# hub$top_n[["Score9_S3_NP"]]
# 
# # Which counties appear across the most scores?
# head(hub$robustness, 10)
# 
# # Full master table with all metrics and scores
# names(hub$master)
# nrow(hub$master)   # one row per county

## ----map----------------------------------------------------------------------
# # Static map -- save to file
# map_hub_sites(hub,
#               score     = "Score9_S3_NP",
#               save_path = "hub_flagship.png",
#               width     = 17,
#               height    = 8,
#               dpi       = 300)
# 
# # Return the ggplot object without saving
# p <- map_hub_sites(hub, score = "Score3_S1_NP")

## ----export-------------------------------------------------------------------
# export_hub_results(hub, output_dir = "hub_outputs")

## ----multiyear----------------------------------------------------------------
# years    <- c(2007, 2012, 2016)
# 
# hub_list <- lapply(years, function(yr) {
#   ms_yr <- run_builtin_analysis(
#     scale        = "county",
#     year         = yr,
#     nutrients    = c("nitrogen", "phosphorus"),
#     include_wwtp = TRUE,
#     verbose      = FALSE
#   )
#   score_hub_sites(ms_yr, scores = "Score9_S3_NP", verbose = FALSE)
# })
# names(hub_list) <- paste0("yr_", years)
# 
# # Counties in the flagship top-10 across all three years
# top_by_year  <- lapply(hub_list, function(h) h$top_n[["Score9_S3_NP"]]$FIPS)
# stable_sites <- Reduce(intersect, top_by_year)
# 
# cat("Counties stable across all years:", length(stable_sites), "\n")
# print(stable_sites)

## ----catchment----------------------------------------------------------------
# hub_50  <- score_hub_sites(ms, catchment_miles = 50,
#                             scores = "Score9_S3_NP", verbose = FALSE)
# hub_100 <- score_hub_sites(ms, catchment_miles = 100,
#                             scores = "Score9_S3_NP", verbose = FALSE)
# 
# t50  <- hub_50$top_n[["Score9_S3_NP"]]$FIPS
# t100 <- hub_100$top_n[["Score9_S3_NP"]]$FIPS
# 
# cat("Overlap 50 vs 100 miles:", sum(t50 %in% t100), "/ 10\n")