Hovmoeller dTA, dDIC, and CDR Efficiency Plots
Victoria Froh
03 February, 2025
Last updated: 2025-02-03
Checks: 7 0
Knit directory: oae_ccs_roms/
This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20240919) was run prior to running
the code in the R Markdown file. Setting a seed ensures that any results
that rely on randomness, e.g. subsampling or permutations, are
reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version 639b38d. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for
the analysis have been committed to Git prior to generating the results
(you can use wflow_publish or
wflow_git_commit). workflowr only checks the R Markdown
file, but you know if there are other scripts or data files that it
depends on. Below is the status of the Git repository when the results
were generated:
Ignored files:
Ignored: .Rhistory
Ignored: .Rproj.user/
Untracked files:
Untracked: code/read_data.Rmd
Untracked: oae_ccs_roms.Rproj.RData
Unstaged changes:
Modified: analysis/_site.yml
Deleted: analysis/read_data.Rmd
Modified: code/Workflowr_project_managment.R
Modified: code/regrid_cdr_eff_molar.R
Modified: code/regrid_dTA_dDIC.R
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the repository in which changes were
made to the R Markdown (analysis/hovmo_plots.Rmd) and HTML
(docs/hovmo_plots.html) files. If you’ve configured a
remote Git repository (see ?wflow_git_remote), click on the
hyperlinks in the table below to view the files as they were in that
past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| Rmd | 639b38d | vgfroh | 2025-02-03 | Column integrated plots and hovmoeller plots completed |
Introduction
Plotting dTA, dDIC, and CDR efficiency across depth and time of the OAE addition site. Also plotting their differences from a mean of all 3 phases to compare.
#loading packages
library(tidyverse)
library(data.table)
library(arrow)
library(scales)
# Path to intermediate computation outputs
path_outputs <- "/net/sea/work/vifroh/oae_ccs_roms_data/regrid_2/"
# Path to save practice plots when working on them
path_plots <- "/net/sea/work/vifroh/test_plots/"
# loading in previous saved depth integrated data for hovmo plots
lanina_depthint <- read_feather(
paste0(path_outputs, "lanina_depthintRG2.feather"))
neutral_depthint <- read_feather(
paste0(path_outputs, "neutral_depthintRG2.feather"))
elnino_depthint <- read_feather(
paste0(path_outputs, "elnino_depthintRG2.feather"))Regular Hovmoeller Plots
Change in total (added) alkalinity
# formatting tables
phase_data <- list(lanina_depthint, neutral_depthint, elnino_depthint)
phase_data <- phase_data %>%
lapply(function(table) {
setDT(table)
table[, depth := as.numeric(depth)]
table[, dTA_sum := as.numeric(as.character(dTA_sum))]
table[, CDR_eff := as.numeric(as.character(CDR_eff))]
table[, time := as.Date(paste0(time, "-01"))]
return(table)
})
# return tables to each item
lanina_depthint <- phase_data[[1]]
neutral_depthint <- phase_data[[2]]
elnino_depthint <- phase_data[[3]]
# subset first two years
lanina_depthint_sub <- lanina_depthint[time <= as.Date("2000-05-01")]
neutral_depthint_sub <- neutral_depthint[time <= as.Date("2005-05-01")]
elnino_depthint_sub <- elnino_depthint[time <= as.Date("2017-05-01")]
# # removing negative data points and transform
# lanina_depthint_sub <- lanina_depthint_sub[dTA_sum > 0]
# lanina_depthint_sub[, dTA_log := log10(dTA_sum)]
# setting boundaries of bins
depth_ranges <- unique(lanina_depthint_sub[, .(depth)])
setorder(depth_ranges, depth)
depth_ranges[, ":=" (
depth_upper = (shift(depth, type = "lag", fill = 0) + depth) / 2, # set shallower boundary
depth_lower = (shift(depth, type = "lead", fill = max(depth)) + depth) / 2
)] # deeper bound
# merge depth ranges back to tables and divide into surface and subsurface
phase_data <- list(lanina_depthint_sub, neutral_depthint_sub, elnino_depthint_sub)
phase_data <- phase_data %>%
lapply(function(table) {
table <- merge(table, depth_ranges, by = "depth", all.x = TRUE)
table[, category := fifelse(depth < 100, "Surface", "Subsurface")]
table$category <- factor(table$category, levels = c("Surface", "Subsurface"))
return(table)
})
# return tables to each item
lanina_depthint_sub <- phase_data[[1]]
neutral_depthint_sub <- phase_data[[2]]
elnino_depthint_sub <- phase_data[[3]]
# plotting delta dTA from mean
phase_titles <- c("La Niña", "Neutral", "El Niño")
phases <- c("lanina", "neutral", "elnino")
create_hovmo_dTA <- function(data, title_text, phase) {
plot <- ggplot(data, aes(x = time, y = depth, fill = dTA_sum)) +
geom_rect(aes(ymin = depth_upper, ymax = depth_lower,
xmin = time - 16, xmax = time + 16)) +
scale_y_reverse(expand = c(0, 0)) +
scale_x_date(
date_breaks = "4 months", # Tick every 4 months
date_labels = "%Y-%m", # Show only the year and month
expand = c(0, 0)
) +
scale_fill_viridis_c(guide = guide_colorbar(
barwidth = unit(0.7, "cm"), barheight = unit(6, "cm")),
limits = c(-926000, 3443000000)) +
# limits = c(0, 3443000000)) +
facet_wrap(~category, scales = "free_y", ncol = 1) + # splits into two panels
labs(x = "Time",
y = "Depth (m)",
fill = "Added Alkalinity (mol/m)",
title = paste0(title_text, " Phase (10N to 60N and -170E to -85E)")) +
theme_bw() +
theme(strip.text = element_blank(),
axis.ticks.length = unit(0.15, "inches"),
panel.border = element_blank(), # Remove box around facets
panel.grid = element_blank(),
panel.spacing = unit(0.5, "lines"),
strip.background = element_blank())
print(plot)
# save plot
ggsave(paste0(path_plots, phase, "_hovmo_dTA.png"), plot = plot,
width = 8, height = 6, dpi = 300)
}
lapply(seq_along(phase_data), function(i) {
create_hovmo_dTA(phase_data[[i]], phase_titles[i], phases[i])
})
[[1]]
[1] "/net/sea/work/vifroh/test_plots/lanina_hovmo_dTA.png"
[[2]]
[1] "/net/sea/work/vifroh/test_plots/neutral_hovmo_dTA.png"
[[3]]
[1] "/net/sea/work/vifroh/test_plots/elnino_hovmo_dTA.png"rm(create_hovmo_dTA)
gc() used (Mb) gc trigger (Mb) max used (Mb)
Ncells 1703113 91.0 3354720 179.2 2208121 118.0
Vcells 3118613 23.8 8388608 64.0 6554929 50.1Change in DIC
# building off of data ran in chunk above
# plotting dDIC
create_hovmo_dDIC <- function(data, title_text, phase) {
plot <- ggplot(data, aes(x = time, y = depth, fill = dDIC_sum)) +
geom_rect(aes(ymin = depth_upper, ymax = depth_lower,
xmin = time - 16, xmax = time + 16)) +
scale_y_reverse(expand = c(0, 0)) +
scale_x_date(
date_breaks = "4 months", # Tick every 4 months
date_labels = "%Y-%m", # Show only the year and month
expand = c(0, 0)
) +
scale_fill_viridis_c(guide = guide_colorbar(
barwidth = unit(0.7, "cm"), barheight = unit(6, "cm")),
limits = c(-480, 2008000000)) +
# limits = c(0, 2008000000)) +
facet_wrap(~category, scales = "free_y", ncol = 1) + # splits into two panels
labs(x = "Time",
y = "Depth (m)",
fill = "Change in DIC (mol/m)",
title = paste0(title_text, " Phase (10N to 60N and -170E to -85E)")) +
theme_bw() +
theme(strip.text = element_blank(),
axis.ticks.length = unit(0.15, "inches"),
panel.border = element_blank(), # Remove box around facets
panel.grid = element_blank(),
panel.spacing = unit(0.5, "lines"),
strip.background = element_blank())
print(plot)
# save plot
ggsave(paste0(path_plots, phase, "_hovmo_dDIC.png"), plot = plot,
width = 8, height = 6, dpi = 300)
}
lapply(seq_along(phase_data), function(i) {
create_hovmo_dDIC(phase_data[[i]], phase_titles[i], phases[i])
})
[[1]]
[1] "/net/sea/work/vifroh/test_plots/lanina_hovmo_dDIC.png"
[[2]]
[1] "/net/sea/work/vifroh/test_plots/neutral_hovmo_dDIC.png"
[[3]]
[1] "/net/sea/work/vifroh/test_plots/elnino_hovmo_dDIC.png"rm(create_hovmo_dDIC)
gc() used (Mb) gc trigger (Mb) max used (Mb)
Ncells 1708429 91.3 3354720 179.2 3354720 179.2
Vcells 3132752 24.0 8388608 64.0 6554929 50.1CDR Efficiency (dDIC/dTA)
# building off of data ran in chunk above
# plotting CDReff
create_hovmo_CDReff <- function(data, title_text, phase) {
plot <- ggplot(data, aes(x = time, y = depth, fill = CDR_eff)) +
geom_rect(aes(ymin = depth_upper, ymax = depth_lower,
xmin = time - 16, xmax = time + 16)) +
scale_y_reverse(expand = c(0, 0)) +
scale_x_date(
date_breaks = "4 months", # Tick every 4 months
date_labels = "%Y-%m", # Show only the year and month
expand = c(0, 0)
) +
scale_fill_viridis_c(guide = guide_colorbar(
barwidth = unit(0.7, "cm"), barheight = unit(6, "cm")),
limits = c(0, 1)) + # cutting off weird negative or >1 CDR effs
facet_wrap(~category, scales = "free_y", ncol = 1) + # splits into two panels
labs(x = "Time",
y = "Depth (m)",
fill = "CDR Efficiency (Fraction))",
title = paste0(title_text, " Phase (10N to 60N and -170E to -85E)")) +
theme_bw() +
theme(strip.text = element_blank(),
axis.ticks.length = unit(0.15, "inches"),
panel.border = element_blank(), # Remove box around facets
panel.grid = element_blank(),
panel.spacing = unit(0.5, "lines"),
strip.background = element_blank())
print(plot)
# save plot
ggsave(paste0(path_plots, phase, "_hovmo_CDReff.png"), plot = plot,
width = 8, height = 6, dpi = 300)
}
lapply(seq_along(phase_data), function(i) {
create_hovmo_CDReff(phase_data[[i]], phase_titles[i], phases[i])
})
[[1]]
[1] "/net/sea/work/vifroh/test_plots/lanina_hovmo_CDReff.png"
[[2]]
[1] "/net/sea/work/vifroh/test_plots/neutral_hovmo_CDReff.png"
[[3]]
[1] "/net/sea/work/vifroh/test_plots/elnino_hovmo_CDReff.png"rm(create_hovmo_CDReff)
gc() used (Mb) gc trigger (Mb) max used (Mb)
Ncells 1708833 91.3 3354720 179.2 3354720 179.2
Vcells 3134126 24.0 8388608 64.0 6554929 50.1Hovmoeller plots of the differences from the means of all 3 phases
Difference in added alkalinity from the phases’ mean
# building off of data ran in chunk above
# convert time column to months since addition begins
phase_data <- phase_data %>%
lapply(function(table) {
setorder(table, time)
table[, month := .GRP, by = time] # gives index in order to each unique time
return(table)
})
# return tables to each item
lanina_depthint_sub <- phase_data[[1]]
neutral_depthint_sub <- phase_data[[2]]
elnino_depthint_sub <- phase_data[[3]]
# combining all three tables in to one
phase_list <- Map(function(table, phase) {
table[, phase := phase]
table[, .(month, depth, depth_upper, depth_lower, category, dTA_sum, dDIC_sum,
CDR_eff, phase)]
}, phase_data, phases)
phase_fulldata <- rbindlist(phase_list)
# calculating mean dTA, dDIC, and CDR-eff across all 3 phases
phase_fulldata[, ":=" (
dTA_mean = mean(dTA_sum, na.rm = TRUE),
dDIC_mean = mean(dDIC_sum, na.rm = TRUE),
CDR_effmean = mean(CDR_eff[CDR_eff >= 0 & CDR_eff <= 1] , na.rm = TRUE)
), by = .(month, depth)]
create_hovmo_ddTA <- function(phase_name, title_text) {
phase_dt <- phase_fulldata[phase_fulldata$phase == phase_name,]
plot <- ggplot(phase_dt, aes(x = month, y = depth, fill = dTA_sum - dTA_mean)) +
geom_rect(aes(xmin = month - 1, xmax = month,
ymin = depth_upper, ymax = depth_lower)) +
scale_y_reverse(expand = c(0, 0)) +
scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0,
guide = guide_colorbar(
barwidth = unit(0.7, "cm"), barheight = unit(6, "cm")),
limits = c(-410000000, 430000000)) +
facet_wrap(~category, scales = "free_y", ncol = 1) + # splits into two panels
labs(x = "Months Since Addition Began",
y = "Depth (m)",
fill = "Added Alkalinity Difference (mol/m)",
title = paste0(title_text, " Phase, Added Alkalinity Difference from Mean")) +
theme_bw() +
theme(strip.text = element_blank(),
axis.ticks.length = unit(0.15, "inches"),
panel.border = element_blank(), # Remove box around facets
panel.grid = element_blank(),
panel.spacing = unit(0.5, "lines"),
strip.background = element_blank())
print(plot)
# save plot
ggsave(paste0(path_plots, phase_name, "_hovmo_ddTA.png"), plot = plot,
width = 8, height = 6, dpi = 300)
}
lapply(seq_along(phases), function(i) {
create_hovmo_ddTA(phases[i], phase_titles[i])
})
[[1]]
[1] "/net/sea/work/vifroh/test_plots/lanina_hovmo_ddTA.png"
[[2]]
[1] "/net/sea/work/vifroh/test_plots/neutral_hovmo_ddTA.png"
[[3]]
[1] "/net/sea/work/vifroh/test_plots/elnino_hovmo_ddTA.png"rm(phase_list, create_hovmo_ddTA)
gc() used (Mb) gc trigger (Mb) max used (Mb)
Ncells 1711640 91.5 3354720 179.2 3354720 179.2
Vcells 3189637 24.4 8388608 64.0 6554929 50.1Difference in dDIC from the phases’ mean
# building off of data ran in chunk above
# plotting delta dDIC from mean
create_hovmo_ddDIC <- function(phase_name, title_text) {
phase_dt <- phase_fulldata[phase_fulldata$phase == phase_name,]
plot <- ggplot(phase_dt, aes(x = month, y = depth, fill = dDIC_sum - dDIC_mean)) +
geom_rect(aes(xmin = month - 1, xmax = month,
ymin = depth_upper, ymax = depth_lower)) +
scale_y_reverse(expand = c(0, 0)) +
scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0,
guide = guide_colorbar(
barwidth = unit(0.7, "cm"), barheight = unit(6, "cm")),
limits = c(-200000000, 200000000)) +
facet_wrap(~category, scales = "free_y", ncol = 1) + # splits into two panels
labs(x = "Months Since Addition Began",
y = "Depth (m)",
fill = "dDIC Difference (mol/m)",
title = paste0(title_text, " Phase, dDIC Difference from Mean")) +
theme_bw() +
theme(strip.text = element_blank(),
axis.ticks.length = unit(0.15, "inches"),
panel.border = element_blank(), # Remove box around facets
panel.grid = element_blank(),
panel.spacing = unit(0.5, "lines"),
strip.background = element_blank())
print(plot)
# save plot
ggsave(paste0(path_plots, phase_name, "_hovmo_ddDIC.png"), plot = plot,
width = 8, height = 6, dpi = 300)
}
lapply(seq_along(phases), function(i) {
create_hovmo_ddDIC(phases[i], phase_titles[i])
})
[[1]]
[1] "/net/sea/work/vifroh/test_plots/lanina_hovmo_ddDIC.png"
[[2]]
[1] "/net/sea/work/vifroh/test_plots/neutral_hovmo_ddDIC.png"
[[3]]
[1] "/net/sea/work/vifroh/test_plots/elnino_hovmo_ddDIC.png"rm(create_hovmo_ddDIC)
gc() used (Mb) gc trigger (Mb) max used (Mb)
Ncells 1711801 91.5 3354720 179.2 3354720 179.2
Vcells 3186651 24.4 8388608 64.0 6554929 50.1Difference in CDR efficiency from the phases’ mean; wonky efficiencies outside the range from 0 to 1 were filtered out (these occurred primarily below 150m)
# building off of data ran in chunk above; CDR efficiencies outside 0-1 were
# excluded from the mean calculation since they are casued by errors
# plotting delta dDIC from mean
create_hovmo_dCDReff <- function(phase_name, title_text) {
phase_dt <- phase_fulldata[phase_fulldata$phase == phase_name,]
# filter out negative dTA/dDIC/wrong CDReff as these will make weird/incorrect efficiency
phase_dt <- phase_dt[dTA_sum < 0, dTA_sum := NA]
phase_dt <- phase_dt[dDIC_sum < 0, dDIC_sum := NA]
phase_dt <- phase_dt[CDR_eff < 0 | CDR_eff > 1, CDR_eff := NA]
plot <- ggplot(phase_dt, aes(x = month, y = depth, fill = CDR_eff - CDR_effmean)) +
geom_rect(aes(xmin = month - 1, xmax = month,
ymin = depth_upper, ymax = depth_lower)) +
scale_y_reverse(expand = c(0, 0)) +
scale_fill_gradient2(na.value = "lightgray", low = "blue", mid = "white",
high = "red", midpoint = 0,
guide = guide_colorbar(
barwidth = unit(0.7, "cm"), barheight = unit(6, "cm")),
limits = c(-0.152, 0.152)) +
facet_wrap(~category, scales = "free_y", ncol = 1) + # splits into two panels
labs(x = "Months Since Addition Began",
y = "Depth (m)",
fill = "CDR Efficiency Difference",
title = paste0(title_text, " Phase, CDR Efficiency Difference from Mean")) +
theme_bw() +
theme(strip.text = element_blank(),
axis.ticks.length = unit(0.15, "inches"),
panel.border = element_blank(), # Remove box around facets
panel.grid = element_blank(),
panel.spacing = unit(0.5, "lines"),
strip.background = element_blank())
print(plot)
# save plot
ggsave(paste0(path_plots, phase_name, "_hovmo_dCDReff.png"), plot = plot,
width = 8, height = 6, dpi = 300)
}
lapply(seq_along(phases), function(i) {
create_hovmo_dCDReff(phases[i], phase_titles[i])
})
[[1]]
[1] "/net/sea/work/vifroh/test_plots/lanina_hovmo_dCDReff.png"
[[2]]
[1] "/net/sea/work/vifroh/test_plots/neutral_hovmo_dCDReff.png"
[[3]]
[1] "/net/sea/work/vifroh/test_plots/elnino_hovmo_dCDReff.png"rm(create_hovmo_dCDReff)
gc() used (Mb) gc trigger (Mb) max used (Mb)
Ncells 1713134 91.5 3354720 179.2 3354720 179.2
Vcells 3191967 24.4 8388608 64.0 6554929 50.1
sessionInfo()R version 4.4.2 (2024-10-31)
Platform: x86_64-pc-linux-gnu
Running under: openSUSE Leap 15.6
Matrix products: default
BLAS/LAPACK: /usr/local/OpenBLAS-0.3.28/lib/libopenblas_haswellp-r0.3.28.so; LAPACK version 3.12.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
time zone: Europe/Zurich
tzcode source: system (glibc)
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] scales_1.3.0 arrow_18.1.0.1 data.table_1.16.2 lubridate_1.9.3
[5] forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4 purrr_1.0.2
[9] readr_2.1.5 tidyr_1.3.1 tibble_3.2.1 ggplot2_3.5.1
[13] tidyverse_2.0.0 workflowr_1.7.1
loaded via a namespace (and not attached):
[1] gtable_0.3.6 xfun_0.49 bslib_0.8.0 processx_3.8.4
[5] callr_3.7.6 tzdb_0.4.0 vctrs_0.6.5 tools_4.4.2
[9] ps_1.8.1 generics_0.1.3 fansi_1.0.6 pkgconfig_2.0.3
[13] assertthat_0.2.1 lifecycle_1.0.4 compiler_4.4.2 farver_2.1.2
[17] git2r_0.35.0 textshaping_0.4.0 munsell_0.5.1 getPass_0.2-4
[21] httpuv_1.6.15 htmltools_0.5.8.1 sass_0.4.9 yaml_2.3.10
[25] later_1.4.1 pillar_1.9.0 jquerylib_0.1.4 whisker_0.4.1
[29] cachem_1.1.0 tidyselect_1.2.1 digest_0.6.37 stringi_1.8.4
[33] labeling_0.4.3 rprojroot_2.0.4 fastmap_1.2.0 grid_4.4.2
[37] colorspace_2.1-1 cli_3.6.3 magrittr_2.0.3 utf8_1.2.4
[41] withr_3.0.2 promises_1.3.2 bit64_4.5.2 timechange_0.3.0
[45] rmarkdown_2.29 httr_1.4.7 bit_4.5.0 ragg_1.3.3
[49] hms_1.1.3 evaluate_1.0.1 knitr_1.49 viridisLite_0.4.2
[53] rlang_1.1.4 Rcpp_1.0.13-1 glue_1.8.0 rstudioapi_0.17.1
[57] jsonlite_1.8.9 R6_2.5.1 systemfonts_1.1.0 fs_1.6.5