Simulation

Source Data Download and NWS Case Detection Simulation

Published

May 4, 2026

1 Overview

The script demonstrates how to interactively download model results and supporting information archived on the Open Science Framework (OSF) and then simulate New World Screwworm (NWS) case detections.

The downloaded data includes spatial GeoTiff files depicting potential NWS abundance and estimated spread for each week between January 2022 and February 2026, spatial boundaries for analysis and visualization, and a copy of previously simulated NWS case point locations.


2 Libraries

Load needed libraries and packages.

Show R Code
library(tidyverse) # data wrangling
options(dplyr.summarise.inform = FALSE)
library(here) # directory management
library(pals) # color palettes
library(osfr) # interact with Open Science Framework (OSF, data repository)
library(terra) # spatial data manipulation
library(sf) # spatial data manipulation
library(spatstat) # spatial data manipulation

3 Custom Functions

Load customized functions.

Show R Code
source(here("R/utilities.R"))
source_dir(here("R"))

4 Download OSF Data

Spatial GeoTIFF files (.tif) depicting NWS estimated potential abundance and spread are available on the Open Science Framework. Here we link the location and download the files to a demonstration directory (/demo).

Downloaded Files:
- spread_GeoTIFF_2026-03-24.zip: Zipped folder containing rasters showing estimated NWS spread.
- niche_GeoTIFF_2026-03-24.zip: Zipped folder containing rasters showing potential NWS abundance.
- study_polygons.gpkg: Geopackage containing sf class multipolygons delineating political boundaries.
- sim_locations_2026-03-24.csv: Simulated NWS detections; the same file is reproduced in this script.

Show R Code
# file location on OSF
osf_project_location <- osf_retrieve_node("https://osf.io/uvqmw/")

# Project name on OSF
osf_project_location$name 
[1] "hominivorax-geostat"
Show R Code
# create demo directory
if (!dir.exists("demo")) dir.create("demo", recursive = TRUE)

# get data file names
osf_id <- osf_project_location %>%
  osf_ls_files()

# download all data
osf_download(osf_id,
             path = here("demo"),
             conflicts = "overwrite")
