Temperature-Salinity Diagram
Temperature-salinity diagrams (TS-plots) are used in physical oceanography to identify water masses. Since water masses also have distinct oxygen concentrations and exist along particular isolines of potential density (isopycnals), it is helpful to combine all four parameters into one plot for the purpose of water mass identification.
This function allows you to create TS-plots with isopycnals (using potential density) and dissolved oxygen concentration as the z-axis. This is an adaptation of the ggTS function created by David Kaiser (http://doi.org/10.5281/zenodo.3901308). Here, I've used the viridis palette for the z-axis for accessibility. An extensive list of colour palettes in R, including more information about colour-blind friendly palettes and palettes suitable for black and white displays, can be found here.
The function ggTS requires vectors of equal length of practical salinity (sal), absolute salinity (abs.sal), in-situ temperature (temp), potential temperature (pot.temp), and sea pressure (pres). Optional arguments include reference pressure for the potential density function (default = 0); a vector for the z-axis, specifically oxygen concentration (col.par).
ggTS <- function(sal, abs.sal, temp, pot.temp, pres, reference.p = 0, col.par = NA) {
# load packages
library(gsw)
library(ggplot2) Create a long table, TSpd, which is used to calculate potential density from absolute salinity and in-situ temperature. Create a second long table, TS, which will be used for the plotting of potential temperature and practical salinity with isopycnals and a z-axis. I prefer displaying potential density anomaly so I've subtracted 1000 from the potential density.
# use expand.grid to make long tables
TSpd <- expand.grid(
abs.sal = seq(floor(min(abs.sal, na.rm = TRUE)),ceiling(max(abs.sal, na.rm = TRUE)),length.out = 100),
temp = seq(floor(min(temp, na.rm = TRUE)),ceiling(max(temp, na.rm = TRUE)),length.out = 100),
pres = seq(floor(min(pres, na.rm = TRUE)),ceiling(max(pres, na.rm = TRUE)),length.out = 100))
TS <- expand.grid(
sal = seq(floor(min(sal, na.rm = TRUE)), ceiling(max(sal, na.rm = TRUE)),length.out = 100),
pot.temp = seq(floor(min(pot.temp, na.rm = TRUE)),ceiling(max(pot.temp, na.rm = TRUE)),length.out = 100),
pres = seq(floor(min(pres, na.rm = TRUE)),ceiling(max(pres, na.rm = TRUE)),length.out = 100))
TSpd$density <- gsw_pot_rho_t_exact(SA = TSpd$abs.sal, t = TSpd$temp, p = TSpd$pres, p_ref = reference.p) -
1000
TS$density <- TSpd$density explain next block
h.isopycnals <- subset(TS, TS$sal == ceiling(max(TS$sal)) &
round(TS$density, 1) %in% seq(min(round(TS$density * 2)/2,
na.rm = TRUE), max(round(TS$density * 2)/2, na.rm = TRUE), by = 0.5))
if (nrow(h.isopycnals) > 0) {
h.isopycnals$density <- round(h.isopycnals$density, 1)
h.isopycnals <- aggregate(pot.temp ~ density, h.isopycnals, mean)
}
if (nrow(h.isopycnals) == 0) {
rm(h.isopycnals)
v.isopycnals <- subset(TS, pot.temp == ceiling(max(TS$pot.temp)) &
round(TS$density, 1) %in% seq(min(round(TS$density * 2)/2),
max(round(TS$density * 2)/2), by = 0.5))
v.isopycnals$density <- round(v.isopycnals$density, 1)
v.isopycnals <- aggregate(sal ~ density, v.isopycnals, mean)
}Next, the TS-plot is created using ggplot2. Since col.par is not a required argument, the plot can be created with or without the oxygen concentration (or other required z-axis) overlay.
# create a data frame to relate the z-axis to the other axes
data <- data.frame(sal, pot.temp, col.par)
# plot
p1 <- ggplot() +
# add density contours
geom_contour(data = TS, aes(x = sal, y = pot.temp, z = density), col = "grey",
linetype = "dashed", breaks = seq(min(round(TS$density * 2)/2, na.rm = TRUE),
max(round(TS$density * 2)/2, na.rm = TRUE), by = 0.5)) +
# plot other points with existing oxygen concentration data
geom_point(data = data[!is.na(data$col.par), ], aes(sal, pot.temp, color = col.par), size = 2.5) +
# add in potential density (sigma-theta) symbol
annotate(geom = "text", x = floor(min(TS$sal, na.rm = TRUE)) + 0.6, y = ceiling(max(TS$pot.temp,
na.rm = TRUE)), hjust = "inward", vjust = "inward", color = "grey60",
size = 14, label = expression(italic("σ")[italic("θ")]),
parse = T) +
scale_x_continuous(name = "Practical Salinity [PSU]", expand = c(0, 0), limits = c(floor(min(TS$sal,
na.rm = TRUE)) + 0.5, ceiling(max(TS$sal, na.rm = TRUE)))) +
scale_y_continuous(name = "Potential Temperature [°C]",
expand = c(0, 0), limits = c(floor(min(TS$pot.temp, na.rm = TRUE)) + 0, ceiling(max(TS$pot.temp,
na.rm = TRUE)))) +
scale_colour_gradientn(colours = viridis_pal()(10), # viridis is a colour blind friendly palette
guide = "colorbar", limits = c(3.5,7.5), oob = squish,
na.value = 'black') +
guides(fill = guide_colourbar(barwidth = 0.8, barheight = 20))+
labs(colour = expression("Oxygen \n[ml/l]")) +
theme_classic() +
theme(text = element_text(size = 25))Add potential density labels to the contours at the top and/or right edges of the plot.
if (exists("h.isopycnals")) {
p1 <- p1 + geom_text(data = h.isopycnals, aes(x = ceiling(max(TS$sal)), y = pot.temp, label = density),
hjust = "inward", vjust = -0.65, col = "grey", size = 6)
}
if (exists("v.isopycnals")) {
p1 <- p1 + geom_text(data = v.isopycnals, aes(x = sal, y = ceiling(max(TS$pot.temp)), label = density),
vjust = "inward", hjust = 0, col = "grey", size = 6)
}
return(p1) # display the ggplot
} # close functionAn example TS-plot:
Adapted from : Kaiser, David. (2020, June 19). Davidatlarge/ggTS: ggTS first release (Version v1.0.0). Zenodo. http://doi.org/10.5281/zenodo.3901308
The following adaptations were made:
- Adapted the arguments, long tables, and density calculation to use potential density instead of in-situ denstiy
- Specifically adapted the z-axis for oxygen concentrations
- Changed the colour palette to viridis
- Other minor plotting changes