Source Data Download and NWS Case Detection Simulation
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.
Load needed libraries and packages.
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 manipulationLoad customized functions.
source(here("R/utilities.R"))
source_dir(here("R"))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.
# 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"# 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/ |
Spatial grids are bundled in zip files. Unzip the files.
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" 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 gridsRead polygons defining jurisdictional and administrative boundaries
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_defsSpatial 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.
abundance_stk <- stack_rasters(here("demo/niche")) # potential abundance given niche conditions, limits, or capacity
spread_stk <- stack_rasters(here("demo/spread")) # reflecting invasion/spreadquick_map_plot(grid_stk=abundance_stk, target_year=2024, target_week=10, boundary_polys=states_proj)
quick_map_plot(grid_stk=abundance_stk, target_year=2025, target_week=30, boundary_polys=states_proj)
quick_map_plot(grid_stk=spread_stk, target_year=2024, target_week=10, boundary_polys=states_proj)
quick_map_plot(grid_stk=spread_stk, target_year=2025, target_week=30, boundary_polys=states_proj)
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.
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 = " "
)det_plt
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.
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] TRUEhead(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 |
# 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.
write_csv(nws_point_latlong_df, here("demo/sim_locations_demo.csv"))Visual inspection of simulated locations for select weeks.
quick_point_map(target_year=2022, target_week=10,
boundary_polys=states_proj,
point_data=sim_points)
quick_point_map(target_year=2024, target_week=30,
boundary_polys=states_proj,
point_data=sim_points)
quick_point_map(target_year=2026, target_week=2,
boundary_polys=states_proj,
point_data=sim_points)