name id local_path meta
spread_GeoTIFF_2026-03-24.zip 69c2d3a4e03dca33dc229717 demo/spread_GeoTIFF_2026-03-24.zip spread_GeoTIFF_2026-03-24.zip , file , /69c2d3a4e03dca33dc229717 , 2465016 , osfstorage , /spread_GeoTIFF_2026-03-24.zip , 1774375844 , 1774375844 , 09634b49a5a4cf7152e86cceaee889e0 , d1b3407a37034b84c5e5e62dd39fce63f265b92ec42b0f7e8d1e4119aee1546c , 10 , FALSE , 1 , FALSE , https://api.osf.io/v2/files/69c2d3a4e03dca33dc229717/ , https://files.osf.io/v1/resources/uvqmw/providers/osfstorage/69c2d3a4e03dca33dc229717 , https://files.osf.io/v1/resources/uvqmw/providers/osfstorage/69c2d3a4e03dca33dc229717 , https://files.osf.io/v1/resources/uvqmw/providers/osfstorage/69c2d3a4e03dca33dc229717 , https://osf.io/download/69c2d3a4e03dca33dc229717/ , https://mfr.osf.io/render?url=https%3A%2F%2Fosf.io%2Fdownload%2F69c2d3a4e03dca33dc229717%2F%3Fdirect%26mode%3Drender, https://osf.io/uvqmw/files/osfstorage/69c2d3a4e03dca33dc229717 , https://api.osf.io/v2/files/69c2d3a4e03dca33dc229717/ , https://api.osf.io/v2/files/69c2b365ae03e5ce6d107b09/ , 69c2b365ae03e5ce6d107b09 , files , https://api.osf.io/v2/files/69c2d3a4e03dca33dc229717/versions/ , https://api.osf.io/v2/nodes/uvqmw/ , uvqmw , nodes , https://api.osf.io/v2/nodes/uvqmw/ , nodes , nodes , uvqmw , https://api.osf.io/v2/files/69c2d3a4e03dca33dc229717/cedar_metadata_records/
niche_GeoTIFF_2026-03-24.zip 69c3204a686690fea0927d64 demo/niche_GeoTIFF_2026-03-24.zip niche_GeoTIFF_2026-03-24.zip , file , /69c3204a686690fea0927d64 , 12205660 , osfstorage , /niche_GeoTIFF_2026-03-24.zip , 1774395466 , 1774395466 , eb0bbf520ebbb15df84d02ad2222e6c3 , 28dd79fcc7f42242a3e9b1b68ba6598875452d77de9a5453e1bf5e791f76baaa , 9 , FALSE , 1 , FALSE , https://api.osf.io/v2/files/69c3204a686690fea0927d64/ , https://files.osf.io/v1/resources/uvqmw/providers/osfstorage/69c3204a686690fea0927d64 , https://files.osf.io/v1/resources/uvqmw/providers/osfstorage/69c3204a686690fea0927d64 , https://files.osf.io/v1/resources/uvqmw/providers/osfstorage/69c3204a686690fea0927d64 , https://osf.io/download/69c3204a686690fea0927d64/ , https://mfr.osf.io/render?url=https%3A%2F%2Fosf.io%2Fdownload%2F69c3204a686690fea0927d64%2F%3Fdirect%26mode%3Drender, https://osf.io/uvqmw/files/osfstorage/69c3204a686690fea0927d64 , https://api.osf.io/v2/files/69c3204a686690fea0927d64/ , https://api.osf.io/v2/files/69c2b365ae03e5ce6d107b09/ , 69c2b365ae03e5ce6d107b09 , files , https://api.osf.io/v2/files/69c3204a686690fea0927d64/versions/ , https://api.osf.io/v2/nodes/uvqmw/ , uvqmw , nodes , https://api.osf.io/v2/nodes/uvqmw/ , nodes , nodes , uvqmw , https://api.osf.io/v2/files/69c3204a686690fea0927d64/cedar_metadata_records/
study_polygons.gpkg 69c2c12bcf7cb14d7b48c456 demo/study_polygons.gpkg study_polygons.gpkg , file , /69c2c12bcf7cb14d7b48c456 , 2670592 , osfstorage , /study_polygons.gpkg , 1774371115 , 1774371115 , 08c72aa055cb38b1473caa3a8063527c , 6734b6f66ebf27962f003e81962bd391cf1381fc74660ddabe64b01fc68b8160 , 11 , FALSE , 1 , FALSE , https://api.osf.io/v2/files/69c2c12bcf7cb14d7b48c456/ , https://files.osf.io/v1/resources/uvqmw/providers/osfstorage/69c2c12bcf7cb14d7b48c456 , https://files.osf.io/v1/resources/uvqmw/providers/osfstorage/69c2c12bcf7cb14d7b48c456 , https://files.osf.io/v1/resources/uvqmw/providers/osfstorage/69c2c12bcf7cb14d7b48c456 , https://osf.io/download/69c2c12bcf7cb14d7b48c456/ , https://mfr.osf.io/render?url=https%3A%2F%2Fosf.io%2Fdownload%2F69c2c12bcf7cb14d7b48c456%2F%3Fdirect%26mode%3Drender, https://osf.io/uvqmw/files/osfstorage/69c2c12bcf7cb14d7b48c456 , https://api.osf.io/v2/files/69c2c12bcf7cb14d7b48c456/ , https://api.osf.io/v2/files/69c2b365ae03e5ce6d107b09/ , 69c2b365ae03e5ce6d107b09 , files , https://api.osf.io/v2/files/69c2c12bcf7cb14d7b48c456/versions/ , https://api.osf.io/v2/nodes/uvqmw/ , uvqmw , nodes , https://api.osf.io/v2/nodes/uvqmw/ , nodes , nodes , uvqmw , https://api.osf.io/v2/files/69c2c12bcf7cb14d7b48c456/cedar_metadata_records/
sim_locations_2026-03-24.csv 69c2f06142e811a7159282bc demo/sim_locations_2026-03-24.csv sim_locations_2026-03-24.csv , file , /69c2f06142e811a7159282bc , 17625680 , osfstorage , /sim_locations_2026-03-24.csv , 1774383201 , 1774383201 , 90d4cf68473c797fb0d4b8b007fa9519 , 852f7fc148ea5051bc70337ac34ad51f5f8c59f536d88a6f2a555ab8c9cf956d , 9 , FALSE , 1 , FALSE , https://api.osf.io/v2/files/69c2f06142e811a7159282bc/ , https://files.osf.io/v1/resources/uvqmw/providers/osfstorage/69c2f06142e811a7159282bc , https://files.osf.io/v1/resources/uvqmw/providers/osfstorage/69c2f06142e811a7159282bc , https://files.osf.io/v1/resources/uvqmw/providers/osfstorage/69c2f06142e811a7159282bc , https://osf.io/download/69c2f06142e811a7159282bc/ , https://mfr.osf.io/render?url=https%3A%2F%2Fosf.io%2Fdownload%2F69c2f06142e811a7159282bc%2F%3Fdirect%26mode%3Drender, https://osf.io/uvqmw/files/osfstorage/69c2f06142e811a7159282bc , https://api.osf.io/v2/files/69c2f06142e811a7159282bc/ , https://api.osf.io/v2/files/69c2b365ae03e5ce6d107b09/ , 69c2b365ae03e5ce6d107b09 , files , https://api.osf.io/v2/files/69c2f06142e811a7159282bc/versions/ , https://api.osf.io/v2/nodes/uvqmw/ , uvqmw , nodes , https://api.osf.io/v2/nodes/uvqmw/ , nodes , nodes , uvqmw , https://api.osf.io/v2/files/69c2f06142e811a7159282bc/cedar_metadata_records/

