diff --git a/.Rbuildignore b/.Rbuildignore
index 43df59c..931b37f 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -7,3 +7,4 @@
 ^docs$
 ^pkgdown$
 add_data_temp.R
+data-raw
diff --git a/.github/workflows/check-standard.yaml b/.github/workflows/check-standard.yaml
index a3ac618..9c45ee1 100644
--- a/.github/workflows/check-standard.yaml
+++ b/.github/workflows/check-standard.yaml
@@ -4,7 +4,7 @@ on:
   push:
     branches: [main, master]
   pull_request:
-    branches: [main, master]
+    branches: [main, master, dev]
 
 name: R-CMD-check
 
diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml
index 2c5bb50..23171ad 100644
--- a/.github/workflows/test-coverage.yaml
+++ b/.github/workflows/test-coverage.yaml
@@ -4,7 +4,7 @@ on:
   push:
     branches: [main, master]
   pull_request:
-    branches: [main, master]
+    branches: [main, master, dev]
 
 name: test-coverage
 
diff --git a/DESCRIPTION b/DESCRIPTION
index a88819a..2d67968 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: site
 Type: Package
 Title: Site specific malariasimulation modelling
-Version: 0.2.2
+Version: 1.0.4
 Authors@R: c(
     person("Pete", "Winskill", email = "p.winskill@imperial.ac.uk", role = c("aut", "cre"))
     )
@@ -11,14 +11,20 @@ License: MIT + file LICENSE
 Encoding: UTF-8
 LazyData: true
 Remotes:
-  mrc-ide/malariasimulation
-Suggests: 
-    testthat (>= 3.0.0)
+    mrc-ide/malariasimulation,
+    mrc-ide/orderly2
 Config/testthat/edition: 3
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.2
 Depends: 
     R (>= 2.10)
 Imports: 
+    orderly2,
     dplyr,
-    malariasimulation (>= 1.4.1),
-    tidyr
+    malariasimulation,
+    rappdirs,
+    rlang,
+    tidyr,
+    withr
+Suggests:
+    testthat (>= 3.0.0),
+    fs
diff --git a/NAMESPACE b/NAMESPACE
index 485614d..a453865 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,4 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
-export(single_site)
+export(fetch_files)
+export(fetch_site)
 export(site_parameters)
+export(subset_site)
diff --git a/R/burnin.R b/R/burnin.R
index 38fe53e..66d3376 100644
--- a/R/burnin.R
+++ b/R/burnin.R
@@ -1,6 +1,6 @@
 #' Expand interventions retrospectively for burn in period
 #'
-#' Assumes interventions for burn in period are = to those in year 1.
+#' Assumes interventions for burn in period are = 0 coverage.
 #'
 #' @param interventions Site intervention inputs
 #' @param burnin Burn in period (years)
@@ -10,6 +10,19 @@ burnin_interventions <- function(interventions, burnin){
   start_year <- min(interventions$year) - burnin
   interventions <- interventions |>
     tidyr::complete(year = start_year:(max(interventions$year))) |>
+    tidyr::replace_na(
+      replace = list(
+        itn_use = 0,
+        itn_input_dist = 0,
+        tx_cov = 0,
+        irs_cov = 0,
+        smc_cov = 0,
+        rtss_cov = 0,
+        r21_cov = 0,
+        pmc_cov = 0,
+        lsm_cov = 0
+      )
+    ) |>
     tidyr::fill(dplyr::everything(), .direction = "up")
   return(interventions)
 }
diff --git a/R/data.R b/R/data.R
index 55fbd33..089e31b 100644
--- a/R/data.R
+++ b/R/data.R
@@ -5,14 +5,20 @@
 #'
 #' @format A list with 10 variables:
 #' \describe{
-#'   \item{country}{The iso3c country code}
-#'   \item{level}{The level of subnational disaggregation}
+#'   \item{country}{The country name}
+#'   \item{version}{Site file version}
+#'   \item{admin_level}{The levels of spatial disaggregation}
 #'   \item{sites}{Unique sites}
-#'   \item{epi}{Epidemiological site data}
+#'   \item{shape}{Sptial boundaries (set to NA in testing shapefile)}
+#'   \item{cases_deaths}{Epidemiological site data}
+#'   \item{prevalence}{Prevalence data}
 #'   \item{interventions}{Intervention coverage and specification}
-#'   \item{population}{Population and population at risk projections}
+#'   \item{population}{Population, population at risk mand age-disaggregated population projections}
 #'   \item{demography}{Demographic projections}
-#'   \item{vectors}{Vector proportions}
-#'   \item{seasonality}{Seasonal profile parameters}
+#'   \item{vectors}{Vector specieis and pyrethroid resistance}
+#'   \item{seasonality}{Seasonal profile parameters, monthly rainfall and fourier predictions}
+#'   \item{blood_disorders}{Blood disorder data}
+#'   \item{accessibility}{Accessibility data}
+#'   \item{eir}{Calibrated eir}
 #' }
 "example_site"
