T–S diagram（温盐图解）是研究海洋水团、海水混合时的一种图解，物理海洋学上经常用到，以温度（T）为纵坐标，盐度（S）为横坐标。绘制T–S diagram常用ODV软件，然而ODV实在太......。下面提供R语言绘制法。
#先加载包
library(testthat); library(gsw); library(oce)
library(ggplot2); library(stringr); library(dplyr); library(plyr)
#选择GSW；Oceanographic Analysis with R: The oceanographic community has a choice of two sets of formulae for calculating seawater properties: the UNESCO formulation, popularized in the 1980s, and the Gibbs-SeaWater (GSW) formulation, proposed in 2010.
options(oceEOS = 'gsw')
#读取CTD数据（cnv格式）
setwd("D:/R/Data/Ocean/ctd") #该文件夹共4个文件，"S1.cnv" "S2.cnv" "S3.cnv" "S4.cnv"
temp <- list.files( pattern = ".cnv")
下面提供两种绘制方法
方案一：oce包绘图
ctds<- vector("list", length(temp))
for (i in 1: length(temp)) {
Cast <- read.oce ( temp[i] )
#str(Cast) #查看数据结构, cnv文件由metadata、data、processingLog三部分构成。
#读取文件名和经纬度信息；下面4行代码需根据cnv的不同格式作相应调整。
Filename <- str_replace(temp[i], '\\.cnv', '') #字符串替换
Geo <- unlist (str_split(Cast@metadata$header[20 : 21],pattern = ' '))
Lon <- as.numeric(Geo[3]) + as.numeric(Geo[4])/60 + as.numeric(Geo[5]) / 3600
Lat <- as.numeric(Geo[8]) + as.numeric(Geo[9])/60 + as.numeric(Geo[10]) / 3600
#读取数据
ctds[[i]] <- as.ctd( salinity = Cast@data$salinity,
temperature = Cast@data$temperature,
pressure = Cast@data$pressure,
type = "SBE",
station = Filename,
longitude = Lon,
latitude = Lat,
deploymentType = "profile")
rm( Filename, Geo, Lon, Lat)
}
#合并所有CTD数据
Data0 <- as.section(ctds)
summary(Data0)
#绘图
plotTS(Data0, nlevels = 9, grid = FALSE, bg = "transparent",
pch = 16, cex = 0.5,
col.rho = "lightgray", cex.rho = 1, lty.rho = 1)
方案二：ggplot2包绘图
#建立空白数据框
Data1 <- data.frame(Station = character(0),
Lon = numeric(0),
Lat = numeric(0),
Pressure = numeric(0),
Temperature = numeric(0),
Salinity = numeric(0),
SA = numeric(0),
CT = numeric(0))
for (i in 1: length(temp)) {
Cast <- read.oce (temp[i])
#str(Cast) #查看数据结构, cnv文件由metadata、data、processingLog三部分构成。
#读取文件名和经纬度信息；下面4行代码需根据cnv的不同格式作相应调整。
Filename <- str_replace(temp[i], '\\.cnv', '') #字符串替换
Geo <- unlist (str_split(Cast@metadata$header[20 : 21],pattern = ' '))
Lon <- as.numeric(Geo[3]) + as.numeric(Geo[4])/60 + as.numeric(Geo[5]) / 3600
Lat <- as.numeric(Geo[8]) + as.numeric(Geo[9])/60 + as.numeric(Geo[10]) / 3600
#提取数据到data.frame
data.tem <- data.frame(Station = Filename,
Lon = Lon,
Lat = Lat,
Pressure = Cast@data$pressure,
Temperature = Cast@data$temperature,
Salinity = Cast@data$salinity)
#计算Absolute Salinity和Conservative Temperature
data.tem$SA <- gsw_SA_from_SP(data.tem$Salinity, data.tem$Pressure, longitude = Lon, latitude = Lat)
data.tem$CT <- gsw_CT_from_t (SA = data.tem$SA, t = data.tem$Temperature, p = data.tem$Pressure)
#合并数据
Data1 <- rbind.fill(Data1, data.tem)
rm( Filename, Geo, Lon, Lat)
}
#构建σ0等密度线数据
# make TS long table
TS <- expand.grid(
SA = seq( floor( min( Data1$SA )), ceiling( max( Data1$SA )), length.out = 100),
CT = seq( floor( min( Data1$CT )) - 3, ceiling( max( Data1$CT )) + 3, length.out = 100)
) #为了显示最上和最下的等密度线，CT的范围选择+-3
TS$Density <- gsw_rho(TS$SA, TS$CT , 0) - 1000
#选择拟绘出的等密度线
isopycnals <- subset(TS,
round(SA,1) == (ceiling(max(TS$SA))-0.8) & #等密度线加标注的位置的x轴坐标，数值0.8根据图片效果作相应调整
round(Density,1) %in% seq(min(round(TS$Density*2)/2),
max(round(TS$Density*2)/2),
by = 1)) #选择plot将绘出的等密度线上的数据
isopycnals$Density <- round(isopycnals$Density, 1) #保留一位小数
isopycnals <- aggregate(CT ~ Density, isopycnals, mean) #相同Density的CT值求均值
#绘图
p <- ggplot() +
geom_contour(data =TS, aes(x = SA, y = CT, z = Density), col = "grey", linetype = "dashed",
breaks = seq(min(round(TS$Density*2)/2), max(round(TS$Density*2)/2), by = 1)) +
geom_text(data = isopycnals, aes(x = 35.2, y = CT, label = Density),
hjust = "inward", vjust = 0, col = "grey60", angle = 15 ) +
geom_point(data=Data1, aes(SA, CT, col = Station))
p;
#细节完善
p + scale_x_continuous(name = "Absolute Salinity [g/kg]", limits = c(33.4, 35.2)) +
scale_y_continuous(name = "Conservative Temperature [°C]", limits = c(0, 32),breaks = seq(5, 30, 5)) +
theme_bw() +
annotate(geom = "text", x = 33.7, y = 32, label = "20", hjust = "inward", vjust = 0, col = "grey60", angle = 15) +
annotate(geom = "text", x = 35.1, y = 32, label = "21", hjust = "inward", vjust = 0, col = "grey60", angle = 15) +
theme(text = element_text(size = 14), axis.text = element_text(size = 12),
panel.grid = element_blank())