4.0.1 Unzip Grids

Spatial grids are bundled in zip files. Unzip the files.

Show R Code
osf_id$name # file names
[1] "spread_GeoTIFF_2026-03-24.zip" "niche_GeoTIFF_2026-03-24.zip" 
[3] "study_polygons.gpkg"           "sim_locations_2026-03-24.csv" 
Show R Code
unzip(file.path(here("demo"),"spread_GeoTIFF_2026-03-24.zip"), exdir = here("demo")) # spread grids
unzip(file.path(here("demo"), "niche_GeoTIFF_2026-03-24.zip"), exdir = here("demo")) # potential abundance grids

4.0.2 Spatial Polygons

Read polygons defining jurisdictional and administrative boundaries

Show R Code
states_proj <- st_read(here("demo/study_polygons.gpkg"))
Reading layer `study_states' from data source 
  `D:\Github\hominivorax-geostat\demo\study_polygons.gpkg' using driver `GPKG'
Simple feature collection with 492 features and 122 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -4363.574 ymin: -2203.058 xmax: 2962.346 ymax: 2835.693
Projected CRS: +proj=aea +lat_0=20 +lon_0=-75 +lat_1=10 +lat_2=30 +x_0=0 +y_0=0 +datum=WGS84 +units=km +no_defs

5 Inspect Grids

Spatial grid file names are coded with the year and week they represent. The stack_rasters() function will read the files from the corresponding directory, order them by year and week, and return the data as SpatRaster objects. Values in both objects provide the number of NWS detections per 25x25~km grid cell (approx. 625~km^2). The niche product is the number of potential detections based on environmental conditions, whereas the spread product is an estimate of how many occurred between January 2022 and February 2026.

Organize and stack rasters.

Show R Code
abundance_stk <- stack_rasters(here("demo/niche")) # potential abundance given niche conditions, limits, or capacity
spread_stk <- stack_rasters(here("demo/spread")) # reflecting invasion/spread

5.0.1 Example Potential Abundance Maps

Show R Code
quick_map_plot(grid_stk=abundance_stk, target_year=2024, target_week=10, boundary_polys=states_proj)

Show R Code
quick_map_plot(grid_stk=abundance_stk, target_year=2025, target_week=30, boundary_polys=states_proj)

5.0.2 Example Spread Maps

Show R Code
quick_map_plot(grid_stk=spread_stk, target_year=2024, target_week=10, boundary_polys=states_proj)

Show R Code
quick_map_plot(grid_stk=spread_stk, target_year=2025, target_week=30, boundary_polys=states_proj)

6 Estimated Detections