diff --git a/R/demography.R b/R/demography.R
index cc5a577..c4a310d 100644
--- a/R/demography.R
+++ b/R/demography.R
@@ -8,11 +8,8 @@ add_demography <- function(p, demography){
 
   # Age group upper
   ages <- round(unique(demography$age_upper) * 365)
-  # Single demography currently
   timesteps <- 365 * (unique(demography$year) - p$baseline_year)
-  # Take demography is first year as static demography
-  deathrates <- demography$mortality_rate / 365
-  # Create matrix of death rates
+  deathrates <- demography$adjusted_mortality_rates / 365
   deathrates_matrix <- matrix(deathrates, nrow = length(timesteps), byrow = TRUE)
   # Add parameters
   p <- malariasimulation::set_demography(
diff --git a/R/fetch.R b/R/fetch.R
new file mode 100644
index 0000000..7e33fcd
--- /dev/null
+++ b/R/fetch.R
@@ -0,0 +1,163 @@
+LOCATION_NAME <- "malariaverse-sitefiles"
+
+location_configuration <- function() {
+  token <- Sys.getenv("GITHUB_TOKEN")
+  if (token == "") {
+    token <- NULL
+  }
+
+  getOption("site.orderly_location", list(
+    type = "packit",
+    args = list(
+      url = "https://packit.dide.ic.ac.uk/malariaverse-sitefiles",
+      token = token)
+  ))
+}
+
+#' Add or update an Orderly location.
+#'
+#' If a location with the given name already exists and its configuration does
+#' not match the given parameters, it is removed first before being added with
+#' the new parameters. If it exists and has the same parameters already nothing
+#' happens.
+#'
+#' This functionality could probably be moved to the orderly2 package.
+#' @noRd
+location_add_or_update <- function(name, type, args, root) {
+  locations <- orderly2::orderly_location_list(root = root, verbose = TRUE)
+  locations <- locations[locations$name == name,]
+
+  if (nrow(locations) == 0) {
+    orderly2::orderly_location_add(name, type, args, root = root)
+  } else if (locations[[1, "type"]] != type ||
+             !identical(locations[[1, "args"]], args)) {
+    orderly2::orderly_location_remove(name, root = root)
+    orderly2::orderly_location_add(name, type, args, root = root)
+  }
+}
+
+
+#' Configure the orderly root used to fetch sitefiles.
+#' 
+#' This creates a folder in the user's home directory used to download and cache
+#' site files. The location of the cache folder is determined by
+#' [rappdirs::user_cache_dir()] and depends on the OS.
+#'
+#' A remote location from which the sitefiles will be fetched is configured on
+#' the root. By default this is the malariaverse Packit instance hosted at
+#' `https://packit.dide.ic.ac.uk/malariaverse-sitefiles`. This can be customized
+#' by setting the `site.orderly_location` option.
+#'
+#' Users shouldn't need to call this function, as it is called implicitly by
+#' [fetch_files] already.
+#'
+#' @return the path to the orderly root.
+#' @noRd
+configure_orderly <- function() {
+  root <- file.path(rappdirs::user_cache_dir("malariaverse-sitefiles"), "store")
+
+  orderly2::orderly_init(root, use_file_store = TRUE)
+
+  cfg <- location_configuration()
+  location_add_or_update(LOCATION_NAME, type = cfg$type, args = cfg$args,
+                         root = root)
+
+  root
+}
+
+
+#' Get files from the malariaverse sitefile server.
+#'
+#' @param name The name of the orderly report.
+#' @param parameters A named list of parameters to use when searching for the
+#'   orderly packet. If a query expression `expr` is specified, these parameters
+#'   are substituted into the query using the this: prefix. If no expression is
+#'   specified, the latest packet matching these parameters exactly is selected.
+#' @param dest A directory into which the files should be copied.
+#' @param files An optionally-named character vector of files to copy from the
+#'   packet and into the destination directory. If the vector is named, these
+#'   names are used as the destination file path.
+#' @param expr The query expression to filter packets. This may be an arbitrary
+#'   orderly query, including a literal packet ID. If absent or NULL, the
+#'   specified list of parameters is used and matched exactly.
+#' @return the id of the orderly packet the files were copied from.
+#' @export
+fetch_files <- function(name, parameters, dest, files, expr = NULL) {
+  root <- configure_orderly()
+
+  if (is.null(expr)) {
+    filter <- paste(sprintf("parameter:%1$s == this:%1$s", names(parameters)),
+                    collapse = " && ")
+    expr <- sprintf("latest(%s)", filter)
+  }
+
+  options <- orderly2::orderly_search_options(
+    location = LOCATION_NAME,
+    allow_remote = TRUE,
+    pull_metadata = TRUE)
+
+  plan <- orderly2::orderly_copy_files(
+    name = name,
+    expr = expr,
+    parameters = parameters,
+    dest = dest,
+    files = files,
+    options = options,
+    root = root)
+
+  plan$id
+}
+
+#' Fetch a site file for a country from the malariaverse sitefile server.
+#'
+#' The site file is identified by its country code, and optionally the
+#' admin_level, urban/rural setting and version of the site files. The latest
+#' packet from the server matching these parameters is used.
+#'
+#' Alternatively, a packet ID can be specified in order to pick an exact file
+#' set.
+#'
+#' @param iso3c the ISO country code, a scalar character.
+#' @param version the dataset version, a scalar character.
+#' @param admin_level a scalar number.
+#' @param urban_rural a scalar logical.
+#' @param id a packet ID used to select an exact packet.
+#' @return The contents of the site file.
+#' @examples
+#' \dontrun{
+#' fetch_site("NGA")
+#' fetch_site("NGA", admin_level = 1)
+#' fetch_site(id = "20240801-062621-6f95851a")
+#' }
+#' @export
+fetch_site <- function(iso3c = NULL, version = NULL,
+                       admin_level = NULL, urban_rural = NULL,
+                       id = NULL)
+{
+  dest <- withr::local_tempdir()
+  if (!xor(is.null(iso3c), is.null(id))) {
+    rlang::abort("Exactly one of `iso3c` and `id` must be supplied")
+  }
+
+  if (!is.null(iso3c)) {
+    parameters <- list(
+      iso3c = iso3c,
+      version = version,
+      admin_level = admin_level,
+      urban_rural = urban_rural)
+    parameters <- parameters[!sapply(parameters, is.null)]
+    expr <- NULL
+  } else {
+    parameters <- list()
+    expr <- id
+  }
+
+  fetch_files(name = "calibration_diagnostics",
+              files = "calibrated_scaled_site.rds",
+              expr = expr,
+              parameters = parameters,
+              dest = dest)
+
+
+  readRDS(file.path(dest, "calibrated_scaled_site.rds"))
+}
diff --git a/R/interventions.R b/R/interventions.R
index b55bf4d..4a3b926 100644
--- a/R/interventions.R
+++ b/R/interventions.R
@@ -2,12 +2,13 @@
 #'
 #' @param p parameter list
 #' @param interventions site intervention inputs
-#' @param species Can be falciparum: "pf" or vivax: "pv", for vivax SMC, RTSS
-#'  and PMC are not implemented
+#' @param species Can be falciparum: "pf" or vivax: "pv", for vivax SMC, RTSS,
+#'  R21 and PMC are not implemented
 #'
 #' @return modified parameter list
 add_interventions <- function(p, interventions, species){
 
+  pf <- species ==  "pf"
   # Drug types
   p <- add_drugs(p)
   # Treatment
@@ -18,31 +19,65 @@ add_interventions <- function(p, interventions, species){
   if(sum(interventions$itn_input_dist, na.rm = TRUE) > 0){
     p <- add_itns(
       p = p,
-      interventions = interventions)
+      interventions = interventions
+    )
   }
   # IRS
   if(sum(interventions$irs_cov, na.rm = TRUE) > 0){
     p <- add_irs(
       p = p,
-      interventions = interventions)
+      interventions = interventions
+    )
   }
   # SMC
-  if(sum(interventions$smc_cov, na.rm = TRUE) > 0 &
-     species ==  "pf"){
-    p <- add_smc(p = p,
-                 interventions = interventions)
+  if(sum(interventions$smc_cov, na.rm = TRUE) > 0 & pf){
+    p <- add_smc(
+      p = p,
+      interventions = interventions
+    )
+  }
+
+  if(
+    sum(interventions$rtss_cov, na.rm = TRUE) > 0 &
+    sum(interventions$r21_cov, na.rm = TRUE) > 0
+  ) {
+    warning("Cannot currently model two vaccine types,
+            defaulting to R21 vaccine type implemented at
+            the maximum yearly coverage inputs across rtss
+            and r21")
+    interventions$r21_cov <- pmax(interventions$r21_cov, interventions$rtss_cov)
+    interventions$rtss_cov <- 0
   }
   # RTSS
-  if(sum(interventions$rtss_cov, na.rm = TRUE) > 0 &
-     species ==  "pf"){
-    p <- add_rtss(p = p,
-                  interventions = interventions)
+  if(sum(interventions$rtss_cov, na.rm = TRUE) > 0 & pf){
+    p <- add_rtss(
+      p = p,
+      interventions = interventions
+    )
+  }
+  # R21
+  if(sum(interventions$r21_cov, na.rm = TRUE) > 0 & pf){
+    p <- add_r21(
+      p = p,
+      interventions = interventions
+    )
   }
   # PMC
-  if(sum(interventions$pmc_cov, na.rm = TRUE) > 0 &
-     species ==  "pf"){
-    p <- add_pmc(p = p,
-                 interventions = interventions)
+  if(sum(interventions$pmc_cov, na.rm = TRUE) > 0 & pf){
+    p <- add_pmc(
+      p = p,
+      interventions = interventions
+    )
+  }
+  # Interventions that modify the carrying capacity
+  # The combined carrying capacity scaling must be estimated for all
+  # interventions that modify it, before updating. Currently, only LSM is
+  # implemented here. This could included An. stephensi in the future.
+  if(sum(interventions$lsm_cov, na.rm = TRUE) > 0){
+    p <- adjust_carrying_capacity(
+      p = p,
+      interventions = interventions
+    )
   }
 
   return(p)
@@ -68,13 +103,14 @@ add_drugs <- function(p){
 
   p <- malariasimulation::set_drugs(
     parameters = p,
-    list(tetst = SP_params,
-         malariasimulation::AL_params,
-         malariasimulation::SP_AQ_params,
-         SP_full_efficacy,
-         AL_full_efficacy))
-
-
+    list(
+      SP_params,
+      malariasimulation::AL_params,
+      malariasimulation::SP_AQ_params,
+      SP_full_efficacy,
+      AL_full_efficacy
+    )
+  )
   return(p)
 }
 
@@ -88,15 +124,19 @@ add_treatment <- function(p, interventions){
   timesteps <- 1 + (interventions$year - p$baseline_year) * 365
 
   # Non ACT (SP)
-  p <- malariasimulation::set_clinical_treatment(parameters = p,
-                                                 drug = 4,
-                                                 timesteps = timesteps,
-                                                 coverages = interventions$tx_cov * (1 - interventions$prop_act))
+  p <- malariasimulation::set_clinical_treatment(
+    parameters = p,
+    drug = 4,
+    timesteps = timesteps,
+    coverages = interventions$tx_cov * (1 - interventions$prop_act)
+  )
   # ACT (AL)
-  p <- malariasimulation::set_clinical_treatment(parameters = p,
-                                                 drug = 5,
-                                                 timesteps = timesteps,
-                                                 coverages = interventions$tx_cov * interventions$prop_act)
+  p <- malariasimulation::set_clinical_treatment(
+    parameters = p,
+    drug = 5,
+    timesteps = timesteps,
+    coverages = interventions$tx_cov * interventions$prop_act
+  )
 
 
   return(p)
@@ -108,10 +148,17 @@ add_treatment <- function(p, interventions){
 #'
 #' @return modified parameter list
 add_itns <- function(p, interventions){
-  # Assuming net distribution happens on January 1st
-  timesteps <- 1 + (interventions$year - p$baseline_year) * 365
+
+  # If not specified, assume distribution happens January 1st
+  if(!"itn_distribution_day" %in% colnames(interventions)){
+    interventions$itn_distribution_day <- 1
+  }
+  timesteps <- interventions$itn_distribution_day + (interventions$year - p$baseline_year) * 365
   # Net retention half life does not vary over time (Should match what is used when fitting input dist)
-  retention <- 365 * 5
+  retention <- unique(interventions$mean_retention)
+  if(length(retention) > 1){
+    stop("Time-varying net rentetion is not currently supported")
+  }
   # Net input coverage
   coverages <- interventions$itn_input_dist
   coverages[is.na(coverages)] <- 0
@@ -145,17 +192,28 @@ add_irs <- function(p, interventions){
   month <- 365 / 12
   peak <- malariasimulation::peak_season_offset(p)
   year_start_times <-  1 + (interventions$year - p$baseline_year) * 365
-  peak_season_times <- peak + year_start_times
-  # Assume IRS occurs 3 months before seasonal peak
-  irs_spray_times <- round(peak_season_times - 3 * month)
-  coverages <- interventions$irs_cov
-  ls_theta <- interventions$ls_theta
-  ls_gamma <- interventions$ls_gamma
-  ks_theta <- interventions$ks_theta
-  ks_gamma <- interventions$ks_gamma
-  ms_theta <- interventions$ms_theta
-  ms_gamma <- interventions$ms_gamma
 
+  if(any(interventions$irs_spray_rounds > 2)){
+    stop("Maximum of 2 IRS spray rounds per year supported")
+  }
+
+  peak_season_times <- peak + year_start_times
+  peak_season_times <- rep(peak_season_times, interventions$irs_spray_rounds)
+
+  # Assume multiple spray rounds are 6 months apart (First is 3 months prior to peak)
+  offset <- sapply(interventions$irs_spray_rounds, function(x){
+    c(3, -3)[1:x]
+  }) |>
+    unlist()
+  irs_spray_times <- round(peak_season_times - offset * month)
+
+  coverages <- rep(interventions$irs_cov, interventions$irs_spray_rounds)
+  ls_theta <- rep(interventions$ls_theta, interventions$irs_spray_rounds)
+  ls_gamma <- rep(interventions$ls_gamma, interventions$irs_spray_rounds)
+  ks_theta <- rep(interventions$ks_theta, interventions$irs_spray_rounds)
+  ks_gamma <- rep(interventions$ks_gamma, interventions$irs_spray_rounds)
+  ms_theta <- rep(interventions$ms_theta, interventions$irs_spray_rounds)
+  ms_gamma <- rep(interventions$ms_gamma, interventions$irs_spray_rounds)
 
   index <- irs_spray_times <= 0
   if(sum(index) > 0){
@@ -199,7 +257,6 @@ add_smc <- function(p, interventions){
   }
 
   peak <- malariasimulation::peak_season_offset(p)
-  # Note: min age and max age are not currently time-varying
   rounds <- interventions$smc_n_rounds
   year_start_times <-  1 + (interventions$year - p$baseline_year) * 365
   peak_season_times <- peak + year_start_times
@@ -228,7 +285,8 @@ add_smc <- function(p, interventions){
     timesteps = timesteps,
     coverages = coverages,
     min_age = min_age,
-    max_age = max_age)
+    max_age = max_age
+  )
 
   return(p)
 }
@@ -242,15 +300,73 @@ add_rtss <- function(p, interventions){
   month <- 365 / 12
   timesteps <- 1 + (interventions$year - p$baseline_year) * 365
 
-  p <- malariasimulation::set_rtss_epi(
+  if("n_doses" %in% colnames(interventions)){
+    booster_cov <- rep(0, length(timesteps))
+    booster_cov[interventions$n_doses == 4] <- 0.8
+  } else {
+    booster_cov <- rep(0.8, length(timesteps))
+  }
+
+  p <- malariasimulation::set_pev_epi(
     parameters = p,
-    timesteps = timesteps,
+    profile = malariasimulation::rtss_profile,
     coverages = interventions$rtss_cov,
+    timesteps = timesteps,
+    age = round(6 * month),
+    min_wait = 0,
+    booster_spacing = 12 * month, # The booster is administered 12 months following the third dose.
+    booster_coverage = matrix(booster_cov),
+    booster_profile = list(malariasimulation::rtss_booster_profile) # We will model implementation of the RTSS booster.
+  )
+
+  return(p)
+}
+
+#' Add RTS,S
+#'
+#' @inheritParams add_interventions
+#'
+#' @return modified parameter list
+add_r21 <- function(p, interventions){
+  month <- 365 / 12
+  timesteps <- 1 + (interventions$year - p$baseline_year) * 365
+  if("n_doses" %in% colnames(interventions)){
+    booster_cov <- rep(0, length(timesteps))
+    booster_cov[interventions$n_doses == 4] <- 0.8
+  } else {
+    booster_cov <- rep(0.8, length(timesteps))
+  }
+
+  r21_profile <- malariasimulation::create_pev_profile(
+    vmax = 0.87,
+    alpha = 0.91,
+    beta = 471,
+    cs = c(9.3199677, 0.8387902),
+    rho = c(0.8071676, 0.6010363),
+    ds = c(3.7996007, 0.1618982),
+    dl =c(6.2820200, 0.4549185)
+  )
+  r21_booster_profile <- malariasimulation::create_pev_profile(
+    vmax = 0.87,
+    alpha = 0.91,
+    beta = 471,
+    cs = c(9.2372858, 0.7188541),
+    rho = c(0.07140337, 0.54175154),
+    ds = c(3.7996007, 0.1618982),
+    dl =c(6.2820200, 0.4549185)
+  )
+
+  p <- malariasimulation::set_pev_epi(
+    parameters = p,
+    profile = r21_profile,
+    coverages = interventions$r21_cov,
+    timesteps = timesteps,
+    # TODO: Check R21 timings/ages
     age = round(6 * month),
     min_wait = 0,
-    boosters = round(18 * month),
-    booster_coverage = 0.8,
-    seasonal_boosters = FALSE
+    booster_spacing = 12 * month, # The booster is administered 12 months following the third dose.
+    booster_coverage = matrix(booster_cov),
+    booster_profile = list(r21_booster_profile) # We will model implementation of the RTSS booster.
   )
 
   return(p)
@@ -280,3 +396,26 @@ add_pmc <- function(p, interventions){
   return(p)
 }
 
+#' Adjust carrying capacity
+#'
+#' @inheritParams add_interventions
+#'
+#' @return modified parameter list
+adjust_carrying_capacity <- function(p, interventions){
+  lsm_impact <- rep(1 - interventions$lsm_cov, each = length(p$species))
+  carrying_capacity_scaler <- matrix(
+    data = lsm_impact,
+    ncol = length(p$species),
+    byrow = TRUE
+  )
+  month <- 365 / 12
+  timesteps <- 1 + (interventions$year - p$baseline_year) * 365
+
+  p <- malariasimulation::set_carrying_capacity(
+    parameters = p,
+    carrying_capacity = carrying_capacity_scaler,
+    timesteps = timesteps
+  )
+
+  return(p)
+}
diff --git a/R/outputs.R b/R/outputs.R
index 17fb845..9a40f19 100644
--- a/R/outputs.R
+++ b/R/outputs.R
@@ -5,16 +5,15 @@
 #'
 #' @return modified parameter list
 set_age_outputs <- function(p, min_ages){
-  max_ages <- c(min_ages[-1], 100 * 365) - 1
-
-  p$age_group_rendering_min_ages = min_ages
-  p$age_group_rendering_max_ages = max_ages
-
-  p$clinical_incidence_rendering_min_ages = min_ages
-  p$clinical_incidence_rendering_max_ages = max_ages
-
-  p$severe_incidence_rendering_min_ages = min_ages
-  p$severe_incidence_rendering_max_ages = max_ages
+  # Max age upper bound is the maximimum of 100 years, or maximum of min_ages + 1 year
+  max_upper_bound <- max(100 * 365, min_ages[length(min_ages)] + 365)
+  ages <- c(min_ages, max_upper_bound)
+  p <- p |>
+    malariasimulation::set_epi_outputs(
+      age_group = ages,
+      clinical_incidence = ages,
+      severe_incidence = ages
+    )
 
   # PfPr: 2-10, PvPr: 1-100
   p$prevalence_rendering_min_ages = c(1, 2) * 365
diff --git a/R/single_site.R b/R/single_site.R
index 8ca543c..a468ee7 100644
--- a/R/single_site.R
+++ b/R/single_site.R
@@ -1,23 +1,59 @@
-#' Extract a single site-input from a country site file
+#' Extract a subset from a country site file
 #'
-#' @param site_file Country site file
-#' @param index Index row from site_file$sites
+#' @param site Country site file
+#' @param site_filter Data.frame to filter site file elements by. Filtering will
+#' be conducted on any matched columns for each element.
 #'
-#' @return Single site
+#' @return Filtered site file
 #' @export
-single_site <- function(site_file, index){
-  if(index < 1 | index > nrow(site_file$sites)){
-    stop("Mis-specified site index")
-  }
-  index_site <- site_file$sites[index, ]
+subset_site <- function(site, site_filter){
 
-  to_mod <- c("sites", "interventions", "pyrethroid_resistance", "population",
-              "vectors", "seasonality", "prevalence", "eir")
+  sub_site <- list()
+  sub_site$country <- site$country
+  sub_site$version <- site$version
+  sub_site$admin_level <- site$admin_level
+  sub_site$sites <- site_filter[ , sub_site$admin_level]
+  sub_site$shape <- lapply(
+    site$shape,
+    match,
+    y = site_filter
+  )
+  names(sub_site$shape) <- names(site$shape)
+  sub_site$cases_deaths <- match(site$cases_deaths, site_filter)
+  sub_site$prevalence <- match(site$prevalence, site_filter)
+  sub_site$interventions <- match(site$interventions, sub_site$site)
+  sub_site$population <- list(
+    population_total = match(site$population$population_total, sub_site$site),
+    population_by_age = match(site$population$population_by_age, sub_site$site)
+  )
+  sub_site$demography <- match(site$demography, sub_site$site)
+  sub_site$vectors <-  list(
+    vector_species = match(site$vectors$vector_species, sub_site$site),
+    pyrethroid_resistance = match(site$vectors$pyrethroid_resistance, sub_site$site)
+  )
+  sub_site$seasonality <-  list(
+    seasonality_parameters = match(site$seasonality$seasonality_parameters, sub_site$site),
+    monthly_rainfall = match(site$seasonality$monthly_rainfall, sub_site$site),
+    fourier_prediction = match(site$seasonality$fourier_prediction, sub_site$site)
+  )
+  sub_site$blood_disorders <- match(site$blood_disorders, sub_site$site)
+  sub_site$accessibility <- match(site$accessibility, sub_site$site)
+  sub_site$eir <- match(site$eir, site_filter)
+  return(sub_site)
+}
 
-  site <- site_file
-  for(level in to_mod){
-    mc <- intersect(colnames(index_site), colnames(site[[level]]))
-    site[[level]] <- dplyr::left_join(index_site, site[[level]], by = mc)
+#' Matched join
+#'
+#' @param x Site file element
+#' @param y Data.frame to match for
+#'
+#' @return Site file element filtered by y
+match <- function(x, y){
+  by_names <- names(y)[names(y) %in% names(x)]
+  if(length(by_names) == 0){
+    return(x)
   }
-  return(site)
+  y <- y[, by_names, drop = FALSE]
+  y |>
+    dplyr::left_join(x, by = by_names)
 }
diff --git a/R/site_parameters.R b/R/site_parameters.R
index faf294d..ae23639 100644
--- a/R/site_parameters.R
+++ b/R/site_parameters.R
@@ -18,7 +18,6 @@ site_parameters <- function(interventions, demography, vectors, seasonality,
                             eir = NULL, overrides = list(), burnin = 0){
 
   p <- malariasimulation::get_parameters(overrides = overrides)
-  p$individual_mosquitoes <- FALSE
 
   p$burnin <- 0
   if(burnin > 0){
diff --git a/R/test-helpers.R b/R/test-helpers.R
new file mode 100644
index 0000000..dc3b028
--- /dev/null
+++ b/R/test-helpers.R
@@ -0,0 +1,60 @@
+local_orderly_root <- function(..., .local_envir = parent.frame()) {
+  root <- withr::local_tempdir(.local_envir = .local_envir)
+  suppressMessages(orderly2::orderly_init(root, ...))
+  root
+}
+
+#' Configure the site package to use a temporary cache directory and upstream.
+#'
+#' This setup is suitable for unit testing the package and is automatically torn
+#' down when the caller's frame exits.
+#'
+#' @return the path to the orderly root used by the package as the upstream.
+#' @noRd
+local_test_setup <- function(.local_envir = parent.frame()) {
+  cache <- withr::local_tempdir(.local_envir = .local_envir)
+  upstream <- local_orderly_root(.local_envir = .local_envir)
+
+  withr::local_envvar(R_USER_CACHE_DIR = cache, .local_envir = .local_envir)
+  withr::local_options(
+    .local_envir = .local_envir,
+    "site.orderly_location" = list(type = "path",
+                                   args = list(path = upstream)))
+
+  upstream
+}
+
+#' Create an orderly packet with the given parameters and contents.
+#'
+#' Orderly doesn't have an easy way to craft packets from scratch, so we have
+#' to resort to creating a report directory, writing the files into it,
+#' generating a <name>.R file with a call to `orderly_parameters` in it and
+#' running the report.
+#'
+#' @param name the name used for the packet
+#' @param files a named list of files and their contents to include in the
+#'   report. If the filename ends in `.rds`, the value is written with
+#'   [saveRDS], otherwise [writeLines] is used.
+#' @param parameters a named list of parameters to attach to the packet.
+#' @param root the orderly root in which the packet will be created.
+#' @return the packet id.
+#' @noRd
+create_orderly_packet <- function(name, files, parameters = list(), root) {
+  src <- fs::dir_create(root, "src", name)
+  withr::defer(fs::dir_delete(src))
+
+  args <- paste0(sprintf("%s = NULL", names(parameters)), collapse=",")
+  code <- sprintf("orderly2::orderly_parameters(%s)", args)
+  writeLines(code, fs::path(src, sprintf("%s.R", name)))
+
+  for (i in seq_along(files)) {
+    f <- fs::path(src, names(files)[[i]])
+    if (fs::path_ext(f) == "rds") {
+      saveRDS(files[[i]], f)
+    } else {
+      writeLines(files[[i]], f)
+    }
+  }
+  suppressMessages(orderly2::orderly_run(name, parameters, echo = FALSE,
+                                         root = root))
+}
diff --git a/R/vectors.R b/R/vectors.R
index 8026862..f94b322 100644
--- a/R/vectors.R
+++ b/R/vectors.R
@@ -6,7 +6,7 @@
 #' @return modified parameter list
 add_vectors <- function(p, vectors){
   bionomics <- vectors[ , c("species", "blood_meal_rates", "foraging_time",
-                          "Q0", "phi_bednets", "phi_indoors", "mum")]
+                            "Q0", "phi_bednets", "phi_indoors", "mum")]
 
   species <- list()
   for(s in 1:nrow(bionomics)){
@@ -16,7 +16,8 @@ add_vectors <- function(p, vectors){
   p <- malariasimulation::set_species(
     parameters = p,
     species = species,
-    proportions = vectors$prop)
+    proportions = vectors$prop
+  )
 
   return(p)
 }
diff --git a/README.Rmd b/README.Rmd
index ccac269..4423033 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -20,7 +20,8 @@ knitr::opts_chunk$set(
 [![R-CMD-check](https://github.com/mrc-ide/site/workflows/R-CMD-check/badge.svg)](https://github.com/mrc-ide/site/actions)
 <!-- badges: end -->
 
-Site acts as a translational layer between a site file and malariasimulation.
+Site acts as a translational layer between a site file and malariasimulation: 
+🌍 ➡️ 📉
 
 The site file is the file storing all of the context specific information for a site,
 such as historical intervention coverage, seasonality, vectors etc. To simulate
@@ -31,14 +32,16 @@ If we have a correctly configured site file `example_site`, then all we need to
 do is create the parameter list and pass that to malariasimulation to run:
 
 ```{r, eval=FALSE}
-site <- example_site
+site <- subset_site(example_site, example_site$eir[1,])
 site_par <- site_parameters(
   interventions = site$interventions,
   demography = site$demography,
-  vectors = site$vectors,
-  seasonality = site$seasonality,
-  eir = site$eir,
-  overrides = list(human_population = 10000)
+  vectors = site$vectors$vector_species,
+  seasonality = site$seasonality$seasonality_parameters,
+  eir = site$eir$eir,
+  overrides = list(
+    human_population = 1000
+  )
 )
 site_sim <- malariasimulation::run_simulation(
   timesteps = site_par$timesteps,
diff --git a/README.md b/README.md
index 223f712..1b92104 100644
--- a/README.md
+++ b/README.md
@@ -12,7 +12,7 @@ public.](https://www.repostatus.org/badges/latest/wip.svg)](https://www.repostat
 <!-- badges: end -->
 
 Site acts as a translational layer between a site file and
-malariasimulation.
+malariasimulation: 🌍 ➡️ 📉
 
 The site file is the file storing all of the context specific
 information for a site, such as historical intervention coverage,
@@ -25,14 +25,16 @@ need to do is create the parameter list and pass that to
 malariasimulation to run:
 
 ``` r
-site <- example_site
+site <- subset_site(example_site, example_site$eir[1,])
 site_par <- site_parameters(
   interventions = site$interventions,
   demography = site$demography,
-  vectors = site$vectors,
-  seasonality = site$seasonality,
-  eir = site$eir,
-  overrides = list(human_population = 10000)
+  vectors = site$vectors$vector_species,
+  seasonality = site$seasonality$seasonality_parameters,
+  eir = site$eir$eir,
+  overrides = list(
+    human_population = 1000
+  )
 )
 site_sim <- malariasimulation::run_simulation(
   timesteps = site_par$timesteps,
diff --git a/data-raw/create_example_site.R b/data-raw/create_example_site.R
new file mode 100644
index 0000000..94fb2c8
--- /dev/null
+++ b/data-raw/create_example_site.R
@@ -0,0 +1,11 @@
+example_site <- readRDS("c:/users/pwinskil/desktop/site.rds")
+example_site$eir$eir <- 5
+
+example_site$shape <- lapply(example_site$shape, function(x){
+  data.frame()
+})
+example_site$population$population_by_age <- data.frame()
+
+format(object.size(example_site), "Mb")
+
+usethis::use_data(example_site, overwrite = TRUE)
diff --git a/data/example_site.rda b/data/example_site.rda
index 6f7d172..b485fa7 100644
Binary files a/data/example_site.rda and b/data/example_site.rda differ
diff --git a/man/add_interventions.Rd b/man/add_interventions.Rd
index fcda1e0..5e8b479 100644
--- a/man/add_interventions.Rd
+++ b/man/add_interventions.Rd
@@ -11,8 +11,8 @@ add_interventions(p, interventions, species)
 
 \item{interventions}{site intervention inputs}
 
-\item{species}{Can be falciparum: "pf" or vivax: "pv", for vivax SMC, RTSS
-and PMC are not implemented}
+\item{species}{Can be falciparum: "pf" or vivax: "pv", for vivax SMC, RTSS,
+R21 and PMC are not implemented}
 }
 \value{
 modified parameter list
diff --git a/man/add_r21.Rd b/man/add_r21.Rd
new file mode 100644
index 0000000..0fd8697
--- /dev/null
+++ b/man/add_r21.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/interventions.R
+\name{add_r21}
+\alias{add_r21}
+\title{Add RTS,S}
+\usage{
+add_r21(p, interventions)
+}
+\arguments{
+\item{p}{parameter list}
+
+\item{interventions}{site intervention inputs}
+}
+\value{
+modified parameter list
+}
+\description{
+Add RTS,S
+}
diff --git a/man/adjust_carrying_capacity.Rd b/man/adjust_carrying_capacity.Rd
new file mode 100644
index 0000000..efc2362
--- /dev/null
+++ b/man/adjust_carrying_capacity.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/interventions.R
+\name{adjust_carrying_capacity}
+\alias{adjust_carrying_capacity}
+\title{Adjust carrying capacity}
+\usage{
+adjust_carrying_capacity(p, interventions)
+}
+\arguments{
+\item{p}{parameter list}
+
+\item{interventions}{site intervention inputs}
+}
+\value{
+modified parameter list
+}
+\description{
+Adjust carrying capacity
+}
diff --git a/man/burnin_interventions.Rd b/man/burnin_interventions.Rd
index 0038791..ecdaced 100644
--- a/man/burnin_interventions.Rd
+++ b/man/burnin_interventions.Rd
@@ -15,5 +15,5 @@ burnin_interventions(interventions, burnin)
 Intervention inputs with burn in
 }
 \description{
-Assumes interventions for burn in period are = to those in year 1.
+Assumes interventions for burn in period are = 0 coverage.
 }
diff --git a/man/example_site.Rd b/man/example_site.Rd
index 22c3130..63de602 100644
--- a/man/example_site.Rd
+++ b/man/example_site.Rd
@@ -7,15 +7,21 @@
 \format{
 A list with 10 variables:
 \describe{
-  \item{country}{The iso3c country code}
-  \item{level}{The level of subnational disaggregation}
+  \item{country}{The country name}
+  \item{version}{Site file version}
+  \item{admin_level}{The levels of spatial disaggregation}
   \item{sites}{Unique sites}
-  \item{epi}{Epidemiological site data}
+  \item{shape}{Sptial boundaries (set to NA in testing shapefile)}
+  \item{cases_deaths}{Epidemiological site data}
+  \item{prevalence}{Prevalence data}
   \item{interventions}{Intervention coverage and specification}
-  \item{population}{Population and population at risk projections}
+  \item{population}{Population, population at risk mand age-disaggregated population projections}
   \item{demography}{Demographic projections}
-  \item{vectors}{Vector proportions}
-  \item{seasonality}{Seasonal profile parameters}
+  \item{vectors}{Vector specieis and pyrethroid resistance}
+  \item{seasonality}{Seasonal profile parameters, monthly rainfall and fourier predictions}
+  \item{blood_disorders}{Blood disorder data}
+  \item{accessibility}{Accessibility data}
+  \item{eir}{Calibrated eir}
 }
 }
 \usage{
diff --git a/man/fetch_files.Rd b/man/fetch_files.Rd
new file mode 100644
index 0000000..2cd3936
--- /dev/null
+++ b/man/fetch_files.Rd
@@ -0,0 +1,32 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/fetch.R
+\name{fetch_files}
+\alias{fetch_files}
+\title{Get files from the malariaverse sitefile server.}
+\usage{
+fetch_files(name, parameters, dest, files, expr = NULL)
+}
+\arguments{
+\item{name}{The name of the orderly report.}
+
+\item{parameters}{A named list of parameters to use when searching for the
+orderly packet. If a query expression `expr` is specified, these parameters
+are substituted into the query using the this: prefix. If no expression is
+specified, the latest packet matching these parameters exactly is selected.}
+
+\item{dest}{A directory into which the files should be copied.}
+
+\item{files}{An optionally-named character vector of files to copy from the
+packet and into the destination directory. If the vector is named, these
+names are used as the destination file path.}
+
+\item{expr}{The query expression to filter packets. This may be an arbitrary
+orderly query, including a literal packet ID. If absent or NULL, the
+specified list of parameters is used and matched exactly.}
+}
+\value{
+the id of the orderly packet the files were copied from.
+}
+\description{
+Get files from the malariaverse sitefile server.
+}
diff --git a/man/fetch_site.Rd b/man/fetch_site.Rd
new file mode 100644
index 0000000..cbb0706
--- /dev/null
+++ b/man/fetch_site.Rd
@@ -0,0 +1,44 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/fetch.R
+\name{fetch_site}
+\alias{fetch_site}
+\title{Fetch a site file for a country from the malariaverse sitefile server.}
+\usage{
+fetch_site(
+  iso3c = NULL,
+  version = NULL,
+  admin_level = NULL,
+  urban_rural = NULL,
+  id = NULL
+)
+}
+\arguments{
+\item{iso3c}{the ISO country code, a scalar character.}
+
+\item{version}{the dataset version, a scalar character.}
+
+\item{admin_level}{a scalar number.}
+
+\item{urban_rural}{a scalar logical.}
+
+\item{id}{a packet ID used to select an exact packet.}
+}
+\value{
+The contents of the site file.
+}
+\description{
+The site file is identified by its country code, and optionally the
+admin_level, urban/rural setting and version of the site files. The latest
+packet from the server matching these parameters is used.
+}
+\details{
+Alternatively, a packet ID can be specified in order to pick an exact file
+set.
+}
+\examples{
+\dontrun{
+fetch_site("NGA")
+fetch_site("NGA", admin_level = 1)
+fetch_site(id = "20240801-062621-6f95851a")
+}
+}
diff --git a/man/match.Rd b/man/match.Rd
new file mode 100644
index 0000000..82902df
--- /dev/null
+++ b/man/match.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/single_site.R
+\name{match}
+\alias{match}
+\title{Matched join}
+\usage{
+match(x, y)
+}
+\arguments{
+\item{x}{Site file element}
+
+\item{y}{Data.frame to match for}
+}
+\value{
+Site file element filtered by y
+}
+\description{
+Matched join
+}
diff --git a/man/single_site.Rd b/man/single_site.Rd
deleted file mode 100644
index 82c5269..0000000
--- a/man/single_site.Rd
+++ /dev/null
@@ -1,19 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/single_site.R
-\name{single_site}
-\alias{single_site}
-\title{Extract a single site-input from a country site file}
-\usage{
-single_site(site_file, index)
-}
-\arguments{
-\item{site_file}{Country site file}
-
-\item{index}{Index row from site_file$sites}
-}
-\value{
-Single site
-}
-\description{
-Extract a single site-input from a country site file
-}
diff --git a/man/subset_site.Rd b/man/subset_site.Rd
new file mode 100644
index 0000000..0f1a9ae
--- /dev/null
+++ b/man/subset_site.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/single_site.R
+\name{subset_site}
+\alias{subset_site}
+\title{Extract a subset from a country site file}
+\usage{
+subset_site(site, site_filter)
+}
+\arguments{
+\item{site}{Country site file}
+
+\item{site_filter}{Data.frame to filter site file elements by. Filtering will
+be conducted on any matched columns for each element.}
+}
+\value{
+Filtered site file
+}
+\description{
+Extract a subset from a country site file
+}
diff --git a/tests/testthat/test-add_time.R b/tests/testthat/test-add_time.R
index 0c44897..613eff1 100644
--- a/tests/testthat/test-add_time.R
+++ b/tests/testthat/test-add_time.R
@@ -1,5 +1,5 @@
 test_that("multiplication works", {
-  example_site <- single_site(example_site, 1)
+  example_site <- subset_site(example_site, example_site$eir[1,])
 
   p0 <- list()
   interventions <- example_site$interventions
diff --git a/tests/testthat/test-burnin.R b/tests/testthat/test-burnin.R
index 9e587bd..90e37dd 100644
--- a/tests/testthat/test-burnin.R
+++ b/tests/testthat/test-burnin.R
@@ -1,6 +1,6 @@
 test_that("burnin works", {
   bi <- 10
-  example_site <- single_site(example_site, 1)
+  example_site <- subset_site(example_site, example_site$eir[1,])
 
   interventions <- burnin_interventions(example_site$interventions, bi)
   expect_equal(nrow(interventions), nrow(example_site$interventions) + bi)
@@ -13,8 +13,8 @@ test_that("burnin works", {
   p <- site_parameters(
     interventions = example_site$interventions,
     demography = example_site$demography,
-    vectors = example_site$vectors,
-    seasonality = example_site$seasonality,
+    vectors = example_site$vectors$vector_species,
+    seasonality = example_site$seasonality$seasonality_parameters,
     burnin = bi
   )
   expect_identical(p$burnin, bi)
diff --git a/tests/testthat/test-carrying_capacity.R b/tests/testthat/test-carrying_capacity.R
new file mode 100644
index 0000000..e015335
--- /dev/null
+++ b/tests/testthat/test-carrying_capacity.R
@@ -0,0 +1,16 @@
+test_that("Carrying capacity works", {
+  single_site <- subset_site(example_site, example_site$eir[1,])
+  interventions <- single_site$interventions
+  interventions$lsm_cov[1:10] <- 0.5
+  p0 <- malariasimulation::get_parameters()
+  p0$baseline_year <- 2000
+
+  p1 <- p0 |>
+    adjust_carrying_capacity(
+      interventions = interventions
+    )
+
+  expect_true(p1$carrying_capacity)
+  expect_equal(p1$carrying_capacity_scalers, matrix(1 - interventions$lsm_cov, ncol = 1))
+  expect_equal(p1$carrying_capacity_timesteps, 1 + (interventions$year - p1$baseline_year) * 365)
+})
diff --git a/tests/testthat/test-demography.R b/tests/testthat/test-demography.R
index ec16957..87f2b43 100644
--- a/tests/testthat/test-demography.R
+++ b/tests/testthat/test-demography.R
@@ -1,5 +1,5 @@
 test_that("adding (static) demography works", {
-  example_site <- single_site(example_site, 1)
+  example_site <- subset_site(example_site, example_site$eir[1,])
   demography <- example_site$demography
   p0 <- malariasimulation::get_parameters()
   p0$baseline_year <- 2000
@@ -8,7 +8,7 @@ test_that("adding (static) demography works", {
   expect_equal(p1$deathrate_agegroups, round(unique(demography$age_upper) * 365))
   expect_equal(p1$deathrate_timesteps, 365 * (unique(demography$year) - p0$baseline_year))
   expect_equal(p1$deathrates,
-               matrix(demography$mortality_rate / 365,
+               matrix(demography$adjusted_mortality_rates / 365,
                       nrow = length(unique(demography$year)),
                       byrow = TRUE))
 })
diff --git a/tests/testthat/test-fetch.R b/tests/testthat/test-fetch.R
new file mode 100644
index 0000000..b6a189d
--- /dev/null
+++ b/tests/testthat/test-fetch.R
@@ -0,0 +1,202 @@
+test_that("has a default location", {
+  cfg <- location_configuration()
+  expect_equal(cfg$type, "packit")
+  expect_match(cfg$args$url, "^https?://")
+})
+
+test_that("can override github token of default location", {
+  withr::local_envvar(GITHUB_TOKEN="mytoken")
+  cfg <- location_configuration()
+
+  expect_equal(cfg$type, "packit")
+  expect_match(cfg$args$url, "^https?://")
+  expect_equal(cfg$args$token, "mytoken")
+})
+
+test_that("can override default location", {
+  cfg <- list(type = "path", args = list(path = "path/to/location"))
+  withr::local_options("site.orderly_location" = cfg)
+  expect_equal(location_configuration(), cfg)
+})
+
+test_that("can configure orderly repository", {
+  local_test_setup()
+
+  root <- suppressMessages(configure_orderly())
+  expect_true(fs::dir_exists(root))
+  expect_contains(orderly2::orderly_location_list(root = root), LOCATION_NAME)
+})
+
+test_that("can change orderly location", {
+  local_test_setup()
+  upstream1 <- local_orderly_root()
+  upstream2 <- local_orderly_root()
+
+  withr::local_options("site.orderly_location" = list(
+    type = "path",
+    args = list(path = upstream1)))
+
+  root <- suppressMessages(configure_orderly())
+  res <- orderly2::orderly_location_list(root = root, verbose = TRUE)
+  expect_equal(res[[which(res$name == LOCATION_NAME), "args"]]$path, upstream1)
+
+  withr::local_options("site.orderly_location" = list(
+    type = "path",
+    args = list(path = upstream2)))
+
+  root <- configure_orderly()
+  res <- orderly2::orderly_location_list(root = root, verbose = TRUE)
+
+  expect_equal(res[[which(res$name == LOCATION_NAME), "args"]]$path, upstream2)
+})
+
+
+test_that("can fetch files from upstream", {
+  upstream <- local_test_setup()
+
+  id <- create_orderly_packet("data", list("data.txt" = "Contents"),
+                              root = upstream)
+
+  dest <- withr::local_tempdir()
+  res <- suppressMessages(fetch_files("data", list(), dest, "data.txt"))
+
+  expect_equal(readLines(fs::path(dest, "data.txt")), "Contents")
+  expect_equal(id, res)
+})
+
+
+test_that("files are cached and only downloaded once", {
+  upstream <- local_test_setup()
+
+  id <- create_orderly_packet("data", list("data.txt" = "Contents"),
+                              root = upstream)
+
+  dest <- withr::local_tempdir()
+  suppressMessages({
+    expect_message(fetch_files("data", list(), dest,
+                               c("first.txt" = "data.txt")),
+                   "Need to fetch 1 file")
+
+    expect_no_message(fetch_files("data", list(), dest,
+                                  c("second.txt" = "data.txt")),
+                      message = "Need to fetch 1 file")
+  })
+
+  expect_equal(readLines(fs::path(dest, "first.txt")), "Contents")
+  expect_equal(readLines(fs::path(dest, "second.txt")), "Contents")
+})
+
+
+test_that("can filter based on parameters", {
+  upstream <- local_test_setup()
+
+  create_orderly_packet("data", list("data.txt" = "Hello"),
+                        parameters = list(code = "foo"),
+                        root = upstream)
+
+  create_orderly_packet("data", list("data.txt" = "World"),
+                        parameters = list(code = "bar"),
+                        root = upstream)
+
+  dest <- withr::local_tempdir()
+
+  suppressMessages(fetch_files("data", list(code = "foo"),
+                               dest, c("hello.txt" = "data.txt")))
+
+  suppressMessages(fetch_files("data", list(code = "bar"),
+                               dest, c("world.txt" = "data.txt")))
+
+  suppressMessages(fetch_files("data", list(), dest,
+                               c("latest.txt" = "data.txt")))
+
+  expect_equal(readLines(fs::path(dest, "hello.txt")), "Hello")
+  expect_equal(readLines(fs::path(dest, "world.txt")), "World")
+  expect_equal(readLines(fs::path(dest, "latest.txt")), "World")
+})
+
+
+test_that("can filter based on packet ID", {
+  upstream <- local_test_setup()
+
+  id1 <- create_orderly_packet("data", list("data.txt" = "Hello"),
+                               parameters = list(code = "foo"),
+                               root = upstream)
+
+  id2 <- create_orderly_packet("data", list("data.txt" = "World"),
+                               parameters = list(code = "bar"),
+                               root = upstream)
+
+  dest <- withr::local_tempdir()
+
+  suppressMessages(fetch_files("data", list(), expr = id1,
+                               dest, c("hello.txt" = "data.txt")))
+
+  suppressMessages(fetch_files("data", list(), expr = id2,
+                               dest, c("world.txt" = "data.txt")))
+
+  expect_equal(readLines(fs::path(dest, "hello.txt")), "Hello")
+  expect_equal(readLines(fs::path(dest, "world.txt")), "World")
+})
+
+
+test_that("can filter based on an arbitrary query", {
+  upstream <- local_test_setup()
+
+  id1 <- create_orderly_packet("data", list("data.txt" = "Hello"),
+                               parameters = list(foo = 1, bar = 1),
+                               root = upstream)
+
+  id2 <- create_orderly_packet("data", list("data.txt" = "World"),
+                               parameters = list(foo = 1, bar = 2),
+                               root = upstream)
+
+  dest <- withr::local_tempdir()
+
+  suppressMessages(fetch_files(
+    expr = "single(parameter:foo == parameter:bar)",
+    "data", list(), dest, c("hello.txt" = "data.txt")))
+
+  suppressMessages(fetch_files(
+    expr = "single(parameter:foo != parameter:bar)",
+    "data", list(), dest, c("world.txt" = "data.txt")))
+
+  expect_equal(readLines(fs::path(dest, "hello.txt")), "Hello")
+  expect_equal(readLines(fs::path(dest, "world.txt")), "World")
+})
+
+
+test_that("can fetch site", {
+  upstream <- local_test_setup()
+
+  make <- function(iso3c, admin_level, version, value) {
+    create_orderly_packet(
+      "calibration_diagnostics",
+      list("calibrated_scaled_site.rds" = list(value = value)),
+      parameters = list(iso3c = iso3c, version = version,
+                        admin_level = admin_level),
+      root = upstream)
+  }
+  fetch <- function(...) suppressMessages(fetch_site(...))
+
+  id1 <- make("TGO", 1, "v1", value = 1)
+  id2 <- make("TGO", 2, "v1", value = 2)
+  id3 <- make("NGA", 1, "v1", value = 3)
+  id4 <- make("NGA", 1, "v2", value = 4)
+
+  expect_equal(fetch("TGO"), list(value = 2))
+  expect_equal(fetch("NGA"), list(value = 4))
+
+  expect_equal(fetch("TGO", admin_level = 1), list(value = 1))
+  expect_equal(fetch("TGO", admin_level = 2), list(value = 2))
+
+  expect_equal(fetch("NGA", version = "v1"), list(value = 3))
+  expect_equal(fetch("NGA", version = "v2"), list(value = 4))
+
+  expect_equal(fetch(id = id1), list(value = 1))
+  expect_equal(fetch(id = id2), list(value = 2))
+  expect_equal(fetch(id = id3), list(value = 3))
+  expect_equal(fetch(id = id4), list(value = 4))
+
+  expect_error(fetch("TGO", id = id1),
+               "Exactly one of `iso3c` and `id` must be supplied")
+})
diff --git a/tests/testthat/test-irs.R b/tests/testthat/test-irs.R
index f1a8e3e..6ae5522 100644
--- a/tests/testthat/test-irs.R
+++ b/tests/testthat/test-irs.R
@@ -1,19 +1,20 @@
 test_that("adding irs works", {
-  example_site <- single_site(example_site, 1)
-  interventions <- example_site$interventions
+  single_site <- subset_site(example_site, example_site$eir[1,])
+  interventions <- single_site$interventions
   interventions$irs_cov[1:10] <- 0.5
 
   p0 <- malariasimulation::get_parameters()
   p0$baseline_year <- 2000
-  p0 <- add_seasonality(p0, example_site$seasonality)
+  p0 <- add_seasonality(p0, single_site$seasonality$seasonality_parameters)
   p0$g <- p0$g + 5
   p1 <- add_irs(
     p = p0,
-    interventions = interventions)
+    interventions = interventions
+  )
 
   month <- 365 / 12
   peak <- malariasimulation::peak_season_offset(p1)
-  year_start_times <-  1 + (example_site$interventions$year - p1$baseline_year) * 365
+  year_start_times <-  1 + (single_site$interventions$year - p1$baseline_year) * 365
   peak_season_times <- peak + year_start_times
   # Assume IRS occurs 3 months before seasonal peak
   timesteps <- round(peak_season_times - 3 * month)
@@ -31,20 +32,21 @@ test_that("adding irs works", {
   expect_equal(p1$spraying_ms_theta, matrix(interventions$ms_theta[!index], ncol = 1))
 
   # With negative timestep early relative to peak
-  example_site <- single_site(example_site, 1)
-  interventions <- example_site$interventions
+  single_site <- subset_site(example_site, example_site$eir[1,])
+  interventions <- single_site$interventions
   interventions$irs_cov[1:10] <- 0.5
 
   p0 <- malariasimulation::get_parameters()
   p0$baseline_year <- 2000
-  p0 <- add_seasonality(p0, example_site$seasonality)
+  p0 <- add_seasonality(p0, single_site$seasonality$seasonality_parameters)
   p1 <- add_irs(
     p = p0,
-    interventions = interventions)
+    interventions = interventions
+  )
 
   month <- 365 / 12
   peak <- malariasimulation::peak_season_offset(p1)
-  year_start_times <-  1 + (example_site$interventions$year - p1$baseline_year) * 365
+  year_start_times <-  1 + (single_site$interventions$year - p1$baseline_year) * 365
   peak_season_times <- peak + year_start_times
   # Assume IRS occurs 3 months before seasonal peak
   timesteps <- round(peak_season_times - 3 * month)
@@ -62,3 +64,39 @@ test_that("adding irs works", {
   expect_equal(p1$spraying_ms_theta, matrix(interventions$ms_theta[!index], ncol = 1))
 
 })
+
+test_that("adding irs multiple rounds works", {
+  single_site <- subset_site(example_site, example_site$eir[1,])
+  interventions <- single_site$interventions
+  interventions$irs_spray_rounds[nrow(interventions)] <- 2
+  interventions$irs_cov[1:23] <- 0.5
+
+  p0 <- malariasimulation::get_parameters()
+  p0$baseline_year <- 2000
+  p0 <- add_seasonality(p0, single_site$seasonality$seasonality_parameters)
+  p0$g <- p0$g + 5
+  p1 <- add_irs(
+    p = p0,
+    interventions = interventions
+  )
+
+  month <- 365 / 12
+  peak <- malariasimulation::peak_season_offset(p1)
+  year_start_times <-  1 + (single_site$interventions$year - p1$baseline_year) * 365
+  peak_season_times <- peak + year_start_times
+  # Assume IRS occurs 3 months before seasonal peak
+  timesteps <- round(peak_season_times - 3 * month)
+  timesteps <- round(c(timesteps, peak_season_times[length(peak_season_times)] + 3 * month))
+  index <- timesteps < 0
+  timesteps <- timesteps[!index]
+
+  expect_equal(p1$spraying, TRUE)
+  expect_equal(p1$spraying_timesteps, timesteps)
+  expect_equal(p1$spraying_coverages, rep(interventions$irs_cov, interventions$irs_spray_rounds)[!index])
+  expect_equal(p1$spraying_ls_theta, matrix(rep(interventions$ls_theta, interventions$irs_spray_rounds)[!index], ncol = 1))
+  expect_equal(p1$spraying_ls_gamma, matrix(rep(interventions$ls_gamma, interventions$irs_spray_rounds)[!index], ncol = 1))
+  expect_equal(p1$spraying_ks_theta, matrix(rep(interventions$ks_theta, interventions$irs_spray_rounds)[!index], ncol = 1))
+  expect_equal(p1$spraying_ks_theta, matrix(rep(interventions$ks_theta, interventions$irs_spray_rounds)[!index], ncol = 1))
+  expect_equal(p1$spraying_ms_theta, matrix(rep(interventions$ms_theta, interventions$irs_spray_rounds)[!index], ncol = 1))
+  expect_equal(p1$spraying_ms_theta, matrix(rep(interventions$ms_theta, interventions$irs_spray_rounds)[!index], ncol = 1))
+})
diff --git a/tests/testthat/test-itns.R b/tests/testthat/test-itns.R
index 37ed219..17a661e 100644
--- a/tests/testthat/test-itns.R
+++ b/tests/testthat/test-itns.R
@@ -1,19 +1,30 @@
 test_that("adding itns works", {
-  example_site <- single_site(example_site, 1)
-  interventions <- example_site$interventions
+  single_site <- subset_site(example_site, example_site$eir[1,])
+  interventions <- single_site$interventions
   p0 <- malariasimulation::get_parameters()
   p0$baseline_year <- 2000
   p1 <- add_itns(
     p = p0,
-    interventions = interventions)
+    interventions = interventions
+  )
 
   expect_equal(p1$bednets, TRUE)
-  expect_equal(p1$bednet_timesteps, 1 + (interventions$year - p0$baseline_year) * 365)
+  expect_equal(p1$bednet_timesteps, interventions$itn_distribution_day + (interventions$year - p0$baseline_year) * 365)
   expect_equal(p1$bednet_coverages, interventions$itn_input_dist)
-  expect_equal(p1$bednet_retention, 365 * 5)
+  expect_equal(p1$bednet_retention, unique(interventions$mean_retention))
   expect_equal(p1$bednet_dn0,
                matrix(rep(interventions$dn0, length(p1$species)), ncol = length(p1$species)))
   expect_equal(p1$bednet_rn, matrix(rep(interventions$rn0, length(p1$species)), ncol = length(p1$species)))
   expect_equal(p1$bednet_rnm, matrix(rep(interventions$rnm, length(p1$species)), ncol = length(p1$species)))
   expect_equal(p1$bednet_gamman, interventions$gamman * 365)
+
+
+  interventions$mean_retention[1:2] <- 1:2
+  expect_error(
+    add_itns(
+      p = p0,
+      interventions = interventions
+    ),
+    "Time-varying net rentetion is not currently supported"
+  )
 })
diff --git a/tests/testthat/test-pmc.R b/tests/testthat/test-pmc.R
index 92ba671..aeb2d22 100644
--- a/tests/testthat/test-pmc.R
+++ b/tests/testthat/test-pmc.R
@@ -1,6 +1,6 @@
 test_that("PMC works", {
-  example_site <- single_site(example_site, 1)
-  interventions <- example_site$interventions
+  single_site <- subset_site(example_site, example_site$eir[1,])
+  interventions <- single_site$interventions
   interventions$pmc_cov[1:10] <- 0.5
   p0 <- malariasimulation::get_parameters()
   p0$baseline_year <- 2000
diff --git a/tests/testthat/test-rtss.R b/tests/testthat/test-rtss.R
deleted file mode 100644
index 2c0c556..0000000
--- a/tests/testthat/test-rtss.R
+++ /dev/null
@@ -1,19 +0,0 @@
-test_that("adding rtss works", {
-  example_site <- single_site(example_site, 1)
-  interventions <- example_site$interventions
-  interventions$rtss_cov[1:10] <- 0.5
-  p0 <- malariasimulation::get_parameters()
-  p0$baseline_year <- 2000
-  p1 <- add_rtss(
-    p = p0,
-    interventions = interventions)
-
-  month <- 365 / 12
-  expect_equal(p1$rtss_epi_age, round(6 * month))
-  expect_equal(p1$rtss_epi_booster_coverage, 0.8)
-  expect_equal(p1$rtss_epi_boosters, round(18 * month))
-  expect_equal(p1$rtss_epi_timesteps, 1 + (365 * (interventions$year - p0$baseline_year)))
-  expect_equal(p1$rtss_epi_coverages, interventions$rtss_cov)
-  expect_equal(p1$rtss_epi_min_wait, 0)
-  expect_equal(p1$rtss_epi_seasonal_boosters, FALSE)
-})
diff --git a/tests/testthat/test-seasonality.R b/tests/testthat/test-seasonality.R
index 27ca5ab..7307207 100644
--- a/tests/testthat/test-seasonality.R
+++ b/tests/testthat/test-seasonality.R
@@ -1,6 +1,6 @@
 test_that("adding seasonality works", {
-  example_site <- single_site(example_site, 1)
-  seasonality <- example_site$seasonality
+  single_site <- subset_site(example_site, example_site$eir[1,])
+  seasonality <- single_site$seasonality$seasonality_parameters
   p0 <- malariasimulation::get_parameters()
   p1 <- add_seasonality(p = p0, seasonality = seasonality)
 
diff --git a/tests/testthat/test-site_parameters.R b/tests/testthat/test-site_parameters.R
index 4fe2b24..9c98b5d 100644
--- a/tests/testthat/test-site_parameters.R
+++ b/tests/testthat/test-site_parameters.R
@@ -1,18 +1,18 @@
 test_that("site parameters wrapper works", {
-  example_site <- single_site(example_site, 1)
+  single_site <- subset_site(example_site, example_site$eir[1,])
   p <- site_parameters(
-    interventions = example_site$interventions,
-    demography = example_site$demography,
-    vectors = example_site$vectors,
-    seasonality = example_site$seasonality
+    interventions = single_site$interventions,
+    demography = single_site$demography,
+    vectors = single_site$vectors$vector_species,
+    seasonality = single_site$seasonality$seasonality_parameters
   )
   expect_type(p, "list")
 
   p <- site_parameters(
-    interventions = example_site$interventions,
-    demography = example_site$demography,
-    vectors = example_site$vectors,
-    seasonality = example_site$seasonality,
+    interventions = single_site$interventions,
+    demography = single_site$demography,
+    vectors = single_site$vectors$vector_species,
+    seasonality = single_site$seasonality$seasonality_parameters,
     eir = 10
   )
   expect_type(p, "list")
@@ -20,19 +20,19 @@ test_that("site parameters wrapper works", {
 
 
 test_that("setting vivax works", {
-  example_site <- single_site(example_site, 1)
-  example_site$interventions$rtss_cov <- 0.1
-  example_site$interventions$pmc_cov <- 0.1
-  example_site$interventions$smc_cov <- 0.1
+  single_site <- subset_site(example_site, example_site$eir[1,])
+  single_site$interventions$rtss_cov <- 0.1
+  single_site$interventions$pmc_cov <- 0.1
+  single_site$interventions$smc_cov <- 0.1
 
   p <- site_parameters(
-    interventions = example_site$interventions,
-    demography = example_site$demography,
-    vectors = example_site$vectors,
-    seasonality = example_site$seasonality,
+    interventions = single_site$interventions,
+    demography = single_site$demography,
+    vectors = single_site$vectors$vector_species,
+    seasonality = single_site$seasonality$seasonality_parameters,
     species = "pv"
   )
-  expect_false(p$rtss)
+  expect_false(p$pev)
   expect_false(p$smc)
   expect_false(p$pmc)
 })
diff --git a/tests/testthat/test-smc.R b/tests/testthat/test-smc.R
index 4b72462..44d72d3 100644
--- a/tests/testthat/test-smc.R
+++ b/tests/testthat/test-smc.R
@@ -1,13 +1,14 @@
 test_that("adding smc works", {
-  example_site <- single_site(example_site, 1)
-  interventions <- example_site$interventions
+  single_site <- subset_site(example_site, example_site$eir[1,])
+  interventions <- single_site$interventions
   interventions$smc_cov[1:10] <- 0.5
   p0 <- malariasimulation::get_parameters()
 
   p0$baseline_year <- 2000
   p1 <- add_smc(
     p = p0,
-    interventions = interventions)
+    interventions = interventions
+  )
 
   month <- 365 / 12
 
diff --git a/tests/testthat/test-treatment.R b/tests/testthat/test-treatment.R
index e535f3d..e9333be 100644
--- a/tests/testthat/test-treatment.R
+++ b/tests/testthat/test-treatment.R
@@ -1,13 +1,15 @@
 test_that("adding treatment works", {
-  example_site <- single_site(example_site, 1)
-  interventions <- example_site$interventions
+  single_site <- subset_site(example_site, example_site$eir[1,])
+  interventions <- single_site$interventions
   p0 <- malariasimulation::get_parameters()
   p0$baseline_year <- 2000
   # Drug types
   p0 <- add_drugs(p0)
 
-  p1 <- add_treatment(p = p0,
-                      interventions = interventions)
+  p1 <- add_treatment(
+    p = p0,
+    interventions = interventions
+  )
 
   expect_equal(p1$clinical_treatment_drugs, list(4, 5))
 
@@ -18,7 +20,8 @@ test_that("adding treatment works", {
   expect_equal(p1$clinical_treatment_timesteps,
                list(
                  1 + 365 * (interventions$year - p0$baseline_year),
-                 1 + 365 * (interventions$year - p0$baseline_year))
+                 1 + 365 * (interventions$year - p0$baseline_year)
+               )
   )
 })
 
diff --git a/tests/testthat/test-vaccine.R b/tests/testthat/test-vaccine.R
new file mode 100644
index 0000000..9d11515
--- /dev/null
+++ b/tests/testthat/test-vaccine.R
@@ -0,0 +1,81 @@
+test_that("adding rtss works", {
+  single_site <- subset_site(example_site, example_site$eir[1,])
+  interventions <- single_site$interventions
+  interventions$rtss_cov[1:10] <- 0.5
+  p0 <- malariasimulation::get_parameters()
+  p0$baseline_year <- 2000
+  p1 <- add_rtss(
+    p = p0,
+    interventions = interventions
+  )
+
+  month <- 365 / 12
+  expect_equal(p1$pev_epi_age, round(6 * month))
+  expect_equal(p1$pev_epi_booster_coverage, matrix(rep(0.8, length(p1$pev_epi_timesteps))))
+  expect_equal(p1$pev_epi_booster_spacing, round(12 * month))
+  expect_equal(p1$pev_epi_timesteps, 1 + (365 * (interventions$year - p0$baseline_year)))
+  expect_equal(p1$pev_epi_coverages, interventions$rtss_cov)
+  expect_equal(p1$pev_epi_min_wait, 0)
+  expect_equal(p1$pev_epi_seasonal_boosters, FALSE)
+})
+
+test_that("adding r21 works", {
+  single_site <- subset_site(example_site, example_site$eir[1,])
+  interventions <- single_site$interventions
+  interventions$r21_cov[1:10] <- 0.5
+  p0 <- malariasimulation::get_parameters()
+  p0$baseline_year <- 2000
+  p1 <- add_r21(
+    p = p0,
+    interventions = interventions
+  )
+
+  month <- 365 / 12
+  expect_equal(p1$pev_epi_age, round(6 * month))
+  expect_equal(p1$pev_epi_booster_coverage, matrix(rep(0.8, length(p1$pev_epi_timesteps))))
+  expect_equal(p1$pev_epi_booster_spacing, round(12 * month))
+  expect_equal(p1$pev_epi_timesteps, 1 + (365 * (interventions$year - p0$baseline_year)))
+  expect_equal(p1$pev_epi_coverages, interventions$r21_cov)
+  expect_equal(p1$pev_epi_min_wait, 0)
+  expect_equal(p1$pev_epi_seasonal_boosters, FALSE)
+})
+
+test_that("adding r21 and rtss fails as expected", {
+  single_site <- subset_site(example_site, example_site$eir[1,])
+  interventions <- single_site$interventions
+  interventions$rtss_cov[1:10] <- 0.5
+  interventions$r21_cov[1:10] <- 0.5
+  p0 <- malariasimulation::get_parameters()
+  p0$baseline_year <- 2000
+  expect_warning(
+    p1 <- add_interventions(
+      p = p0,
+      interventions = interventions,
+      species = "pf"
+    )
+  )
+})
+
+test_that("n_doses can be specified", {
+  single_site <- subset_site(example_site, example_site$eir[1,])
+  interventions <- single_site$interventions
+  interventions$rtss_cov[1:10] <- 0.5
+  interventions$n_doses <- c(rep(3,8), 4, 4, rep(0, 15))
+  p0 <- malariasimulation::get_parameters()
+  p0$baseline_year <- 2000
+  p1 <- add_rtss(
+    p = p0,
+    interventions = interventions
+  )
+
+  month <- 365 / 12
+  expect_equal(p1$pev_epi_age, round(6 * month))
+  booster_cov <- rep(0, nrow(single_site$interventions))
+  booster_cov[9:10] <- 0.8
+  expect_equal(p1$pev_epi_booster_coverage, matrix(booster_cov))
+  expect_equal(p1$pev_epi_booster_spacing, round(12 * month))
+  expect_equal(p1$pev_epi_timesteps, 1 + (365 * (interventions$year - p0$baseline_year)))
+  expect_equal(p1$pev_epi_coverages, interventions$rtss_cov)
+  expect_equal(p1$pev_epi_min_wait, 0)
+  expect_equal(p1$pev_epi_seasonal_boosters, FALSE)
+})
diff --git a/tests/testthat/test-vectors.R b/tests/testthat/test-vectors.R
index 78addd7..73e0514 100644
--- a/tests/testthat/test-vectors.R
+++ b/tests/testthat/test-vectors.R
@@ -1,6 +1,6 @@
 test_that("adding vectors works", {
-  example_site <- single_site(example_site, 1)
-  vectors <- example_site$vectors
+  single_site <- subset_site(example_site, example_site$eir[1,])
+  vectors <- single_site$vectors$vector_species
   p0 <- malariasimulation::get_parameters()
   p1 <- add_vectors(p = p0, vectors = vectors)