Prototype Ocean Indicators
Sarah Gaichas
2026-04-28
1 Introduction
Ocean data is available for the Southeast US coast, from multiple sources. Regional ocean models and global ocean models, including the NOAA MOM6 model [1] (https://psl.noaa.gov/cefi_portal/) and the Copernicus products (https://data.marine.copernicus.eu/products) including the Global Oceanic Reanalysis (GLORYS, [2], https://data.marine.copernicus.eu/product/GLOBAL_MULTIYEAR_PHY_001_030/description) already provide both hindcasts and forecasts of ocean temperature and salinity for the South Atlantic. This project will identify specific resources necessary to translate available data into ecosystem indicators tailored for the South Atlantic and its fishery management needs. For example, the NEFSC has publicly available code to derive indicators from these data sources for the region north of Cape Hatteras that could be modified to derive prototype indicators for the South Atlantic region. Here, I develop prototype ocean indicators to illustrate both opportunities and challenges for operational application in the region.
A prototype habitat indicator characterizing bottom temperature in space over time was developed using GLORYS data with a reproducible workflow based on workflows used in the Northeast US by NEFSC. The objective is to produce maps and time series indicating where and when bottom temperature conditions approach or exceed stressful levels input particular species ranging from corals to targeted fish species. Stressful temperature levels would need to be determined from literature or other studies if unknown.
Validation of this model-based temperature indicator could be achieved with an expanded Citizen Science program; establishment of a formal Study Fleet such as that in the Northeast US [3] will be outlined if the indicator shows promise.
2 Methods
2.1 Pull GLORYS bottom temperature
To use GLORYS information, users must register for a free account on the Copernicus website by clicking the link at the top of the page. Users supply an email and optional information on their organization (one can select Others) and then receive an email with a link to establish a password.
I am working in R [4], but Python code is used to access GLORYS information (there is not a way to access in R). Python can be installed using the R package reticulate [5] and the command reticulate::install_miniconda() if you don’t already have a Python installation. This installs Miniconda, a lightweight distribution of Python. Because it is a minimal distribution, one or two packages need to be added as below.
For this script you need the python pandas package, which is installed using reticulate::py_install("pandas")
These lines are left here for the record but I was unable to use themotuclient package
You also need the motuclient package, from https://github.com/clstoulouse/motu-client-python
Try reticulate::py_install("motuclient", pip = TRUE) which worked. Almost. reticulate::py_install("setuptools", pip = TRUE)
Maybe motuclient doesnt work anymore?
I also didn’t need to install with pip, which was recommended for motuclient but not needed for other packages, I kept using it to stay consistent
The lines below would work, possibly better, without the pip=“TRUE” argument
We also need to install the package from Copernicus that makes downloading the information easy:
Install the copernicusmarine package using reticulate::py_install("copernicusmarine", pip = TRUE)
Install h5py. reticulate::py_install("h5py", pip = TRUE)
The dataset id we want is cmems_mod_glo_phy_my_0.083deg_P1D-m
With the username and password, the following code can be used to download data locally for processing. It is important to specify the latitude and longitude bounds of the data as well as the year(s). Downloading one year at a time makes for manageable files to process for spatial or time series indicators. The output directory is hardcoded below, this needs to change for use elsewhere.
This code for downloading is courtesy Joseph Caracappa, NOAA NMFS NEFSC, and is a test downloading NEUS information so I can compare my results with published indicators to ensure my code is working correctly:
import datetime
import pandas
import os
from os import path
#Set start and end year
year_start=2020
year_end=2024
years=range(year_start,year_end+1)
#Input CMEMS User and Password (case sensitive)
USER = ""
PASSWORD = ""
#Set Lat/lon bounds
min_lon = str(-82.5)
max_lon = str(-51.5)
min_lat = str(22.5)
max_lat = str(48.5)
for y in range(len(years)):
out_dir = "...ADD YOUR PATH HERE.../SAFMCindicators/GLORYS/testNE/"+str(years[y])+"/"
if(not os.path.exists(out_dir)):
os.makedirs(out_dir)
dt = datetime.datetime(years[y],1,1)
end = datetime.datetime(years[y],12,31)
step = datetime.timedelta(days = 1)
all_days = []
while dt <= end:
all_days.append(dt.strftime('%Y-%m-%d'))
dt += step
for d in range(len(all_days)):
t1 = all_days[d]+"T00:00:00"
t2 = all_days[d]+"T23:59:59"
#Change to desired filename prefix
new_name = "GLORYS_REANALYSIS_"+all_days[d]+".nc"
#Change to appropriate output path
if(path.exists(out_dir+new_name)):
print(new_name+" EXISTS")
continue
#print(new_name)
#If additional variables are desired: need to add "--variable var.name"
command = "copernicusmarine subset -i cmems_mod_glo_phy_my_0.083deg_P1D-m -x "+min_lon+" -X "+max_lon+" -y "+min_lat+" -Y "+max_lat+" -z 0. -Z 5000. -t "+t1+" -T "+t2+" -v sea_water_potential_temperature_at_sea_floor -o "+out_dir+" -f "+new_name+" --force-download --username "+USER+" --password "+PASSWORD
#command = "copernicusmarine subset -i cmems_mod_glo_phy_myint_0.083deg_P1D-m -x "+min_lon+" -X "+max_lon+" -y "+min_lat+" -Y "+max_lat+" -z 0. -Z 5000. -t "+t1+" -T "+t2+" -v thetao -o "+out_dir+" -f "+new_name+" --force-download"
#print(command)
os.system(command)There are separate files for each day of each year.
A single file called GLORYS_REANALYSIS_2025-12-31.nc is now downloaded to the local GLORYS/testNE/2025 folder.
Take a look at it
tidync::tidync(here::here("GLORYS/testNE/2025/GLORYS_REANALYSIS_2025-12-31.nc"))##
## Data Source (1): GLORYS_REANALYSIS_2025-12-31.nc ...
##
## Grids (4) <dimension family> : <associated variables>
##
## [1] D2,D1,D0 : bottomT **ACTIVE GRID** ( 116749 values per variable)
## [2] D0 : time
## [3] D1 : latitude
## [4] D2 : longitude
##
## Dimensions 3 (all active):
##
## dim name length min max start count dmin dmax unlim coord_dim
## <chr> <chr> <dbl> <dbl> <dbl> <int> <int> <dbl> <dbl> <lgl> <lgl>
## 1 D0 time 1 6.66e5 6.66e5 1 1 6.66e5 6.66e5 FALSE TRUE
## 2 D1 lati… 313 2.25e1 4.85e1 1 313 2.25e1 4.85e1 FALSE TRUE
## 3 D2 long… 373 -8.25e1 -5.15e1 1 373 -8.25e1 -5.15e1 FALSE TRUEPlot it?
ncfile <- here::here("GLORYS/testNE/2025/GLORYS_REANALYSIS_2025-12-31.nc")
# find the unit from https://ropensci.org/blog/2019/11/05/tidync/
tunit <- ncmeta::nc_atts(ncfile, "time") %>% tidyr::unnest(cols = c(value)) %>% dplyr::filter(name == "units")
# its hours since 1950-01-01
glorysnctotibble <- function(ncfile = ncfile){
# time is hours since 1950-01-01
origin <- as.Date("1950-01-01")
# add the time variable into each
time <- tidync::tidync(ncfile) |>
tidync::activate("D0") |> tidync::hyper_array()
# google says 86400000 milliseconds in a day
# irrelevant because this unit is hours in a day
bttib <- tidync::tidync(ncfile) |>
tidync::hyper_tibble(force = TRUE) |>
dplyr::mutate(time = unname(time$time)) |>
dplyr::mutate(date = as.Date((time/24), origin = origin),
year = lubridate::year(date),
month = lubridate::month(date),
day = lubridate::day(date))|>
dplyr::mutate(longitude = as.numeric(longitude),
latitude = as.numeric(latitude)) |>
dplyr::rename(mod_bt = bottomT)
return(bttib)
}
dailybtplot <- function(oneday){
ggplot() +
geom_tile(data = oneday, aes(x = longitude, y = latitude, fill = mod_bt)) +
geom_sf(data = ecodata::coast) +
#geom_point(data = FishStatsUtils::northwest_atlantic_grid, aes(x = Lon, y = Lat), size=0.05, alpha=0.1) +
scale_fill_gradientn(name = "Temp C",
limits = c(0.5, 31),
colours = c(scales::muted("blue"), "white",
scales::muted("red"), "black")) +
coord_sf(xlim = c(-77, -65), ylim = c(35, 45)) +
ecodata::theme_map() +
ggtitle(paste("Bottom temp, mm-dd-yyyy:", unique(oneday$month),
unique(oneday$day), unique(oneday$year), sep = " "))
}
dec3125 <- glorysnctotibble(ncfile)
dailybtplot(dec3125)
To highlight a box for the South Atlantic region, I visited https://data.marine.copernicus.eu/product/GLOBAL_MULTIYEAR_PHY_001_030/download then selected variables, drew a box around area on the map, selected a year range, then clicked “automate” to generate a command line script. This script was modified and put into the same framework as used above for the Northeast region.
This script pulled initial bottom temperature files for the South Atlantic, 2020-2025:
import datetime
import pandas
import os
from os import path
#Set start and end year
year_start=2020
year_end=2025
years=range(year_start,year_end+1)
#Input CMEMS User and Password (case sensitive)
USER = ""
PASSWORD = ""
#Set Lat/lon bounds
min_lon = str(-82.8)
max_lon = str(-70.6)
min_lat = str(23.4)
max_lat = str(36.1)
for y in range(len(years)):
out_dir = "[....]/GLORYS/SAtl/"+str(years[y])+"/"
if(not os.path.exists(out_dir)):
os.makedirs(out_dir)
dt = datetime.datetime(years[y],1,1)
end = datetime.datetime(years[y],12,31)
step = datetime.timedelta(days = 1)
all_days = []
while dt <= end:
all_days.append(dt.strftime('%Y-%m-%d'))
dt += step
for d in range(len(all_days)):
t1 = all_days[d]+"T00:00:00"
t2 = all_days[d]+"T23:59:59"
#Change to desired filename prefix
new_name = "GLORYS_REANALYSIS_"+all_days[d]+".nc"
#Change to appropriate output path
if(path.exists(out_dir+new_name)):
print(new_name+" EXISTS")
continue
#print(new_name)
#If additional variables are desired: need to add "--variable var.name"
command = "copernicusmarine subset -i cmems_mod_glo_phy_my_0.083deg_P1D-m -x "+min_lon+" -X "+max_lon+" -y "+min_lat+" -Y "+max_lat+" -z 0. -Z 5000. -t "+t1+" -T "+t2+" -v sea_water_potential_temperature_at_sea_floor -o "+out_dir+" -f "+new_name+" --force-download --username "+USER+" --password "+PASSWORD
#command = "copernicusmarine subset -i cmems_mod_glo_phy_myint_0.083deg_P1D-m -x "+min_lon+" -X "+max_lon+" -y "+min_lat+" -Y "+max_lat+" -z 0. -Z 5000. -t "+t1+" -T "+t2+" -v thetao -o "+out_dir+" -f "+new_name+" --force-download"
#print(command)
os.system(command)Pulling 6 years of data took my laptop and home connection about 8 hours.
The larger footprint took closer to 11 hours.
Two dates did not download which makes the later code break: I downloaded them by hand using the web interface, which may cause problems?
Original output of bottom_temp_model_gridded once I got it running for SAtl
“Missing Dates in 2023: 2023-07-13”
“Missing Dates in 2025: 2025-07-09”
Take a look at a South Atlantic file:
tidync::tidync(here::here("GLORYS/SAtl/2025/GLORYS_REANALYSIS_2025-12-31.nc"))##
## Data Source (1): GLORYS_REANALYSIS_2025-12-31.nc ...
##
## Grids (4) <dimension family> : <associated variables>
##
## [1] D2,D1,D0 : bottomT **ACTIVE GRID** ( 22338 values per variable)
## [2] D0 : time
## [3] D1 : latitude
## [4] D2 : longitude
##
## Dimensions 3 (all active):
##
## dim name length min max start count dmin dmax unlim coord_dim
## <chr> <chr> <dbl> <dbl> <dbl> <int> <int> <dbl> <dbl> <lgl> <lgl>
## 1 D0 time 1 6.66e5 6.66e5 1 1 6.66e5 6.66e5 FALSE TRUE
## 2 D1 lati… 153 2.34e1 3.61e1 1 153 2.34e1 3.61e1 FALSE TRUE
## 3 D2 long… 146 -8.28e1 -7.07e1 1 146 -8.28e1 -7.07e1 FALSE TRUECan we plot this file?
# plotting function modified for the South Atlantic just shifts coord_sf
dailybtplotSA <- function(oneday){
ggplot() +
geom_tile(data = oneday, aes(x = longitude, y = latitude, fill = mod_bt)) +
geom_sf(data = ecodata::coast) +
#geom_point(data = FishStatsUtils::northwest_atlantic_grid, aes(x = Lon, y = Lat), size=0.05, alpha=0.1) +
scale_fill_gradientn(name = "Temp C",
limits = c(0.5, 31),
colours = c(scales::muted("blue"), "white",
scales::muted("red"), "black")) +
coord_sf(xlim = c(-82.5, -71), ylim = c(23.5, 36)) +
#ecodata::theme_map() +
ggtitle(paste("Bottom temp, mm-dd-yyyy:", unique(oneday$month),
unique(oneday$day), unique(oneday$year), sep = " "))
}
ncfile <- here::here("GLORYS/SAtl/2025/GLORYS_REANALYSIS_2025-12-31.nc")
dec3125 <- glorysnctotibble(ncfile)
dailybtplotSA(dec3125)
The NE pulls already include most of the South Atlantic, which is why they are bigger files…
2.2 Process GLORYS bottom temperature
Joe Caracappa, NEFSC, wrote nearly all of the code used to produce these indicators for the Northeast US State of the Ecosystem (SOE) reports.
R packages used in this workflow include tidyverse [6] for general data organization and visualization, terra [7] and ncdf4 [8] for processing netCDF files downloaded from Copernicus, and knitr [9–11], rmarkdown [12–14], and bookdown [15,16] for producing reports in both .html and .pdf formats.
Code for reading the file into time series or maps is found at the NEFSC GitHub
https://github.com/NEFSC/READ-EDAB-GLORYS/tree/main
This repo has a fuller version of the pulling script that adds variables like salinity or mixed layer depth and pulls them into variable specific folders. It could be modified simply by changing the min and max lat and lon as above.
We will focus on the spatial indicator of days with bottom temperature over a certain threshold in a given year. This indicator is called thermal_habitat_gridded in the SOE.
https://github.com/NEFSC/READ-EDAB-GLORYS/blob/main/R/make_thermal_habitat_gridded.R
The full set of SOE indicators using GLORYS bottom temperature are made by the script make_GLORYS_bt_soe.R but we don’t want to make all of them right now for this demonstration.
It also has functions to process the downloaded files into datasets for presentation as either time series or spatial datasets. It calls utility functions in this repository:
https://github.com/NEFSC/READ_EDAB_Utilities/tree/main
The Utilities should be installed so the functions can be used by other functions in READ-EDAB-GLORYS using code devtools::install_github("https://github.com/NEFSC/READ_EDAB_Utilities")
This installed properly so we can use all the utilities functions. One function we will use is make_2d_deg_day_gridded
Also lets try installing the EDAB-GLORYS functions:
devtools::install_github("https://github.com/NEFSC/READ-EDAB-GLORYS")
That doesn’t work due to a missing local file EDAB.GLORYS/data/cold_pool/ROMS/roms_debiased_cold_pool_1959_1992.nc This ROMS information is from a debiased Regional Ocean Modeling System (ROMS) model run for the Northeast US that provides bottom temerpature data from the 1950s to 1993, because the GLORYS dataset starts in 1993. We would not have this historical information available in the South Atlantic so we need to proceed without this dependency no matter what.
I cloned the GitHub repo and built locally by commenting out lines 74-76 in make_cold_pool_extent_shp.R, a function we will not need for the South Atlantic.
Now the other EDAB-GLORYS functions should work for me locally. We’ll take the initial processing from make_GLORYS_bt_soe.R and skip to making only the gridded plots.
Finally, functions to plot indicators are found in the ecodata [17] package and in the NOAA-EDAB GitHub ecodata repository, such as plot_thermal_habitat_gridded : https://github.com/NOAA-EDAB/ecodata/blob/pre-production/R/plot_thermal_habitat_gridded.R
First we try producing the NEUS thermal_habitat_gridded indicator to validate the workflow. All of the hardcoded file paths in the original code need to be changed for this anyway. This code is modified from make_GLORYS_bt_soe.R:
input.dir = '/Users/sarahgaichas/Documents/Work/SAFMCindicators/SAFMCindicators/GLORYS/testNE/'
output.dir = '/Users/sarahgaichas/Documents/Work/SAFMCindicators/SAFMCindicators/GLORYS/testNE-out/'
supp.dir = '/Users/sarahgaichas/Documents/Work/SAFMCindicators/READ-EDAB-GLORYS/data-raw/'
#print('Using default arguments')
shp.file = paste0(supp.dir,'geometry/EPU_NOESTUARIES.shp')
#Get year range
year.start = 2020
year.end = 2024 #format(Sys.time(), '%Y') this would go to now, our data end in 2024
run.years = year.start:year.end
y=1
input.prefix = 'GLORYS_REANALYSIS_'
check.dir = function(file){
if(!dir.exists(dirname(file))){dir.create(dirname(file),recursive =T)}
}
#Produce each year's indices
for(y in 1:length(run.years)){
input.year = run.years[y]
input.diryr = paste0(input.dir, input.year,"/")
input.files = list.files(input.diryr,input.prefix)
input.files.year = as.numeric(gsub(".*(\\d{4}).*", '\\1', input.files))
no.threshold =F
which.files.year = which(input.files.year == run.years[y])
if(length(which.files.year) == 0){
next()
}
this.year.files = paste0(input.diryr,input.files[which.files.year])
# if(any(file.size(this.year.files)==0)){
# next()
# }
#
# if(all(!file.exists(this.year.files))){
# print(paste0('File does not exist: ',this.year.file))
# next()
# }
year.dates = seq.Date(as.Date(paste0(run.years[y],'-01-01')), as.Date(paste0(run.years[y],'-12-31')), by = '1 day')
file.dates <- sub(".*(\\d{4}-\\d{2}-\\d{2}).*", "\\1", this.year.files)
if(!all(year.dates %in% file.dates)){
print(paste0('Missing Dates in ',run.years[y],': ', paste0(year.dates[which(!(year.dates %in% file.dates))],collapse = ', ')))
next()
}
#Make bottom_temp_model_gridded
print(paste0('Starting bottom_temp_model_gridded for year: ',run.years[y], ' at ', Sys.time()))
output.file3 = paste0(output.dir,'data/bottom_temp_model_gridded/GLORYS_bottom_temp_model_gridded_',run.years[y],'.csv')
check.dir(output.file3)
if(file.exists(output.file3)){
print(paste0('File already exists: ',output.file3))
}else{
EDAB.GLORYS::make_bottom_temp_model_gridded(input.file = this.year.files,
output.file = output.file3,
shp.file = shp.file,
file.year = run.years[y],
write.out = T)
}
# Missing output.file1 seasonal anomaly
# #Make bottom_temp_model_anom
# print(paste0('Starting bottom_temp_model_anom for year: ',run.years[y], ' at ', Sys.time()))
#
# output.file4 = paste0(output.dir,'data/bottom_temp_model_anom/GLORYS_bottom_temp_model_anom_',run.years[y],'.csv')
# check.dir(output.file4)
# if(file.exists(output.file4)){
# print(paste0('File already exists: ',output.file4))
# }else{
# EDAB.GLORYS::make_bottom_temp_model_anom(input.file = this.year.files,
# output.file = output.file4,
# shp.file =shp.file,
# file.year = run.years[y],
# climatology.file = output.file1,
# write.out =T)
# }
#Make bottom_temp_model_annual
# add this one already in the repo
output.file2 = paste0(supp.dir, "GLORYS_bottom_temp_clim_1990_2020.csv")
print(paste0('Starting bottom_temp_model_annual for year: ',run.years[y], ' at ', Sys.time()))
output.file5 = paste0(output.dir,'data/bottom_temp_model_annual/GLORYS_bottom_temp_model_annual_',run.years[y],'.csv')
check.dir(output.file5)
if(file.exists(output.file5)){
print(paste0('File already exists: ',output.file5))
}else{
EDAB.GLORYS::make_bottom_temp_model_annual(input.file = this.year.files,
output.file = output.file5,
shp.file = shp.file,
file.year = run.years[y],
climatology.file = output.file2,
write.out =T
)
}
# print(paste0('Starting bottom_temp_daily_epu for year: ',run.years[y], ' at ', Sys.time()))
#
# output.file5b = paste0(output.dir,'data/bottom_temp_daily_epu/GLORYS_bottom_temp_daily_epu_',run.years[y],'.csv')
# check.dir(output.file5b)
# if(file.exists(output.file5b)){
# print(paste0('File already exists: ',output.file5b))
# }else{
# EDAB.GLORYS::make_bottom_temp_daily_epu(input.file = this.year.files,
# output.file = output.file5b,
# shp.file = shp.file,
# file.year = run.years[y],
# write.out =T
# )
#
# }
#Make thermal_habitat_gridded
print(paste0('Starting thermal_habitat_gridded for year: ',run.years[y], ' at ', Sys.time()))
output.file6 = paste0(output.dir,'data/thermal_habitat_gridded/GLORYS_thermal_habitat_gridded_',run.years[y],'.nc')
check.dir(output.file6)
if(file.exists(output.file6) | no.threshold){
print(paste0('File already exists: ',output.file6))
}else{
EDAB.GLORYS::make_thermal_habitat_gridded(input.file = this.year.files,
output.file = output.file6,
supp.dir = supp.dir,
shp.file = shp.file,
file.year = run.years[y],
write.out =T,
t.max.seq = seq(0,30,1)
)
}
#Make thermal_habitat_area
print(paste0('Starting thermal_habitat_area for year: ',run.years[y], ' at ', Sys.time()))
output.file7a = paste0(output.dir,'data/thermal_habitat_area/GLORYS_thermal_habitat_area_',run.years[y],'.csv')
output.file7b = paste0(output.dir,'data/thermal_habitat_gridded/GLORYS_thermal_habitat_gridded_',run.years[y],'.csv')
check.dir(output.file7a)
if(file.exists(output.file7a) & file.exists(output.file7b) | no.threshold){
print(paste0('File already exists: ',output.file7a, ' and ', output.file7b))
}else{
if(any(file.size(this.year.files)==0))
EDAB.GLORYS::make_thermal_habitat_area(input.file = this.year.files,
output.file.area = output.file7a,
output.file.gridded = output.file7b ,
shp.file = shp.file,
file.year = run.years[y],
write.area =T,
write.gridded =T,
t.max.seq = seq(0,30,1)
)
}
}
# glorys.heatwave.files = c(list.files(paste0(output.dir,'data/bottom_temp_daily_epu/'),full.names = T),'/home/jcaracappa/EDAB_Dev/jcaracappa/ROMS_NWA/ROMS_daily_epu_1959_1992.csv' )
# glorys.heatwave.data = lapply(glorys.heatwave.files, read.csv) |>
# dplyr::bind_rows() |>
# dplyr::mutate(source_m = dplyr::coalesce(Source, source)) |>
# dplyr::select(-Source,-source) |>
# dplyr::rename(source = 'source_m')
# #check all dates
# all(1959:2025 %in% sort(unique(format(as.Date(glorys.heatwave.data$date),format = '%Y')) ))
# glorys.heatwaves.out = write.csv(glorys.heatwave.data,paste0('/home/jcaracappa/EDAB_Dev/jcaracappa/ROMS_GLORYS_bottom_temp_model_daily_epu_',1959,'_',format(Sys.time(), '%Y'),'.csv'), row.names =F)
#
# season.index = data.frame(month = 1:12, season = rep(c('Winter','Spring','Summer','Fall'),each = 3))
# data = glorys.heatwave.data %>%
# dplyr::mutate(month = as.numeric(format(as.Date(date),format = '%m')),
# year = as.numeric(format(as.Date(date),format = '%Y'))) %>%
# dplyr::left_join(season.index) %>%
# dplyr::group_by(source,year,season,EPU) %>%
# dplyr::summarise(Value = mean(BottomT.mean,na.rm=T))
#
# library(ggplot2)
# ggplot(data, aes(x = year, y = Value, color = source))+
# geom_line()+
# facet_grid(EPU~season)
# ggsave(here::here('GLORYS_bottom_temp_2025.png'),width =12, height = 12)Did this work for bottom temperature gridded? No only getting winter something is wrong with the dates.
I had to change the EDAB.GLORYS::make_bottom_temp_model_gridded function code to call EDABUtilities::make_2d_summary_gridded with an additional argument, file.time = "daily" to make it work with the files in their original downloaded format. Then running the code above works for bottom temperature grided. Alternatively, all files could be combined into one folder but I am unsure what that looked like so I am modifying the function to have all the steps here.
A similar change was made in EDAB.GLORYS::make_bottom_temp_model_annual to change the call to EDABUtilities::make_2d_summary_ts to `file.time = “daily” instead of annual.
Now I am getting the plots I expected for seasonal gridded bottom temperature.
plot_gridded_bt_NE <- function(bottom_temp_model_gridded_csv, year = NULL, scale = "celsius"){
bottom_temp_model_gridded <- read.csv(bottom_temp_model_gridded_csv) |>
tibble::as_tibble()
xmin = -77
xmax = -65
ymin = 36
ymax = 45
xlims <- c(xmin, xmax)
ylims <- c(ymin, ymax)
if (is.null(year)) year <- max(bottom_temp_model_gridded$Time)
fix <- dplyr::mutate(dplyr::select(dplyr::filter(bottom_temp_model_gridded,
Time == year), -Time), Var = factor(Var, levels = c("winter",
"spring", "summer", "fall")))
if (scale == "fahrenheit") {
fix <- dplyr::mutate(fix, Value = (9/5) * Value + 32)
label <- "Temp. (°F)"
breaks <- c(32, 41, 50, 59, 68, 77)
labelLegend <- c("32", "41", "50", "59", "68", "77")
limits <- c(31.6, 84.2)
midpoint <- 50
} else {
label <- "Temp. (°C)"
breaks <- c(0, 5, 10, 15, 20, 25)
labelLegend <- c("0", "5", "10", "15", "20", "25")
limits <- c(-0.2, 29)
midpoint <- 10
}
p <- ggplot2::ggplot(data = fix) +
ggplot2::geom_tile(ggplot2::aes(x = Longitude,
y = Latitude, fill = sqrt(Value + 1))) +
ggplot2::geom_sf(data = ecodata::coast,
size = 0.4) +
#ggplot2::geom_sf(data = ne_epu_sf, fill = "transparent", size = 0.4) +
ggplot2::coord_sf(xlim = xlims,
ylim = ylims) +
ggplot2::facet_wrap(Var ~.) +
ecodata::theme_map() +
ggplot2::scale_fill_gradient2(name = label, low = "#000004FF",
mid = "#BB3754FF", high = "#FCFFA4FF",
breaks = sqrt(breaks + 1), limits = sqrt(limits + 1),
labels = labelLegend,
midpoint = sqrt(midpoint + 1)) +
# ggplot2::scale_fill_gradient(name = "Temp C",
# limits = c(0.5, 31),
# colours = c(scales::muted("blue"), "white",
# scales::muted("red"), "black")) +
ggplot2::ggtitle(paste("Seasonal Mean Bottom Temperature", year)) +
ggplot2::xlab("Longitude") + ggplot2::ylab("Latitude") +
ggplot2::theme(panel.border = ggplot2::element_rect(colour = "black",
fill = NA, linewidth = 0.75),
legend.key = ggplot2::element_blank(),
axis.title = ggplot2::element_text(size = 11),
strip.background = ggplot2::element_blank(),
strip.text = ggplot2::element_text(hjust = 0),
axis.text = ggplot2::element_text(size = 8),
axis.title.y = ggplot2::element_text(angle = 90)) +
ecodata::theme_title() +
ecodata::theme_ts()
return(p)
}
csv <- here::here("GLORYS/testNE-out/data/bottom_temp_model_gridded/GLORYS_bottom_temp_model_gridded_2024.csv")
plot_gridded_bt_NE(bottom_temp_model_gridded_csv = csv, scale = "celsius")
csv <- here::here("GLORYS/testNE-out/data/bottom_temp_model_gridded/GLORYS_bottom_temp_model_gridded_2023.csv")
plot_gridded_bt_NE(bottom_temp_model_gridded_csv = csv, scale = "celsius")
csv <- here::here("GLORYS/testNE-out/data/bottom_temp_model_gridded/GLORYS_bottom_temp_model_gridded_2022.csv")
plot_gridded_bt_NE(bottom_temp_model_gridded_csv = csv, scale = "celsius")
Produced output for thermal habitat area, but it looks wrong, only has 1 or NA as values when it should be number of days over that temperature; maybe a daily file issue again? I think so, dev branch has a function to put all files in one raster to process.
Diagnose here line by line in make_thermal_habitat_gridded
#from the function call above
# try the dev branch make_thermal_habitat_area which appears to do gridded too
output.file6 = paste0(output.dir,'data/thermal_habitat_gridded/GLORYS_thermal_habitat_gridded_',run.years[y],'.nc')
output.file7a = paste0(output.dir,'data/thermal_habitat_area/GLORYS_thermal_habitat_area_',run.years[y],'.csv')
output.file7b = paste0(output.dir,'data/thermal_habitat_gridded/GLORYS_thermal_habitat_gridded_',run.years[y],'.csv')
check.dir(output.file7a)
input.file = this.year.files
output.file = output.file6
supp.dir = supp.dir
shp.file = shp.file
file.year = run.years[y]
write.out =T
t.max.seq = seq(0,30,1)
output.file.area = output.file7a
output.file.gridded = output.file7b
# make_thermal_habitat_gridded = function(input.file, output.file,supp.dir, shp.file, file.year,t.max.seq, write.out =F){
make_thermal_habitat_area = function(input.file, output.file.area = NA, output.file.gridded =NA, shp.file, file.year,t.max.seq, write.area = F, write.gridded = F){
#EPU.names = c('MAB','GB','GOM','SS')
#EPU.names = c('SAtl') # modify for single large South Atlantic area
#EPU.names = c('Off NC', "Off SC", "Off GA", "Off FL")
#depth.df = data.frame(id = 1:4,
EPU.names = c('MAB','GB','GOM','SS','all')
depth.df = data.frame(depth.min = c(0,25,100, 0),
depth.max = c(25,100,300, 2000),
depth.name = c('0-25m','25-100m','100-300m', 'AllDepths'))
combs = expand.grid(year = file.year,t.max = t.max.seq,depth.name = depth.df$depth.name,EPU = EPU.names,stringsAsFactors = F)%>%
dplyr::left_join(depth.df)
bathy.shp = terra::rast(paste0(supp.dir,'GLORYS/GLORYS_bathymetry_east_coast_crop.nc'),subds = 'deptho')
# --- Refactor Start: Handle Input Files ---
# Convert input to a single SpatRaster.
# terra::rast() handles both a single multi-layer file OR a vector of single-layer files.
processed_rast = terra::rast(input.file)
# Ensure time dimension is set if missing (common when stacking individual daily files)
# This logic matches layers to days in the specified file.year
if (!all(terra::has.time(processed_rast))) {
# Generate daily sequence for the given year
dates = seq(as.Date(paste0(file.year, "-01-01")), as.Date(paste0(file.year, "-12-31")), by="day")
# If layer count matches the days in the year, assign the time
if (terra::nlyr(processed_rast) == length(dates)) {
terra::time(processed_rast) = dates
}
}
# --- Refactor End ---
out.area.ls = list()
out.gridded.ls = list()
for(i in 1:nrow(combs)){
if(combs$EPU[i] == 'all'){
area.names = c('MAB','GB','GOM','SS')
}else{
area.names = combs$EPU[i]
}
if(i ==1){
# Refactor: use processed_rast instead of loading from file again
neus.shp = terra::crop(bathy.shp, processed_rast)
}
depth.rast = terra::clamp(neus.shp, lower = combs$depth.min[i], upper = combs$depth.max[i],values =F)
EPU.vect = terra::vect(shp.file)
area.vect = EPU.vect[which(EPU.vect$EPU %in% area.names)]
area.mask = terra::mask(depth.rast,area.vect)
#Mask of area over t.max
# Refactor: Pass the SpatRaster object (processed_rast) instead of filename
area.i = EDABUtilities::mask_nc_2d(data.in = processed_rast,
write.out =F,
shp.file = area.mask,
var.name = 'BottomT',
min.value = combs$t.max[i],
max.value = Inf,
binary = F,
area.names =NA
)
# Refactor: Subset from the main raster object
this.rast = terra::subset(processed_rast, 1)
nd.i = EDABUtilities::make_2d_deg_day_gridded_nc(data.in = area.i,
shp.file = area.mask,
var.name = 'BottomT',
statistic = 'nd',
type = 'above',
ref.value = combs$t.max[i],
area.names = NA
)
shp.area = terra::expanse(area.mask)$area
area.df = terra::expanse(area.i[[1]]) %>%
as.data.frame()%>%
dplyr::mutate(Time = terra::time(area.i[[1]]),
EPU = combs$EPU[i],
Depth = combs$depth.name[i],
Var = paste0('>',combs$t.max[i],'\u00B0C'),
Value = area/ shp.area,
Source = 'GLORYS',
year = combs$year[i],
temp.threshold = combs$t.max[i],
Units = 'Proportion'
)
out.area.ls[[i]] = area.df
out.gridded.ls[[i]] = as.data.frame(nd.i[[1]],cells =T, xy = T) %>%
dplyr::mutate(Time = combs$year[i], EPU = combs$EPU[i], Depth = combs$depth.name[i], Var = combs$t.max[i], Source = 'GLORYS',Units = 'Number of Days')%>%
dplyr::rename(Latitude = 'y', Longitude = 'x', Value = 'sum')%>%
dplyr::select(Time,EPU, Depth, Var,Value,Latitude,Longitude,Source,Units)
print(signif(i/nrow(combs)*100,2))
}
out.area.df = dplyr::bind_rows(out.area.ls)%>%
dplyr::select(Time, EPU, Depth, Var, Value, Source, year, temp.threshold, Units)%>%
dplyr::mutate(Year = format(as.Date(Time),format = '%Y'))%>%
dplyr::group_by(Year,EPU, Depth, Var, temp.threshold, Units,Source)%>%
dplyr::summarise(Value = mean(Value))%>%
dplyr::rename(Time = Year)
out.gridded.df = dplyr::bind_rows(out.gridded.ls)
if(write.area == T){
write.csv(out.area.df, output.file.area,row.names = F)
}
if(write.gridded == T){
write.csv(out.gridded.df,output.file.gridded,row.names = F)
}
if(write.area == F |write.gridded == F){
return(list(thermal.area = out.area.df,thermal.gridded = out.gridded.df))
}
}
# #create EPU mask layer
# EPU.vect = terra::vect(shp.file)
# neus.shp = terra::crop(bathy.shp,bathy.shp)
#
# epu.ls = lapply(1:length(EPU.names), function(x){
# e.vect = EPU.vect[which(EPU.vect$EPU %in% EPU.names[x])]
# e = terra::mask(neus.shp,e.vect)
# terra::values(e)[!is.na(terra::values(e))] = x
# return(e)
# })
# epu.mask = terra::merge(terra::sprc(epu.ls))
# # epu.mask = as.factor(epu.mask)
# levels(epu.mask) = data.frame(1:4, EPU.names)
# epu.mask.binary = epu.mask
# terra::values(epu.mask.binary)[!is.na(terra::values(epu.mask.binary))] = 1
#
# i=1
# y=32
#
# t=i=1
# out.ind =1
# out.var.names = character()
# year.out.ls = list()
# out.df.ls = list()
# out.ls = list()
# depth.rast.ls = list()
# epu.rast.ls = list()
#
# for(t in 1:length(t.max.seq)){
#
# #Mask of area over t.max
# area.i = EDABUtilities::mask_nc_2d(data.in = input.file,
# write.out =F,
# shp.file = EPU.vect,
# var.name = 'BottomT',
# min.value = t.max.seq[t],
# max.value = Inf,
# binary = F,
# area.names =NA
# )
#
# out.ls[[t]] = EDABUtilities::make_2d_deg_day_gridded_nc(data.in = area.i,
# shp.file = EPU.vect,
# var.name = 'BottomT',
# statistic = 'nd',
# type = 'above',
# ref.value = t.max.seq[t],
# area.names = NA
# )#[[1]]
#
# out.ls[[t]] = out.ls[[t]] *terra::crop(epu.mask.binary, out.ls[[t]])
#
# year.time = as.POSIXct(paste0(file.year, '-01-01 00:00:00UTC'),origin = '1970-01-01 00:00:00',tz = 'UTC')
# terra::time(out.ls[[t]]) = as.numeric(year.time)
#
# print(signif(t/length(t.max.seq) ,2))
#
# out.var.names[t] = paste0('Number of Days per Year Above ',t.max.seq[t],' degrees C')
#
# }
#
# out.sds =terra::sds(out.ls)
#
# out.var = out.sds
#
# #concatenate names
# out.names = c(paste0('nday_',t.max.seq))
# out.names = gsub('\\.','_',out.names)
# names(out.var) = out.names
# terra::longnames(out.var) = out.var.names
# terra::units(out.var) = c(rep('n days',length(t.max.seq)))
#
# terra::writeCDF(out.var,filename = output.file,overwrite =T,missval = 0)
#
# #Format netcdf
# var.atts = read.csv(paste0(supp.dir,'GLORYS/thermal_habitat_gridded_variable_attributes.csv')) %>%
# filter(!is.na(Value) & Attribute.Name != '_FillValue')
# global.atts = read.csv(paste0(supp.dir,'GLORYS/thermal_habitat_gridded_global_attributes.csv'))%>%
# filter(!is.na(Value))
#
# file.nc = ncdf4::nc_open(output.file,write =T)
# file.var = names(file.nc$var)
#
# for(g in 1:nrow(global.atts)){
#
# ncdf4::ncatt_put(file.nc,0,global.atts$Attribute.Name[g],global.atts$Value[g])
#
# }
# ncdf4::ncatt_put(file.nc,0,'time_coverage_start',paste0(file.year,'-01-01'))
# ncdf4::ncatt_put(file.nc,0,'time_coverage_end',paste0(file.year,'-12-31'))
# ncdf4::ncatt_put(file.nc,'time','units','seconds since 1970-01-01T00:00:00Z')
# ncdf4::ncatt_put(file.nc,'latitude','standard_name','latitude')
# ncdf4::ncatt_put(file.nc,'latitude','coverage_content_type','coordinates')
#
# ncdf4::ncatt_put(file.nc,'longitude','standard_name','longitude')
# ncdf4::ncatt_put(file.nc,'longitude','coverage_content_type','coordinates')
#
# for(v in 1:length(file.var)){
#
#
# ncdf4::ncatt_put(file.nc,file.var[v],'standard_name','number_of_days_with_bottom_temperature_above_threshold',prec ='text' )
#
# for(va in 1:nrow(var.atts)){
#
# ncdf4::ncatt_put(file.nc,file.var[v],var.atts$Attribute.Name[va],var.atts$Value[va],prec ='text' )
#
# }
#
# }
#
# ncdf4::nc_close(file.nc)
#
# ##Write into test
# # file.names = list.files(here::here('data','SOE','thermal_habitat_gridded_V2'),'thermal_', full.names = T)
# # test.nc = ncdf4::nc_open(file.names[32])
# # # test.nc = ncdf4::nc_open('C:/Users/joseph.caracappa/Downloads/thermal_habitat_gridded_1993_SOE2025.nc')
# # ncdf4::ncatt_get(test.nc,'nday_10')
# # ncdf4::nc_close(test.nc)
# # x = terra::rast(file.names[1])
# # plot(x)
# # test.sds = terra::sds(file.names[1],1)
# # writeCDF(test.sds,here::here('data','SOE','thermal_habitat_checker.nc'))
#
# }Issue appears to be in EDAButilities::make_2d_deg_day_gridded_nc and I think it is because these are daily files
data.in = area.i
shp.file = EPU.vect
var.name = 'BottomT'
statistic = 'nd'
type = 'above'
ref.value = t.max.seq[t]
area.names = NA
write.out = F
make_2d_deg_day_gridded_nc <- function(data.in,write.out = F,output.files,shp.file,var.name,statistic,ref.value,type,area.names){
if(class(shp.file) %in% c('SpatVector','SpatRaster')){
shp.vect = shp.file
use.shp =T
}else if(!is.na(shp.file)){
shp.vect = terra::vect(shp.file)
use.shp =T
}else{
use.shp = F
}
if(all(!is.na(area.names))){
shp.str = as.data.frame(shp.vect)
which.att = which(apply(shp.str,2,function(x) all(area.names %in% x)))
which.area = match(area.names,shp.str[,which.att])
shp.vect = shp.vect[which.area]
}
out.ls = list()
for(i in 1:length(data.in)){
if(is.character(data.in)){
data = terra::rast(data.in[i])
}else if(class(data.in[[i]])[1] == 'SpatRaster'){
data = data.in[[i]]
}else{
stop('data.in needs to be either file names or spatRasters')
}
data = EDABUtilities::convert_longitude(data)
if(use.shp){
data = terra::mask(data,shp.vect)
}
data.mask = terra::subset(data,1) * 0
if(type == 'raw'){
data.stat = sum(data,na.rm=T)
}else if (type == 'above'){
if(statistic == 'dd'){
data.temp = terra::clamp(data,lower = ref.value, upper = Inf,value = F)
data.stat = sum(data.temp,na.rm=T)
}else if (statistic == 'nd'){
data.temp = (terra::clamp(data,lower = ref.value, upper = Inf,value =F)*0)+1
data.stat = sum(data.temp,na.rm=T)
}else if(statistic == 'nd.con'){
data.temp = (terra::clamp(data,lower = ref.value, upper = Inf,value =F)*0)+1
data.stat = terra::app(data.temp,fun = function(x){
l = rle(x)
m = l$lengths[which(l$values == 1)]
return(ifelse(length(m) == 0, 0,max(m,na.rm=T)))
})
}else{
warning('statistic needs to be "dd" or "nd"')
}
}else if (type == 'below'){
if(statistic == 'dd'){
data.temp = terra::clamp(data,lower = -Inf, upper = ref.value, value =F)
data.stat = sum(data.temp,na.rm=T)
}else if (statistic == 'nd'){
data.temp = (terra::clamp(data,lower = -Inf, upper = ref.value, value =F)*0)+1
data.stat = sum(data.temp,na.rm=T)
}else if(statistic == 'nd.con'){
data.temp = (terra::clamp(data,lower = -Inf, upper = ref.value, value =F)*0)+1
data.stat = terra::app(data.temp,fun = function(x){
l = rle(x)
m = l$lengths[which(l$values == 1)]
return(ifelse(length(m) == 0, 0,max(m,na.rm=T)))
})
}else{
warning('statistic needs to be "dd" or "nd"')
}
}
data.out = sum(data.stat,data.mask,na.rm=T)
if(write.out){
terra::writeCDF(data.out, output.files[i],varname = paste0(var.name,'_',type,'_',ref.value,'_',statistic),overwrite =T)
}else{
out.ls[[i]] = data.out
}
}
if(write.out ==F){
return(out.ls)
}
}Got it, use the dev branch EDAB.GLORYS::make_thermal_habitat_area function
Need to know how it gets from nc to an R dataset, perhaps a clue is in the test code for make_thermal_habitat_gridded
This is NE 2020 testing the dev branch function, which appears to work
file.names = list.files(here::here("GLORYS/testNE-out/data/thermal_habitat_gridded"), full.names = T)
# test.nc = ncdf4::nc_open(file.names[1])
# # # test.nc = ncdf4::nc_open('C:/Users/joseph.caracappa/Downloads/thermal_habitat_gridded_1993_SOE2025.nc')
# ncdf4::ncatt_get(test.nc,'nday_10')
# ncdf4::ncatt_get(test.nc,'nday_5')
# ncdf4::nc_close(test.nc)
# #x = terra::rast(file.names[1])
# #x = terra::rast(test.nc,'nday_5')
# #plot(x)
# # test.sds = terra::sds(file.names[1],1)
# # writeCDF(test.sds,here::here('data','SOE','thermal_habitat_checker.nc'))
#
# # convert to tibble and plot as above modifying ecodata::plot_thermal_habitat_gridded code
# ncfile <- file.names[1] # 2020, retrieve the year from the filename to add to the plot
#
# bttib <- tidync::tidync(ncfile) |>
# tidync::hyper_tibble(force = TRUE) |>
# dplyr::select(nday_5, nday_10, nday_15, nday_20, longitude, latitude) |>
# dplyr::mutate(longitude = as.numeric(longitude),
# latitude = as.numeric(latitude))
btthermalgrid <- read.csv(file.names[1])
fix <- btthermalgrid |>
dplyr::filter(Var %in% c(5, 10, 15, 20))
xmin = -77
xmax = -65
ymin = 36
ymax = 45
xlims <- c(xmin, xmax)
ylims <- c(ymin, ymax)
p <- ggplot2::ggplot(fix) +
ggplot2::geom_tile(ggplot2::aes(x = Longitude,
y = Latitude, color = Value, width = 1/12, height = 1/12),
linewidth = 2) +
ggplot2::geom_sf(data = ecodata::coast,
size = 0.4) +
ggplot2::facet_grid(Depth ~ Var) +
ggplot2::scale_color_viridis_c() +
ggplot2::coord_sf(xlim = c(xmin, xmax), ylim = c(ymin, ymax)) +
ggplot2::xlab("") +
ggplot2::ylab("") +
ecodata::theme_ts() +
ecodata::theme_facet() +
ecodata::theme_title() +
#ggplot2::ggtitle((unique(Time))) +
ggplot2::theme(legend.position = "bottom")
#return(p)
p
Next we try the same code on the South Atlantic GLORYS bottom temperature.
Problem: I need a shape file to run this. Lets try SA_EEZ_off_states from https://www.fisheries.noaa.gov/resource/map/defined-fishery-management-areas-south-atlantic-states-map-gis-data
Which appears too wide for the extent I pulled so lets try pulling again with a bigger area (changed extent in code above).
Next problem: EDAB-GLORYS functions have the Northeast EPUs hardcoded in them, will need to rewrite them. If I get rid of the EPUs and try to use the full South Atlantic area in that shapefile maybe that will work.
First step was changing the code naming EPUs to name only one “SAtl”, see if that works–no there are named regions in the shapefile so use those. Off NC Off SC Off GA and Off FL?
The bottom temp model gridded runs but the thermal habitat fails with Error: [*] extents do not match in the first iteration for 2020
input.dir = '/Users/sarahgaichas/Documents/Work/SAFMCindicators/SAFMCindicators/GLORYS/SAtl/'
output.dir = '/Users/sarahgaichas/Documents/Work/SAFMCindicators/SAFMCindicators/GLORYS/SAtl-out/'
supp.dir = '/Users/sarahgaichas/Documents/Work/SAFMCindicators/READ-EDAB-GLORYS/data-raw/'
#print('Using default arguments')
shp.file = paste0(supp.dir,'geometry/SA_EEZ_off_states.shp')
#shp.file = NA
#Get year range
year.start = 2020
year.end = 2025 #format(Sys.time(), '%Y') this would go to now, our data end in 2024
run.years = year.start:year.end
y=1
input.prefix = 'GLORYS_REANALYSIS_'
check.dir = function(file){
if(!dir.exists(dirname(file))){dir.create(dirname(file),recursive =T)}
}
#Produce each year's indices
for(y in 1:length(run.years)){
input.year = run.years[y]
input.diryr = paste0(input.dir, input.year,"/")
input.files = list.files(input.diryr,input.prefix)
input.files.year = as.numeric(gsub(".*(\\d{4}).*", '\\1', input.files))
no.threshold =F
which.files.year = which(input.files.year == run.years[y])
if(length(which.files.year) == 0){
next()
}
this.year.files = paste0(input.diryr,input.files[which.files.year])
# if(any(file.size(this.year.files)==0)){
# next()
# }
#
# if(all(!file.exists(this.year.files))){
# print(paste0('File does not exist: ',this.year.file))
# next()
# }
year.dates = seq.Date(as.Date(paste0(run.years[y],'-01-01')), as.Date(paste0(run.years[y],'-12-31')), by = '1 day')
file.dates <- sub(".*(\\d{4}-\\d{2}-\\d{2}).*", "\\1", this.year.files)
if(!all(year.dates %in% file.dates)){
print(paste0('Missing Dates in ',run.years[y],': ', paste0(year.dates[which(!(year.dates %in% file.dates))],collapse = ', ')))
next()
}
#Make bottom_temp_model_gridded
print(paste0('Starting bottom_temp_model_gridded for year: ',run.years[y], ' at ', Sys.time()))
output.file3 = paste0(output.dir,'data/bottom_temp_model_gridded/GLORYS_bottom_temp_model_gridded_',run.years[y],'.csv')
check.dir(output.file3)
if(file.exists(output.file3)){
print(paste0('File already exists: ',output.file3))
}else{
EDAB.GLORYS::make_bottom_temp_model_gridded(input.file = this.year.files,
output.file = output.file3,
shp.file = shp.file,
file.year = run.years[y],
write.out = T)
}
# Missing output.file1 seasonal anomaly
# #Make bottom_temp_model_anom
# print(paste0('Starting bottom_temp_model_anom for year: ',run.years[y], ' at ', Sys.time()))
#
# output.file4 = paste0(output.dir,'data/bottom_temp_model_anom/GLORYS_bottom_temp_model_anom_',run.years[y],'.csv')
# check.dir(output.file4)
# if(file.exists(output.file4)){
# print(paste0('File already exists: ',output.file4))
# }else{
# EDAB.GLORYS::make_bottom_temp_model_anom(input.file = this.year.files,
# output.file = output.file4,
# shp.file =shp.file,
# file.year = run.years[y],
# climatology.file = output.file1,
# write.out =T)
# }
# #Make bottom_temp_model_annual
# # add this one already in the repo
# output.file2 = paste0(supp.dir, "GLORYS_bottom_temp_clim_1990_2020.csv")
# print(paste0('Starting bottom_temp_model_annual for year: ',run.years[y], ' at ', Sys.time()))
#
# output.file5 = paste0(output.dir,'data/bottom_temp_model_annual/GLORYS_bottom_temp_model_annual_',run.years[y],'.csv')
# check.dir(output.file5)
# if(file.exists(output.file5)){
# print(paste0('File already exists: ',output.file5))
# }else{
# EDAB.GLORYS::make_bottom_temp_model_annual(input.file = this.year.files,
# output.file = output.file5,
# shp.file = shp.file,
# file.year = run.years[y],
# climatology.file = output.file2,
# write.out =T
# )
#
# }
# print(paste0('Starting bottom_temp_daily_epu for year: ',run.years[y], ' at ', Sys.time()))
#
# output.file5b = paste0(output.dir,'data/bottom_temp_daily_epu/GLORYS_bottom_temp_daily_epu_',run.years[y],'.csv')
# check.dir(output.file5b)
# if(file.exists(output.file5b)){
# print(paste0('File already exists: ',output.file5b))
# }else{
# EDAB.GLORYS::make_bottom_temp_daily_epu(input.file = this.year.files,
# output.file = output.file5b,
# shp.file = shp.file,
# file.year = run.years[y],
# write.out =T
# )
#
# }
# #Make thermal_habitat_gridded
# print(paste0('Starting thermal_habitat_gridded for year: ',run.years[y], ' at ', Sys.time()))
#
# output.file6 = paste0(output.dir,'data/thermal_habitat_gridded/GLORYS_thermal_habitat_gridded_',run.years[y],'.nc')
# check.dir(output.file6)
# if(file.exists(output.file6) | no.threshold){
# print(paste0('File already exists: ',output.file6))
# }else{
# EDAB.GLORYS::make_thermal_habitat_gridded(input.file = this.year.files,
# output.file = output.file6,
# supp.dir = supp.dir,
# shp.file = shp.file,
# file.year = run.years[y],
# write.out =T,
# t.max.seq = seq(0,30,1)
# )
# }
#
# #Make thermal_habitat_area
# print(paste0('Starting thermal_habitat_area for year: ',run.years[y], ' at ', Sys.time()))
#
# output.file7a = paste0(output.dir,'data/thermal_habitat_area/GLORYS_thermal_habitat_area_',run.years[y],'.csv')
# output.file7b = paste0(output.dir,'data/thermal_habitat_gridded/GLORYS_thermal_habitat_gridded_',run.years[y],'.csv')
# check.dir(output.file7a)
#
# if(file.exists(output.file7a) & file.exists(output.file7b) | no.threshold){
# print(paste0('File already exists: ',output.file7a, ' and ', output.file7b))
# }else{
#
# if(any(file.size(this.year.files)==0))
# EDAB.GLORYS::make_thermal_habitat_area(input.file = this.year.files,
# output.file.area = output.file7a,
# output.file.gridded = output.file7b ,
# shp.file = shp.file,
# file.year = run.years[y],
# write.area =T,
# write.gridded =T,
# t.max.seq = seq(0,30,1)
# )
# }
}
# glorys.heatwave.files = c(list.files(paste0(output.dir,'data/bottom_temp_daily_epu/'),full.names = T),'/home/jcaracappa/EDAB_Dev/jcaracappa/ROMS_NWA/ROMS_daily_epu_1959_1992.csv' )
# glorys.heatwave.data = lapply(glorys.heatwave.files, read.csv) |>
# dplyr::bind_rows() |>
# dplyr::mutate(source_m = dplyr::coalesce(Source, source)) |>
# dplyr::select(-Source,-source) |>
# dplyr::rename(source = 'source_m')
# #check all dates
# all(1959:2025 %in% sort(unique(format(as.Date(glorys.heatwave.data$date),format = '%Y')) ))
# glorys.heatwaves.out = write.csv(glorys.heatwave.data,paste0('/home/jcaracappa/EDAB_Dev/jcaracappa/ROMS_GLORYS_bottom_temp_model_daily_epu_',1959,'_',format(Sys.time(), '%Y'),'.csv'), row.names =F)
#
# season.index = data.frame(month = 1:12, season = rep(c('Winter','Spring','Summer','Fall'),each = 3))
# data = glorys.heatwave.data %>%
# dplyr::mutate(month = as.numeric(format(as.Date(date),format = '%m')),
# year = as.numeric(format(as.Date(date),format = '%Y'))) %>%
# dplyr::left_join(season.index) %>%
# dplyr::group_by(source,year,season,EPU) %>%
# dplyr::summarise(Value = mean(BottomT.mean,na.rm=T))
#
# library(ggplot2)
# ggplot(data, aes(x = year, y = Value, color = source))+
# geom_line()+
# facet_grid(EPU~season)
# ggsave(here::here('GLORYS_bottom_temp_2025.png'),width =12, height = 12)The bottom temperature gridded by season produced results for the South Atlantic. Here is a plot following the same code used in ecodata for the Northeast US.
plot_gridded_bt <- function(bottom_temp_model_gridded_csv, year = NULL, scale = "celsius"){
bottom_temp_model_gridded <- read.csv(bottom_temp_model_gridded_csv) |>
tibble::as_tibble()
xmin = -82.5
xmax = -71
ymin = 23.5
ymax = 36
xlims <- c(xmin, xmax)
ylims <- c(ymin, ymax)
if (is.null(year)) year <- max(bottom_temp_model_gridded$Time)
fix <- dplyr::mutate(dplyr::select(dplyr::filter(bottom_temp_model_gridded,
Time == year), -Time), Var = factor(Var, levels = c("winter",
"spring", "summer", "fall")))
if (scale == "fahrenheit") {
fix <- dplyr::mutate(fix, Value = (9/5) * Value + 32)
label <- "Temp. (°F)"
breaks <- c(32, 41, 50, 59, 68, 77)
labelLegend <- c("32", "41", "50", "59", "68", "77")
limits <- c(31.6, 84.2)
midpoint <- 50
} else {
label <- "Temp. (°C)"
breaks <- c(0, 5, 10, 15, 20, 25)
labelLegend <- c("0", "5", "10", "15", "20", "25")
limits <- c(-0.2, 29)
midpoint <- 10
}
p <- ggplot2::ggplot(data = fix) +
ggplot2::geom_tile(ggplot2::aes(x = Longitude,
y = Latitude, fill = sqrt(Value + 1))) +
ggplot2::geom_sf(data = ecodata::coast,
size = 0.4) +
#ggplot2::geom_sf(data = ne_epu_sf, fill = "transparent", size = 0.4) +
ggplot2::coord_sf(xlim = xlims,
ylim = ylims) +
ggplot2::facet_wrap(Var ~.) +
ecodata::theme_map() +
ggplot2::scale_fill_gradient2(name = label, low = "#000004FF",
mid = "#BB3754FF", high = "#FCFFA4FF",
breaks = sqrt(breaks + 1), limits = sqrt(limits + 1),
labels = labelLegend,
midpoint = sqrt(midpoint + 1)) +
ggplot2::ggtitle(paste("Seasonal Mean Bottom Temperature", year)) +
ggplot2::xlab("Longitude") + ggplot2::ylab("Latitude") +
ggplot2::theme(panel.border = ggplot2::element_rect(colour = "black",
fill = NA, linewidth = 0.75),
legend.key = ggplot2::element_blank(),
axis.title = ggplot2::element_text(size = 11),
strip.background = ggplot2::element_blank(),
strip.text = ggplot2::element_text(hjust = 0),
axis.text = ggplot2::element_text(size = 8),
axis.title.y = ggplot2::element_text(angle = 90)) +
ecodata::theme_title() +
ecodata::theme_ts()
return(p)
}
csv <- here::here("GLORYS/SAtl-out/data/bottom_temp_model_gridded/GLORYS_bottom_temp_model_gridded_2025.csv")
plot_gridded_bt(bottom_temp_model_gridded_csv = csv, scale = "celsius")
Now plotting all seasons after correcting the function for daily GLORYS file input.
csv <- here::here("GLORYS/SAtl-out/data/bottom_temp_model_gridded/GLORYS_bottom_temp_model_gridded_2024.csv")
plot_gridded_bt(bottom_temp_model_gridded_csv = csv, scale = "celsius")
csv <- here::here("GLORYS/SAtl-out/data/bottom_temp_model_gridded/GLORYS_bottom_temp_model_gridded_2023.csv")
plot_gridded_bt(bottom_temp_model_gridded_csv = csv, scale = "celsius")
csv <- here::here("GLORYS/SAtl-out/data/bottom_temp_model_gridded/GLORYS_bottom_temp_model_gridded_2022.csv")
plot_gridded_bt(bottom_temp_model_gridded_csv = csv, scale = "celsius")
csv <- here::here("GLORYS/SAtl-out/data/bottom_temp_model_gridded/GLORYS_bottom_temp_model_gridded_2021.csv")
plot_gridded_bt(bottom_temp_model_gridded_csv = csv, scale = "celsius")
csv <- here::here("GLORYS/SAtl-out/data/bottom_temp_model_gridded/GLORYS_bottom_temp_model_gridded_2020.csv")
plot_gridded_bt(bottom_temp_model_gridded_csv = csv, scale = "celsius")
Now to investigate why thermal habitat is breaking for the South Atlantic. First see if the updated dev branch function runs. It runs and produces gridded files with no outputs. the area files appear to be wrong.
line by line…
It produces output when I only use a single area and don’t attempt to use the areas in the South Atlantic shapefile.
make_thermal_habitat_area = function(input.file, output.file.area = NA, output.file.gridded =NA, shp.file, file.year,t.max.seq, write.area = F, write.gridded = F){
supp.dir = '/Users/sarahgaichas/Documents/Work/SAFMCindicators/READ-EDAB-GLORYS/data-raw/'
#EPU.names = c('MAB','GB','GOM','SS')
EPU.names = c('SAtl') # modify for single large South Atlantic area
#EPU.names = c('Off NC', "Off SC", "Off GA", "all")
#depth.df = data.frame(id = 1:4,
#EPU.names = c('MAB','GB','GOM','SS','all')
depth.df = data.frame(depth.min = c(0,25,100, 0),
depth.max = c(25,100,300, 2000),
depth.name = c('0-25m','25-100m','100-300m', 'AllDepths'))
combs = expand.grid(year = file.year,t.max = t.max.seq,depth.name = depth.df$depth.name,EPU = EPU.names,stringsAsFactors = F)%>%
dplyr::left_join(depth.df)
bathy.shp = terra::rast(paste0(supp.dir,'GLORYS/GLORYS_bathymetry_east_coast_crop.nc'),subds = 'deptho')
# --- Refactor Start: Handle Input Files ---
# Convert input to a single SpatRaster.
# terra::rast() handles both a single multi-layer file OR a vector of single-layer files.
processed_rast = terra::rast(input.file)
# Ensure time dimension is set if missing (common when stacking individual daily files)
# This logic matches layers to days in the specified file.year
if (!all(terra::has.time(processed_rast))) {
# Generate daily sequence for the given year
dates = seq(as.Date(paste0(file.year, "-01-01")), as.Date(paste0(file.year, "-12-31")), by="day")
# If layer count matches the days in the year, assign the time
if (terra::nlyr(processed_rast) == length(dates)) {
terra::time(processed_rast) = dates
}
}
# --- Refactor End ---
out.area.ls = list()
out.gridded.ls = list()
for(i in 1:nrow(combs)){
if(combs$EPU[i] == 'all'){
#area.names = c('MAB','GB','GOM','SS')
area.names = c('Off NC', "Off SC", "Off GA")
}else{
area.names = combs$EPU[i]
}
if(i ==1){
# Refactor: use processed_rast instead of loading from file again
neus.shp = terra::crop(bathy.shp, processed_rast)
}
depth.rast = terra::clamp(neus.shp, lower = combs$depth.min[i], upper = combs$depth.max[i],values =F)
EPU.vect = terra::vect(shp.file)
#area.vect = EPU.vect[which(EPU.vect$EPU %in% area.names)]
area.mask = terra::mask(depth.rast,EPU.vect)#area.vect)
#Mask of area over t.max
# Refactor: Pass the SpatRaster object (processed_rast) instead of filename
area.i = EDABUtilities::mask_nc_2d(data.in = processed_rast,
write.out =F,
shp.file = area.mask,
var.name = 'BottomT',
min.value = combs$t.max[i],
max.value = Inf,
binary = F,
area.names =NA
)
# Refactor: Subset from the main raster object
this.rast = terra::subset(processed_rast, 1)
nd.i = EDABUtilities::make_2d_deg_day_gridded_nc(data.in = area.i,
shp.file = area.mask,
var.name = 'BottomT',
statistic = 'nd',
type = 'above',
ref.value = combs$t.max[i],
area.names = NA
)
shp.area = terra::expanse(area.mask)$area
area.df = terra::expanse(area.i[[1]]) %>%
as.data.frame()%>%
dplyr::mutate(Time = terra::time(area.i[[1]]),
EPU = combs$EPU[i],
Depth = combs$depth.name[i],
Var = paste0('>',combs$t.max[i],'\u00B0C'),
Value = area/ shp.area,
Source = 'GLORYS',
year = combs$year[i],
temp.threshold = combs$t.max[i],
Units = 'Proportion'
)
out.area.ls[[i]] = area.df
out.gridded.ls[[i]] = as.data.frame(nd.i[[1]],cells =T, xy = T) %>%
dplyr::mutate(Time = combs$year[i], EPU = combs$EPU[i], Depth = combs$depth.name[i], Var = combs$t.max[i], Source = 'GLORYS',Units = 'Number of Days')%>%
dplyr::rename(Latitude = 'y', Longitude = 'x', Value = 'sum')%>%
dplyr::select(Time,EPU, Depth, Var,Value,Latitude,Longitude,Source,Units)
print(signif(i/nrow(combs)*100,2))
}
out.area.df = dplyr::bind_rows(out.area.ls)%>%
dplyr::select(Time, EPU, Depth, Var, Value, Source, year, temp.threshold, Units)%>%
dplyr::mutate(Year = format(as.Date(Time),format = '%Y'))%>%
dplyr::group_by(Year,EPU, Depth, Var, temp.threshold, Units,Source)%>%
dplyr::summarise(Value = mean(Value))%>%
dplyr::rename(Time = Year)
out.gridded.df = dplyr::bind_rows(out.gridded.ls)
if(write.area == T){
write.csv(out.area.df, output.file.area,row.names = F)
}
if(write.gridded == T){
write.csv(out.gridded.df,output.file.gridded,row.names = F)
}
if(write.area == F |write.gridded == F){
return(list(thermal.area = out.area.df,thermal.gridded = out.gridded.df))
}
}
input.dir = '/Users/sarahgaichas/Documents/Work/SAFMCindicators/SAFMCindicators/GLORYS/SAtl/'
output.dir = '/Users/sarahgaichas/Documents/Work/SAFMCindicators/SAFMCindicators/GLORYS/SAtl-out/'
supp.dir = '/Users/sarahgaichas/Documents/Work/SAFMCindicators/READ-EDAB-GLORYS/data-raw/'
#print('Using default arguments')
shp.file = paste0(supp.dir,'geometry/SA_EEZ_off_states.shp')
#shp.file = NA
#Get year range
year.start = 2020
year.end = 2025 #format(Sys.time(), '%Y') this would go to now, our data end in 2024
run.years = year.start:year.end
y=1
input.prefix = 'GLORYS_REANALYSIS_'
check.dir = function(file){
if(!dir.exists(dirname(file))){dir.create(dirname(file),recursive =T)}
}
#Produce each year's indices
for(y in 1:length(run.years)){
input.year = run.years[y]
input.diryr = paste0(input.dir, input.year,"/")
input.files = list.files(input.diryr,input.prefix)
input.files.year = as.numeric(gsub(".*(\\d{4}).*", '\\1', input.files))
no.threshold =F
which.files.year = which(input.files.year == run.years[y])
if(length(which.files.year) == 0){
next()
}
this.year.files = paste0(input.diryr,input.files[which.files.year])
# if(any(file.size(this.year.files)==0)){
# next()
# }
#
# if(all(!file.exists(this.year.files))){
# print(paste0('File does not exist: ',this.year.file))
# next()
# }
year.dates = seq.Date(as.Date(paste0(run.years[y],'-01-01')), as.Date(paste0(run.years[y],'-12-31')), by = '1 day')
file.dates <- sub(".*(\\d{4}-\\d{2}-\\d{2}).*", "\\1", this.year.files)
if(!all(year.dates %in% file.dates)){
print(paste0('Missing Dates in ',run.years[y],': ', paste0(year.dates[which(!(year.dates %in% file.dates))],collapse = ', ')))
next()
}
#Make thermal_habitat_area
print(paste0('Starting thermal_habitat_area for year: ',run.years[y], ' at ', Sys.time()))
output.file7a = paste0(output.dir,'data/thermal_habitat_area/GLORYS_thermal_habitat_area_',run.years[y],'.csv')
output.file7b = paste0(output.dir,'data/thermal_habitat_gridded/GLORYS_thermal_habitat_gridded_',run.years[y],'.csv')
check.dir(output.file7a)
if(file.exists(output.file7a) & file.exists(output.file7b) | no.threshold){
#if(file.exists(output.file7b) | no.threshold){
print(paste0('File already exists: ',output.file7a, ' and ', output.file7b))
#print(paste0('File already exists: ', output.file7b))
}else{
#if(any(file.size(this.year.files)==0))
# using local function above from dev branch of EDAB.GLORYS
make_thermal_habitat_area(input.file = this.year.files,
output.file.area = output.file7a,
output.file.gridded = output.file7b,
shp.file = shp.file,
file.year = run.years[y],
write.area =T,
write.gridded =T,
t.max.seq = seq(0,30,1)
)
}
}Test plot, I think this worked
file.names = list.files(here::here("GLORYS/SAtl-out/data/thermal_habitat_gridded"), full.names = T)
# test.nc = ncdf4::nc_open(file.names[1])
# # # test.nc = ncdf4::nc_open('C:/Users/joseph.caracappa/Downloads/thermal_habitat_gridded_1993_SOE2025.nc')
# ncdf4::ncatt_get(test.nc,'nday_10')
# ncdf4::ncatt_get(test.nc,'nday_5')
# ncdf4::nc_close(test.nc)
# #x = terra::rast(file.names[1])
# #x = terra::rast(test.nc,'nday_5')
# #plot(x)
# # test.sds = terra::sds(file.names[1],1)
# # writeCDF(test.sds,here::here('data','SOE','thermal_habitat_checker.nc'))
#
# # convert to tibble and plot as above modifying ecodata::plot_thermal_habitat_gridded code
# ncfile <- file.names[1] # 2020, retrieve the year from the filename to add to the plot
#
# bttib <- tidync::tidync(ncfile) |>
# tidync::hyper_tibble(force = TRUE) |>
# dplyr::select(nday_5, nday_10, nday_15, nday_20, longitude, latitude) |>
# dplyr::mutate(longitude = as.numeric(longitude),
# latitude = as.numeric(latitude))
btthermalgrid <- read.csv(file.names[1])
fix <- btthermalgrid |>
dplyr::filter(Var %in% c(5, 10, 15, 20, 25, 30),
Depth %in% c("AllDepths"))
xmin = -82.5
xmax = -71
ymin = 23.5
ymax = 36
p <- ggplot2::ggplot(fix) +
ggplot2::geom_tile(ggplot2::aes(x = Longitude,
y = Latitude, color = Value, width = 1/12, height = 1/12),
linewidth = 2) +
ggplot2::geom_sf(data = ecodata::coast,
size = 0.4) +
#ggplot2::facet_grid(Depth ~ Var) +
ggplot2::facet_wrap(~Var, nrow = 2) +
ggplot2::scale_color_viridis_c() +
ggplot2::coord_sf(xlim = c(xmin, xmax), ylim = c(ymin, ymax)) +
ggplot2::xlab("") +
ggplot2::ylab("") +
ecodata::theme_ts() +
ecodata::theme_facet() +
ecodata::theme_title() +
ggplot2::ggtitle(paste(unique(fix$Time))) +
ggplot2::theme(legend.position = "bottom")
#return(p)
p
Lets take a look at bottom temps across years in a certain band–days over 25 degrees C
btthermalcomp <- function(filename = filename, daysabove = 0, depth = "AllDepths"){
btthermalgrid <- read.csv(filename)
fix <- btthermalgrid |>
dplyr::filter(Var == daysabove,
Depth == depth)
xmin = -82.5
xmax = -71
ymin = 23.5
ymax = 36
p <- ggplot2::ggplot(fix) +
ggplot2::geom_tile(ggplot2::aes(x = Longitude,
y = Latitude, color = Value, width = 1/12, height = 1/12),
linewidth = 2) +
ggplot2::geom_sf(data = ecodata::coast,
size = 0.4) +
#ggplot2::facet_grid(Depth ~ Var) +
ggplot2::facet_wrap(~Var) +
ggplot2::scale_color_viridis_c() +
ggplot2::coord_sf(xlim = c(xmin, xmax), ylim = c(ymin, ymax)) +
ggplot2::xlab("") +
ggplot2::ylab("") +
ecodata::theme_ts() +
ecodata::theme_facet() +
ecodata::theme_title() +
ggplot2::ggtitle(paste(unique(fix$Time))) +
ggplot2::theme(legend.position = "bottom")
return(p)
}
library(patchwork)
p1 <- btthermalcomp(filename = file.names[1], daysabove = 25) + ggplot2::theme(legend.position = "none")
p2 <- btthermalcomp(filename = file.names[2], daysabove = 25)
p3 <- btthermalcomp(filename = file.names[3], daysabove = 25) + ggplot2::theme(legend.position = "none")
p4 <- btthermalcomp(filename = file.names[4], daysabove = 25) + ggplot2::theme(legend.position = "none")
p5 <- btthermalcomp(filename = file.names[5], daysabove = 25)
p6 <- btthermalcomp(filename = file.names[6], daysabove = 25) + ggplot2::theme(legend.position = "none")
p1 + p2 + p3
p4 + p5 + p6
References
1.
Ross AC, Stock CA, Koul V, Delworth TL, Lu F, Wittenberg A, et al. Dynamically downscaled seasonal ocean forecasts for North American east coast ecosystems. Ocean Science. 2024;20: 1631–1656. doi:10.5194/os-20-1631-2024
2.
Jean-Michel L, Eric G, Romain B-B, Gilles G, Angélique M, Marie D, et al. The Copernicus Global 1/12° Oceanic and Sea Ice GLORYS12 Reanalysis. Frontiers in Earth Science. 2021;9. doi:10.3389/feart.2021.698876
3.
Jones AW, Burchard KA, Mercer AM, Hoey JJ, Morin MD, Gianesin GL, et al. Learning From the Study Fleet: Maintenance of a Large-Scale Reference Fleet for Northeast U.S. Fisheries. Frontiers in Marine Science. 2022;9. Available: https://www.frontiersin.org/articles/10.3389/fmars.2022.869560
4.
R Core Team. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing; 2023. Available: https://www.R-project.org/
5.
Ushey K, Allaire JJ, Tang Y. Reticulate: Interface to ’Python’. 2025. Available: https://CRAN.R-project.org/package=reticulate
6.
Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, et al. Welcome to the tidyverse. Journal of Open Source Software. 2019;4: 1686. doi:10.21105/joss.01686
7.
8.
Pierce D. ncdf4: Interface to Unidata netCDF (Version 4 or Earlier) Format Data Files. 2025. Available: https://CRAN.R-project.org/package=ncdf4
9.
Xie Y. Knitr: A General-Purpose Package for Dynamic Report Generation in R. 2025. Available: https://yihui.org/knitr/
10.
Xie Y. Dynamic Documents with R and knitr. 2nd ed. Boca Raton, Florida: Chapman; Hall/CRC; 2015. Available: https://yihui.org/knitr/
11.
Xie Y. Knitr: A Comprehensive Tool for Reproducible Research in R. In: Stodden V, Leisch F, Peng RD, editors. Implementing Reproducible Computational Research. Chapman; Hall/CRC; 2014.
12.
Allaire JJ, Xie Y, Dervieux C, McPherson J, Luraschi J, Ushey K, et al. Rmarkdown: Dynamic Documents for R. 2024. Available: https://github.com/rstudio/rmarkdown
13.
Xie Y, Allaire JJ, Grolemund G. R Markdown: The Definitive Guide. Boca Raton, Florida: Chapman; Hall/CRC; 2018. Available: https://bookdown.org/yihui/rmarkdown
14.
Xie Y, Dervieux C, Riederer E. R Markdown Cookbook. Boca Raton, Florida: Chapman; Hall/CRC; 2020. Available: https://bookdown.org/yihui/rmarkdown-cookbook
15.
Xie Y. Bookdown: Authoring Books and Technical Documents with R Markdown. 2025. Available: https://github.com/rstudio/bookdown
16.
Xie Y. Bookdown: Authoring Books and Technical Documents with R Markdown. Boca Raton, Florida: Chapman; Hall/CRC; 2016. Available: https://bookdown.org/yihui/bookdown
17.
Beltz B, Beet A, Bastille K, Hardison S. Ecodata: Documentation of Ecosystem Indicator Reporting. 2025. Available: https://noaa-edab.github.io/ecodata/