The sum_intensity_grid() function, when applied to the estimated spread grids, will sum the estimated number of detections for each week. This approximates the mean number of cases as predicted by the model described in the manuscript. Note that the function includes an input option min_est= that specifies the minimum cell value to consider a detection.

The dataframe returned by the function includes a column count that will be used to simulate spatial points in the next step.

Show R Code
weekly_counts <- sum_intensity_grid(spread_stk, min_est = 1)

det_plt <- ggplot(data = weekly_counts, aes(x = date, y = count)) +
  geom_col(fill = "steelblue", alpha = 0.6, width = 7) + 
  theme_minimal() +
  scale_x_date(date_labels = "%b %Y", date_breaks = "3 months") + 
  scale_y_continuous(breaks=scales::pretty_breaks(6)) + 
  theme(
    plot.margin = unit(c(5, 0.01, 5, 0.01), "cm"),
    legend.position = " none",
    legend.background = element_rect(fill = "white", color = "grey90"),
    plot.title = element_text(face = "bold", size = 16),
    axis.title   = element_text(size = 22, face = "bold"),
    axis.text.x    = element_text(size = 16, face = "bold", angle = 45, hjust = 1),
    axis.text.y    = element_text(size = 16, face = "bold")
  ) +
  labs(
    title = " ",
    subtitle = " ",
    x = " ",
    y = "Estimated Case Count",
    caption = " "
  )
Show R Code
det_plt

Figure 1: Estimated detections by week. Counts represent the weekly sum of counts in all raster cells.

7 Simulate Spatial Points

The simulate_detection_points() function will use the spatstat package to generate random spatial points based on target weekly counts given in the weekly_counts dataframe created in the prior step, and estimated intensity as given in the spread grids (spread_stk). Values in spread_stk serve as weights for probabilistic assignment of point locations.

The returned result is a dataframe with x and y coordinates, year, and week. Note that the projected coordinate system is used by default. Conversion to latitude and longitude can be performed as shown below.

Show R Code
sim_points <- simulate_detection_points(raster_stack=spread_stk, # weights for location assignment
                                        obs_df=weekly_counts, # how many to simulate per week
                                        rseed=1976 # change to alter assignment probability
                                        )

n_sim_points <- dim(sim_points)[1] # number simulated
n_est_detctions <- sum(weekly_counts$count) # number estimated by model
n_sim_points == n_est_detctions # same
[1] TRUE
Show R Code
head(sim_points)
x y year week
-536.4614 -1212.174 2022 1
-580.4354 -1296.023 2022 1
-510.0831 -1226.140 2022 1
-374.1382 -1268.373 2022 1
-579.8311 -1248.473 2022 1
-618.1229 -1369.127 2022 1
Show R Code
# convert to SpatVector points
nws_points <- vect(sim_points, geom = c("x", "y"), crs = crs(spread_stk)) # projected 
nws_points <- project(nws_points, "EPSG:4326") # convert to lat/long


# back to dataframe, if needed
nws_point_latlong_df <- as.data.frame(nws_points, geom = "xy")
head(nws_point_latlong_df)
year week x y
2022 1 -79.86693 9.073774
2022 1 -80.24219 8.301078
2022 1 -79.62429 8.953852
2022 1 -78.38461 8.599727
2022 1 -80.24999 8.732931
2022 1 -80.56074 7.625281

Save a copy of simulated point locations. This file will be needed to run the Tessellation step.

Note: A similar file is included in initial download from OSF above.

Show R Code
write_csv(nws_point_latlong_df, here("demo/sim_locations_demo.csv"))

7.0.1 Map Points

Visual inspection of simulated locations for select weeks.

Show R Code
quick_point_map(target_year=2022, target_week=10, 
                boundary_polys=states_proj, 
                point_data=sim_points)

Figure 2: Simulated NWS case detections.

Show R Code
quick_point_map(target_year=2024, target_week=30, 
                boundary_polys=states_proj, 
                point_data=sim_points)

Figure 3: Simulated NWS case detections.

Show R Code
quick_point_map(target_year=2026, target_week=2, 
                boundary_polys=states_proj, 
                point_data=sim_points)

Figure 4: Simulated NWS case detections.

Next Step

Once supporting data has been downloaded and NWS case detections have been simulated, proceed to the Tessellation page of this site.