##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
## Registered S3 methods overwritten by 'stars':
## method from
## st_bbox.SpatRaster sf
## st_crs.SpatRaster sf
##
## Attaching package: 'matrixStats'
## The following object is masked from 'package:dplyr':
##
## count
setwd("C:/Work/MAST467/Code/R")
load("C:/Work/MAST467/Code/R/OMG_ENV.RData")
# Quality Control.
# Removes bad Temperature and Salinity data. Removes all data <= 2m. Replaces
# all -99s with R's native missing data "NA". Rebuilds yearly data.frames with
# QA checked data values.
# Thresholds for bad data are, any -99, T <= -2, and S >= 50.
# TO DO: Make thresholds in to variables in the Edit Section.
# EDIT SECTION -----------------------------------------------------------------
QA = 1 # SET TO 1 IF QA PROCESSING IS DESIRED.
# ------------------------------------------------------------------------------
if (QA == 1) {
OMGALL["D"][OMGALL["D"] <= 2] <- -99
TOSS_idx = which(OMGALL$D == -99,OMGALL$T <= -2 && OMGALL$T >= -98,OMGALL$S >= 50)
OMGALL <- OMGALL[-c(TOSS_idx), ]
OMGALL["T"][OMGALL["T"] == -99] <- NA
OMGALL["S"][OMGALL["S"] == -99] <- NA
OMG2021 = OMGALL[OMGALL$DATE > "2021-01-01" & OMGALL$DATE < "2021-12-31",]
OMG2020 = OMGALL[OMGALL$DATE > "2020-01-01" & OMGALL$DATE < "2020-12-31",]
OMG2019 = OMGALL[OMGALL$DATE > "2019-01-01" & OMGALL$DATE < "2019-12-31",]
OMG2018 = OMGALL[OMGALL$DATE > "2018-01-01" & OMGALL$DATE < "2018-12-31",]
OMG2017 = OMGALL[OMGALL$DATE > "2017-01-01" & OMGALL$DATE < "2017-12-31",]
rm(TOSS_idx,QA)
}
# Singular Data.
# Records the deepest depth of each profile and re-indexes the OMGALL data set to retain useful information. New data set is a placeholder for any value unique to one cast that doesn't change with depth.
a = diff(OMGALL$D)
b = which(a < 0)
BTMALL_idx = append(b,length(a))
BTMALL = OMGALL[BTMALL_idx,]
rm(a,b)
# MAPS. I am using the ggOceanMaps package for these maps.
study_region = data.frame(lon = c(-90, -90, -50, -50), lat = c(76, 72, 72, 76))
GREENLAND = basemap(limits = c(-70,-20,58,85),rotate = TRUE, bathymetry = TRUE ,glaciers = TRUE) +
ggtitle("Greenland","Disko Island to Melville Bay (Region 3)") +
theme(plot.title = element_text(hjust = 0.5)) +
geom_spatial_polygon(data = study_region, aes(x = lon, y = lat), alpha = 0.5, fill = "red", color = "red")
## Using lon and lat as longitude and latitude columns, respectively.
## projection transformed from +init=epsg:4326 to +proj=stere +lat_0=90 +lat_ts=71 +lon_0=-45 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
GREENLAND
## Assuming `crs = 4326` in stat_spatial_identity()
REGION3 = basemap(limits = c(-71,-51,71,77),rotate = TRUE, bathymetry = TRUE, glaciers = TRUE) +
ggtitle("Spatial Distribution of CTD Measurements in Region 3 (2017 - 2021)") +
theme(plot.title = element_text(hjust = 0.5)) +
geom_spatial_point(data = BTMALL, aes(x = LON, y = LAT, color = D))
## Using lon and lat as longitude and latitude columns, respectively.
## projection transformed from +init=epsg:4326 to +proj=stere +lat_0=90 +lat_ts=71 +lon_0=-61 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
REGION3
## Assuming `crs = 4326` in stat_spatial_identity()
# Statistics Code
# TO DO: Figure out a way to determine even spatial distribution of profiles by year.
TotalCasts = length(unique(OMGALL$CAST))
CastDistALL = as.numeric(format(OMGALL$DATE[BTMALL_idx],'%Y'))
CastHistYEAR = ggplot(BTMALL,aes(x = DATE)) +
geom_histogram(bins = 9, color="darkblue", fill="lightblue") +
xlab("Year") +
ylab("Number of Casts") +
ggtitle("Annual Distribution of Casts in Region 3","Total Number of Casts = 310") +
theme(plot.title = element_text(hjust = 0.5))
CastHistYEAR
dt = BTMALL
dt$DATE = as.numeric(format(dt$DATE,'%m'))
CastHistMONTH = ggplot(dt,aes(x = DATE)) +
geom_histogram(bins = 30, color="darkblue", fill="lightblue") +
xlab("Month") +
ylab("Number of Casts") +
ggtitle("Seasonal Distribution of Casts in Region 3 (2017 - 2021)","Total Number of Casts = 310") +
theme(plot.title = element_text(hjust = 0.5))+
xlim(0.5,11.5)
CastHistMONTH
## Warning: Removed 2 rows containing missing values (geom_bar).
BTMALL_vec = OMGALL$D[BTMALL_idx]
DeepDistAll = ggplot(NULL) +
geom_histogram(bins = 40, color="darkblue", fill="lightblue") +
aes(BTMALL_vec) +
xlab("Deepest Depth of CTD") +
ylab("Number of Casts") +
ggtitle("Max Depth Distribution of all Casts","Total Number of Casts = 310") +
theme(plot.title = element_text(hjust = 0.5))
DeepDistAll
# Presets for following code.
temp_crit = 0.5 # Obata et al. Temperature threshold for MLD_T
pden_crit = 0.125 # Miller (1976) Density threshold for MLD_D
s_ref = 34.8 # Spall (2013) Salinity Reference for Freshwater Calculation
# Single Profile Plots.
# All values are re-run so individual tests can be run without adjusting values for every cast in the previous loops. Use for in-depth analysis of a single profile.
# TEST is the cast number from OMGALL for the desired profiles.
# NOTE: Cast 31.2 T and S sensor turned on late. Not good for MLD.
# Cast 202.4 T and S sensor intermittent. Not good for MLD.
TEST = 186.1
rm(data,Depth,Temp,Salt,maxDepth,binDepth,binSalt,binTemp,binDepth,binPTemp,binSA,binPres,binCT,binPDen,data_p,data_SA,data_pt,data_CT,data_r0,data,r1,data_N2,data_N2_min,SST,SSPDen,MLD_T_1,MLD_T_2,MLD_D_1,MLD_D_2,a)
data = filter(OMGALL,CAST == TEST)
Depth <- as.numeric(data[,5])
Temp <- as.numeric(data[,6])
Salt <- as.numeric(data[,7])
maxDepth = as.integer(tail (Depth,1))
binDepth = seq(2,maxDepth,1)
binSalt = binMeans(Salt,x=Depth,bx=binDepth,na.rm=TRUE)
binTemp = binMeans(Temp,x=Depth,bx=binDepth,na.rm=TRUE)
binDepth = head(binDepth,-1)
binPTemp = gsw_pt_from_t(binSalt,binTemp,binDepth,0)
binSA = gsw_SA_from_SP(binSalt, binDepth, data$LON, data$LAT)
binPres = gsw_p_from_z(binDepth*-1,latitude = data$LAT)
binCT = gsw_CT_from_t(binSA,binTemp,binPres)
binPDen = gsw_sigma0(binSA,binCT)
data_p = gsw_p_from_z(data$D*-1,latitude = data$LAT)
data_SA = gsw_SA_from_SP(data$S, data_p, data$LON, data$LAT)
data_pt = gsw_pt_from_t(data_SA, data$T, data_p,0)
data_CT = gsw_CT_from_t(data_SA, data$T,data_p)
data_r0 = gsw_sigma0(data_SA, data_CT)
data_r1 = gsw_sigma1(data_SA, data_CT)
data_N2 = gsw_Nsquared(data_SA,data_CT,data_p,latitude = data$LAT)
data_N2_min = data$D[which(data_N2$N2 == min(data_N2$N2, na.rm = TRUE))]
SST = binTemp[1]
SSPDen = binPDen[1]
a = binDepth[which(binTemp <= (SST - temp_crit))]
MLD_T_TEST_1 = a[1]
a = binDepth[which(binTemp >= (SST + temp_crit))]
MLD_T_TEST_2 = a[1]
MLD_T_TEST = min(c(MLD_T_TEST_1,MLD_T_TEST_2),na.rm = TRUE)
if (MLD_T_TEST == Inf) {
MLD_T_TEST = NA
}
if (is.na(MLD_T_TEST) == FALSE) {
MLD_T_TEST_Ptemp_Avg = mean(binPTemp[1:MLD_T_TEST])
MLD_T_TEST_Ptemp_Sd = sd(binPTemp[1:MLD_T_TEST])
MLD_T_TEST_SA_Avg = mean(binSA[1:MLD_T_TEST])
MLD_T_TEST_SA_Sd = sd(binSA[1:MLD_T_TEST])
MLD_T_TEST_Pden_Avg = mean(binPDen[1:MLD_T_TEST])
MLD_T_TEST_Pden_Sd = sd(binPDen[1:MLD_T_TEST])
}
if (is.na(MLD_T_TEST) == TRUE) {
MLD_T_TEST_Ptemp_Avg = NA
MLD_T_TEST_Ptemp_Sd = NA
MLD_T_TEST_SA_Avg = NA
MLD_T_TEST_SA_Sd = NA
MLD_D_TEST_Pden_Avg = NA
MLD_D_TEST_Pden_Sd = NA
}
a = binDepth[which(binPDen <= SSPDen*pden_crit - SSPDen)]
MLD_D_TEST_1 = a[1]
a = binDepth[which(binPDen >= SSPDen*pden_crit + SSPDen)]
MLD_D_TEST_2 = a[1]
MLD_D_TEST = min(c(MLD_D_TEST_1,MLD_D_TEST_2), na.rm = TRUE)
if (MLD_D_TEST == Inf) {
MLD_D_TEST = NA
}
if (is.na(MLD_D_TEST) == FALSE) {
MLD_D_TEST_Ptemp_Avg = mean(binPTemp[1:MLD_D_TEST])
MLD_D_TEST_Ptemp_Sd = sd(binPTemp[1:MLD_D_TEST])
MLD_D_TEST_SA_Avg = mean(binSA[1:MLD_D_TEST])
MLD_D_TEST_SA_Sd = sd(binSA[1:MLD_D_TEST])
MLD_D_TEST_Pden_Avg = mean(binPDen[1:MLD_D_TEST])
MLD_D_TEST_Pden_Sd = sd(binPDen[1:MLD_D_TEST])
}
if (is.na(MLD_D_TEST) == TRUE) {
MLD_D_TEST_Ptemp_Avg = NA
MLD_D_TEST_Ptemp_Sd = NA
MLD_D_TEST_SA_Avg = NA
MLD_D_TEST_SA_Sd = NA
MLD_D_TEST_Pden_Avg = NA
MLD_D_TEST_Pden_Sd = NA
}
tmax_TEST = max(binPTemp[(max(data_N2$N2,na.rm = TRUE)):length(binPTemp)],na.rm = TRUE)
tmin_TEST = min(binPTemp,na.rm=TRUE)
dtmax_TEST = max(binPTemp[na.omit(binDepth[binPTemp == tmin_TEST]):length(binDepth)],na.rm=TRUE)
tmax_depth_TEST = binDepth[which(binPTemp == tmax_TEST)]
tmax_SA_TEST = data_SA[tmax_depth_TEST]
tmin_depth_TEST = binDepth[which(binPTemp == tmin_TEST)]
tmin_SA_TEST = data_SA[tmin_depth_TEST]
dtmax_depth_TEST = binDepth[which(binPTemp == dtmax_TEST)]
dtmax_SA_TEST = data_SA[dtmax_depth_TEST]
ggplot(data,aes(x=T,y=-D)) +
geom_path() +
ggtitle("Temperature Profile",data$CAST) +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Temp, C") +
ylab("Depth, m")
ggplot(data,aes(x=S,y=-D)) +
geom_path() +
ggtitle("Salinity Profile",data$CAST) +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Salinity, psu") +
ylab("Depth, m")
ggplot(data,mapping=aes(x=S,y=T,color=D)) +
geom_point() +
ggtitle("T(S) Profile",data$CAST) +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Salinity, psu") +
ylab("Temperature, C") +
scale_color_viridis_c(guide=guide_colorbar(barheight=15))
ggplot(NULL) +
geom_path(data = data, aes(x = T, color = "Insitu Temp")) +
geom_path(data = NULL, aes(x = data_pt, color = "Potential Temp")) +
aes(y = data$D*-1) +
ggtitle("In-situ Temp vs Potential Temperature", data$CAST) +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("T and PT, C") +
ylab("Depth, m") +
labs(colour = "Unit")
ggplot(NULL) +
geom_path(data = NULL, aes(x = data_SA, y = data_r0)) +
ggtitle("Absolute Salinity and Potential Density Profile", data$CAST) +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("SA g/kg") +
ylab("Potential Density kg/m^3, ref = 0")
ggplot(NULL) +
geom_path(data = NULL, aes(x=binPTemp,y=-binDepth)) +
ggtitle("Potential Temperature Profile with MLD_T",data$CAST) +
annotation_spatial_hline(intercept = -MLD_T_TEST, color="red") +
annotation_spatial_vline(intercept = MLD_T_TEST_Ptemp_Avg,color = "blue") +
annotation_spatial_hline(intercept = -data_N2_min, color = "green") +
ylim(-(MLD_T_TEST + 50), 0)+
xlim(-2, 6) +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Potential Temperature [C]") +
ylab("Depth, [m]") +
labs(caption = "Red = MLD from Temperature Threshold, Green = N2 Minimum, Blue = Average Potential Temperature above the MLD")
## Assuming `crs = 4326` in annotation_spatial_(h|v)line()
## Assuming `crs = 4326` in annotation_spatial_(h|v)line()
## Assuming `crs = 4326` in annotation_spatial_(h|v)line()
ggplot(NULL) +
geom_path(data = NULL,aes(x=binSA,y=-binDepth)) +
ggtitle("Absolute Salinity Profile with MLD_T",data$CAST) +
annotation_spatial_hline(intercept = -MLD_T_TEST, color = "red") +
annotation_spatial_vline(intercept = MLD_T_TEST_SA_Avg,color = "blue") +
annotation_spatial_hline(intercept = -data_N2_min, color = "green")+
ylim(-(MLD_T_TEST + 50), 0)+
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Absolute Salinity [g/kg]") +
ylab("Depth [m]") +
labs(caption = "Red = MLD from Temperature Threshold, Green = N2 Minimum, Blue = Average Absolute Salinity above the MLD")
## Assuming `crs = 4326` in annotation_spatial_(h|v)line()
## Assuming `crs = 4326` in annotation_spatial_(h|v)line()
## Assuming `crs = 4326` in annotation_spatial_(h|v)line()
ggplot(NULL) +
geom_path(data = NULL,aes(x=binPDen,y=-binDepth)) +
ggtitle("Potential Density Profile with MLD_T",data$CAST) +
annotation_spatial_hline(intercept = -MLD_T_TEST, color = "red") +
annotation_spatial_vline(intercept = MLD_T_TEST_Pden_Avg,color = "blue") +
annotation_spatial_hline(intercept = -data_N2_min, color = "green") +
ylim(-(MLD_T_TEST + 50), 0)+
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Potential Density [kg/m^3]") +
ylab("Depth [m]") +
labs(caption = "Red = MLD from Temperature Threshold, Green = N2 Minimum, Blue = Average Potential Density above the MLD")
## Assuming `crs = 4326` in annotation_spatial_(h|v)line()
## Assuming `crs = 4326` in annotation_spatial_(h|v)line()
## Assuming `crs = 4326` in annotation_spatial_(h|v)line()
ggplot(NULL) +
geom_path(data = NULL, aes(x=binPTemp,y=-binDepth)) +
ggtitle("Potential Temperature Profile with MLD_D",data$CAST) +
annotation_spatial_hline(intercept = -MLD_D_TEST, color = "red") +
annotation_spatial_vline(intercept = MLD_D_TEST_Ptemp_Avg,color = "blue") +
annotation_spatial_hline(intercept = -data_N2_min, color = "green") +
ylim(-(MLD_D_TEST + 50), 0)+
xlim(-2, 6) +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Potential Temperature [C]") +
ylab("Depth, [m]") +
labs(caption = "Red = MLD from Density Threshold, Green = N2 Minimum, Blue = Average Potential Temperature above the MLD")
## Assuming `crs = 4326` in annotation_spatial_(h|v)line()
## Assuming `crs = 4326` in annotation_spatial_(h|v)line()
## Assuming `crs = 4326` in annotation_spatial_(h|v)line()
ggplot(NULL) +
geom_path(data = NULL,aes(x=binSA,y=-binDepth)) +
ggtitle("Absolute Salinity Profile with MLD_D",data$CAST) +
annotation_spatial_hline(intercept = -MLD_D_TEST, color = "red") +
annotation_spatial_vline(intercept = MLD_D_TEST_SA_Avg,color = "blue") +
annotation_spatial_hline(intercept = -data_N2_min, color = "green") +
ylim(-(MLD_D_TEST + 50), 0)+
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Absolute Salinity [g/kg]") +
ylab("Depth [m]") +
labs(caption = "Red = MLD from Density Threshold, Green = N2 Minimum, Blue = Average Absolute Salinity above the MLD")
## Assuming `crs = 4326` in annotation_spatial_(h|v)line()
## Assuming `crs = 4326` in annotation_spatial_(h|v)line()
## Assuming `crs = 4326` in annotation_spatial_(h|v)line()
ggplot(NULL) +
geom_path(data = NULL,aes(x=binPDen,y=-binDepth)) +
ggtitle("Potential Density Profile with MLD_D",data$CAST) +
annotation_spatial_hline(intercept = -MLD_D_TEST, color = "red") +
annotation_spatial_vline(intercept = MLD_D_TEST_Pden_Avg,color = "blue") +
annotation_spatial_hline(intercept = -data_N2_min, color = "green") +
ylim(-(MLD_D_TEST + 50), 0)+
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Potential Density [kg/m^3]") +
ylab("Depth [m]") +
labs(caption = "Red = MLD from Density Threshold, Green = N2 Minimum, Blue = Average Potential Density above the MLD")
## Assuming `crs = 4326` in annotation_spatial_(h|v)line()
## Assuming `crs = 4326` in annotation_spatial_(h|v)line()
## Assuming `crs = 4326` in annotation_spatial_(h|v)line()
ggplot(NULL) +
geom_path(data = NULL, aes(x = data_N2$N2[1:length(which(data$D < 200))], y = -data$D[1:length(which(data$D < 200))])) +
ggtitle("Near Surface (<200m) N2",data$CAST) +
theme(plot.title = element_text(hjust = 0.5)) +
ylab("Depth [m]") +
xlab("N2")
#rm(data,Depth,Temp,Salt,maxDepth,binDepth,binSalt,binTemp,binDepth,binPTemp,binSA,binPres,binCT,binPDen,data_p,data_SA,data_pt,data_CT,data_r0,data,r1,data_N2,data_N2_min,SST,SSPDen,MLD_T_1,MLD_T_2,MLD_D_1,MLD_D_2,a)
# Mixed Layer Depth And Fresh Water Loop.
# For loop that calculates the MLD and MLD properties from both a temperature and density threshold method. Also measures freshwater content. Then updates the results to the BTMALL data.frame and splits the data.frame by year.
# Additional values can be recorded with this structure by adding a blank vector before the loop, adding '[i]' to the value to be recorded, and removing the variable from the rm() command.
MLD_T = vector()
MLD_T_Ptemp_Avg = vector()
MLD_T_Ptemp_Sd = vector()
MLD_T_SA_Avg = vector()
MLD_T_SA_Sd = vector()
MLD_T_Pden_Avg = vector()
MLD_T_Pden_Sd = vector()
MLD_D = vector()
MLD_D_Ptemp_Avg = vector()
MLD_D_Ptemp_Sd = vector()
MLD_D_SA_Avg = vector()
MLD_D_SA_Sd = vector()
MLD_D_Pden_Avg = vector()
MLD_D_Pden_Sd = vector()
N2_min = vector()
N2_max = vector()
tmax = vector()
tmin = vector()
dtmax = vector()
tmax_depth = vector()
tmax_SA = vector()
tmin_depth = vector()
tmin_SA = vector()
dtmax_depth = vector()
dtmax_SA = vector()
fresh = vector()
for (i in BTMALL_idx) {
rm(data,Depth,Temp,Salt,maxDepth,binDepth,binSalt,binTemp,binDepth,binPTemp,binSA,binPres,binCT,binPDen,data_p,data_SA,data_pt,data_CT,data_r0,data,r1,data_N2,data_N2_min,SST,SSPDen,MLD_T_1,MLD_T_2,MLD_D_1,MLD_D_2,a)
data = filter(OMGALL,CAST == OMGALL$CAST[i])
Depth <- as.numeric(data[,5])
Temp <- as.numeric(data[,6])
Salt <- as.numeric(data[,7])
maxDepth = as.integer(tail (Depth,1))
binDepth = seq(2,maxDepth,1)
binSalt = binMeans(Salt,x=Depth,bx=binDepth,na.rm=TRUE)
binTemp = binMeans(Temp,x=Depth,bx=binDepth,na.rm=TRUE)
binDepth = head(binDepth,-1)
binPTemp = gsw_pt_from_t(binSalt,binTemp,binDepth,0)
binSA = gsw_SA_from_SP(binSalt, binDepth, data$LON, data$LAT)
binPres = gsw_p_from_z(binDepth*-1,latitude = data$LAT)
binCT = gsw_CT_from_t(binSA,binTemp,binPres)
binPDen = gsw_sigma0(binSA,binCT)
data_p = gsw_p_from_z(data$D*-1,latitude = data$LAT)
data_SA = gsw_SA_from_SP(data$S, data_p, data$LON, data$LAT)
data_pt = gsw_pt_from_t(data_SA, data$T, data_p,0)
data_CT = gsw_CT_from_t(data_SA, data$T,data_p)
data_r0 = gsw_sigma0(data_SA, data_CT)
data_r1 = gsw_sigma1(data_SA, data_CT)
data_N2 = gsw_Nsquared(data_SA,data_CT,data_p,latitude = data$LAT)
data_N2_min = data$D[which(data_N2$N2 == min(data_N2$N2, na.rm = TRUE))]
data_N2_max = data$D[which(data_N2$N2 == max(data_N2$N2, na.rm = TRUE))]
SST = binTemp[1]
SSPDen = binPDen[1]
a = binDepth[which(binTemp <= (SST - temp_crit))]
MLD_T_1 = a[1]
a = binDepth[which(binTemp >= (SST + temp_crit))]
MLD_T_2 = a[1]
MLD_T[i] = min(c(MLD_T_1,MLD_T_2),na.rm = TRUE)
if (MLD_T[i] == Inf) {
MLD_T[i] = NA
}
if (is.na(MLD_T[i]) == FALSE) {
MLD_T_Ptemp_Avg[i] = mean(binPTemp[1:MLD_T[i]])
MLD_T_Ptemp_Sd[i] = sd(binPTemp[1:MLD_T[i]])
MLD_T_SA_Avg[i] = mean(binSA[1:MLD_T[i]])
MLD_T_SA_Sd[i] = sd(binSA[1:MLD_T[i]])
MLD_T_Pden_Avg[i] = mean(binPDen[1:MLD_T[i]])
MLD_T_Pden_Sd[i] = sd(binPDen[1:MLD_T[i]])
}
if (is.na(MLD_T[i]) == TRUE) {
MLD_T_Ptemp_Avg[i] = NA
MLD_T_Ptemp_Sd[i] = NA
MLD_T_SA_Avg[i] = NA
MLD_T_SA_Sd[i] = NA
MLD_T_Pden_Avg[i] = NA
MLD_T_Pden_Sd[i] = NA
}
a = binDepth[which(binPDen <= SSPDen*pden_crit - SSPDen)]
MLD_D_1 = a[1]
a = binDepth[which(binPDen >= SSPDen*pden_crit + SSPDen)]
MLD_D_2 = a[1]
MLD_D[i] = min(c(MLD_D_1,MLD_D_2), na.rm = TRUE)
if (MLD_D[i] == Inf) {
MLD_D[i] = NA
}
if (is.na(MLD_D[i]) == FALSE) {
MLD_D_Ptemp_Avg[i] = mean(binPTemp[1:MLD_D[i]])
MLD_D_Ptemp_Sd[i] = sd(binPTemp[1:MLD_D[i]])
MLD_D_SA_Avg[i] = mean(binSA[1:MLD_D[i]])
MLD_D_SA_Sd[i] = sd(binSA[1:MLD_D[i]])
MLD_D_Pden_Avg[i] = mean(binPDen[1:MLD_D[i]])
MLD_D_Pden_Sd[i] = sd(binPDen[1:MLD_D[i]])
}
if (is.na(MLD_D[i]) == TRUE) {
MLD_D_Ptemp_Avg[i] = NA
MLD_D_Ptemp_Sd[i] = NA
MLD_D_SA_Avg[i] = NA
MLD_D_SA_Sd[i] = NA
MLD_D_Pden_Avg[i] = NA
MLD_D_Pden_Sd[i] = NA
}
tmax[i] = max(binPTemp[data_N2_max:length(binPTemp)],na.rm = TRUE)
tmin[i] = min(binPTemp,na.rm = TRUE)
dtmax[i] = max(binPTemp[na.omit(binDepth[binPTemp == tmin[i]]):length(binDepth)],na.rm=TRUE)
tmax_depth[i] = binDepth[which(binPTemp == tmax[i])]
tmax_SA[i] = data_SA[tmax_depth[i]]
tmin_depth[i] = binDepth[which(binPTemp == tmin[i])]
tmin_SA[i] = data_SA[tmin_depth[i]]
dtmax_depth[i] = binDepth[which(binPTemp == dtmax[i])]
dtmax_SA[i] = data_SA[dtmax_depth[i]]
N2_min[i] = data_N2_min
N2_max[i] = data_N2_max
s_ref_vec = c(rep(s_ref,length(binDepth)))
fresh[i] = round(sum(((s_ref_vec - binSalt)/s_ref_vec),na.rm = TRUE),digits=4)
}
MLD_T_Ptemp_Sd = MLD_T_Ptemp_Sd[BTMALL_idx]
MLD_T_SA_Sd = MLD_T_SA_Sd[BTMALL_idx]
MLD_T_Pden_Sd = MLD_T_Pden_Sd[BTMALL_idx]
MLD_D_Ptemp_Sd = MLD_D_Ptemp_Sd[BTMALL_idx]
MLD_D_Pden_Sd = MLD_D_Pden_Sd[BTMALL_idx]
N2_min = N2_min[BTMALL_idx]
N2_max = N2_max[BTMALL_idx]
BTMALL$MLD_T = MLD_T[BTMALL_idx]
BTMALL$MLD_T_Ptemp_Avg = MLD_T_Ptemp_Avg[BTMALL_idx]
BTMALL$MLD_T_SA_Avg = MLD_T_SA_Avg[BTMALL_idx]
BTMALL$MLD_T_Pden_Avg = MLD_T_Pden_Avg[BTMALL_idx]
BTMALL$MLD_D = MLD_D[BTMALL_idx]
BTMALL$MLD_D_Ptemp_Avg = MLD_D_Ptemp_Avg[BTMALL_idx]
BTMALL$MLD_D_SA_Avg = MLD_D_SA_Avg[BTMALL_idx]
BTMALL$MLD_D_Pden_Avg = MLD_D_Pden_Avg[BTMALL_idx]
BTMALL$TMAX = tmax[BTMALL_idx]
BTMALL$TMAXD = tmax_depth[BTMALL_idx]
BTMALL$TMAXSA = tmax_SA[BTMALL_idx]
BTMALL$TMIN = tmin[BTMALL_idx]
BTMALL$TMIND = tmin_depth[BTMALL_idx]
BTMALL$TMINSA = tmin_SA[BTMALL_idx]
BTMALL$DTMAX = dtmax[BTMALL_idx]
BTMALL$DTMAXD = dtmax_depth[BTMALL_idx]
BTMALL$DTMAXSA = dtmax_SA[BTMALL_idx]
BTMALL$FRESH = fresh[BTMALL_idx]
BTM2021 = BTMALL[BTMALL$DATE > "2021-01-01" & BTMALL$DATE < "2021-12-31",]
BTM2020 = BTMALL[BTMALL$DATE > "2020-01-01" & BTMALL$DATE < "2020-12-31",]
BTM2019 = BTMALL[BTMALL$DATE > "2019-01-01" & BTMALL$DATE < "2019-12-31",]
BTM2018 = BTMALL[BTMALL$DATE > "2018-01-01" & BTMALL$DATE < "2018-12-31",]
BTM2017 = BTMALL[BTMALL$DATE > "2017-01-01" & BTMALL$DATE < "2017-12-31",]
# MLD Analysis
# DISCUSSION: I decided to calculate MLD from both a temperature threshold and a density threshold. From my observations, the MLD_T is a much better indication of the near surface mixed layer present in my region. I wanted to retain the MLD_D though, as sigma*0.125 is derived from the subtropical mode water of the North Atlantic (Levitus 1982). I figured this value could be useful for analyzing changes in the intruding Atlantic water as any value below MLD_D should be water primarily from this source.
# I should note, MLD_D is not available for about half of the profiles (154/310). This is not a bug, the reason is that there is no value at which pden == pden[1]*1.125 for these profiles. This could be a finding as well, as I think that it means that Atlantic water is not present, or at least not significantly strong in the region. This could also indicate that the profiles are very shallow. I intend to plot these unique profiles on a map to find correlation. As the MLD_D is incomplete and not indicative of the near surface mixed layer in my region, I will only be plotting MLD property trends from MLD_T data.
# UPDATE TO DISCUSSION: I find no obvious visual correlation in the spatial distribution or the depth distribution for the casts with no MLD_D value. Perhaps this is just a bad definition for MLD in my region. Or perhaps the water in region is not prominently from the subtropical mode water.
MLD_T_Avg_2021 = mean(BTM2021$MLD_T,na.rm = TRUE)
MLD_T_Avg_2020 = mean(BTM2020$MLD_T,na.rm = TRUE)
MLD_T_Avg_2019 = mean(BTM2019$MLD_T,na.rm = TRUE)
MLD_T_Avg_2018 = mean(BTM2018$MLD_T,na.rm = TRUE)
MLD_T_Avg_2017 = mean(BTM2017$MLD_T,na.rm = TRUE)
MLD_T_Sd_2021 = sd(BTM2021$MLD_T,na.rm = TRUE)
MLD_T_Sd_2020 = sd(BTM2020$MLD_T,na.rm = TRUE)
MLD_T_Sd_2019 = sd(BTM2019$MLD_T,na.rm = TRUE)
MLD_T_Sd_2018 = sd(BTM2018$MLD_T,na.rm = TRUE)
MLD_T_Sd_2017 = sd(BTM2017$MLD_T,na.rm = TRUE)
MLD_T_Med_2021 = median(BTM2021$MLD_T,na.rm = TRUE)
MLD_T_Med_2020 = median(BTM2020$MLD_T,na.rm = TRUE)
MLD_T_Med_2019 = median(BTM2019$MLD_T,na.rm = TRUE)
MLD_T_Med_2018 = median(BTM2018$MLD_T,na.rm = TRUE)
MLD_T_Med_2017 = median(BTM2017$MLD_T,na.rm = TRUE)
date_vec = c(mean(BTM2017$DATE),mean(BTM2018$DATE),mean(BTM2019$DATE),mean(BTM2020$DATE),mean(BTM2021$DATE))
MLD_T_Avg_Vec = c(MLD_T_Avg_2017,MLD_T_Avg_2018,MLD_T_Avg_2019,MLD_T_Avg_2020,MLD_T_Avg_2021)
MLD_T_Med_Vec = c(MLD_T_Med_2017,MLD_T_Med_2018,MLD_T_Med_2019,MLD_T_Med_2020,MLD_T_Med_2021)
MLD_T_Sd_Vec = c(MLD_T_Sd_2017,MLD_T_Sd_2018,MLD_T_Sd_2019,MLD_T_Sd_2020,MLD_T_Sd_2021)
MLD_T_DATA = data.frame(MLD_T_Avg_Vec,MLD_T_Med_Vec,MLD_T_Sd_Vec)
MLD_T_Ptemp_Avg_2021 = mean(BTM2021$MLD_T_Ptemp_Avg,na.rm = TRUE)
MLD_T_Ptemp_Avg_2020 = mean(BTM2020$MLD_T_Ptemp_Avg,na.rm = TRUE)
MLD_T_Ptemp_Avg_2019 = mean(BTM2019$MLD_T_Ptemp_Avg,na.rm = TRUE)
MLD_T_Ptemp_Avg_2018 = mean(BTM2018$MLD_T_Ptemp_Avg,na.rm = TRUE)
MLD_T_Ptemp_Avg_2017 = mean(BTM2017$MLD_T_Ptemp_Avg,na.rm = TRUE)
MLD_T_Ptemp_Sd_2021 = sd(BTM2021$MLD_T_Ptemp_Avg,na.rm = TRUE)
MLD_T_Ptemp_Sd_2020 = sd(BTM2020$MLD_T_Ptemp_Avg,na.rm = TRUE)
MLD_T_Ptemp_Sd_2019 = sd(BTM2019$MLD_T_Ptemp_Avg,na.rm = TRUE)
MLD_T_Ptemp_Sd_2018 = sd(BTM2018$MLD_T_Ptemp_Avg,na.rm = TRUE)
MLD_T_Ptemp_Sd_2017 = sd(BTM2017$MLD_T_Ptemp_Avg,na.rm = TRUE)
MLD_T_Ptemp_Med_2021 = median(BTM2021$MLD_T_Ptemp_Avg,na.rm = TRUE)
MLD_T_Ptemp_Med_2020 = median(BTM2020$MLD_T_Ptemp_Avg,na.rm = TRUE)
MLD_T_Ptemp_Med_2019 = median(BTM2019$MLD_T_Ptemp_Avg,na.rm = TRUE)
MLD_T_Ptemp_Med_2018 = median(BTM2018$MLD_T_Ptemp_Avg,na.rm = TRUE)
MLD_T_Ptemp_Med_2017 = median(BTM2017$MLD_T_Ptemp_Avg,na.rm = TRUE)
MLD_T_Ptemp_Avg_Vec = c(MLD_T_Ptemp_Avg_2017,MLD_T_Ptemp_Avg_2018,MLD_T_Ptemp_Avg_2019,MLD_T_Ptemp_Avg_2020,MLD_T_Ptemp_Avg_2021)
MLD_T_Ptemp_Sd_Vec = c(MLD_T_Ptemp_Sd_2017,MLD_T_Ptemp_Sd_2018,MLD_T_Ptemp_Sd_2019,MLD_T_Ptemp_Sd_2020,MLD_T_Ptemp_Sd_2021)
MLD_T_Ptemp_Med_Vec = c(MLD_T_Ptemp_Med_2017,MLD_T_Ptemp_Med_2018,MLD_T_Ptemp_Med_2019,MLD_T_Ptemp_Med_2020,MLD_T_Ptemp_Med_2021)
MLD_T_SA_Avg_2021 = mean(BTM2021$MLD_T_SA_Avg,na.rm = TRUE)
MLD_T_SA_Avg_2020 = mean(BTM2020$MLD_T_SA_Avg,na.rm = TRUE)
MLD_T_SA_Avg_2019 = mean(BTM2019$MLD_T_SA_Avg,na.rm = TRUE)
MLD_T_SA_Avg_2018 = mean(BTM2018$MLD_T_SA_Avg,na.rm = TRUE)
MLD_T_SA_Avg_2017 = mean(BTM2017$MLD_T_SA_Avg,na.rm = TRUE)
MLD_T_SA_Sd_2021 = sd(BTM2021$MLD_T_SA_Avg,na.rm = TRUE)
MLD_T_SA_Sd_2020 = sd(BTM2020$MLD_T_SA_Avg,na.rm = TRUE)
MLD_T_SA_Sd_2019 = sd(BTM2019$MLD_T_SA_Avg,na.rm = TRUE)
MLD_T_SA_Sd_2018 = sd(BTM2018$MLD_T_SA_Avg,na.rm = TRUE)
MLD_T_SA_Sd_2017 = sd(BTM2017$MLD_T_SA_Avg,na.rm = TRUE)
MLD_T_SA_Med_2021 = median(BTM2021$MLD_T_SA_Avg,na.rm = TRUE)
MLD_T_SA_Med_2020 = median(BTM2020$MLD_T_SA_Avg,na.rm = TRUE)
MLD_T_SA_Med_2019 = median(BTM2019$MLD_T_SA_Avg,na.rm = TRUE)
MLD_T_SA_Med_2018 = median(BTM2018$MLD_T_SA_Avg,na.rm = TRUE)
MLD_T_SA_Med_2017 = median(BTM2017$MLD_T_SA_Avg,na.rm = TRUE)
MLD_T_SA_Avg_Vec = c(MLD_T_SA_Avg_2017,MLD_T_SA_Avg_2018,MLD_T_SA_Avg_2019,MLD_T_SA_Avg_2020,MLD_T_SA_Avg_2021)
MLD_T_SA_Sd_Vec = c(MLD_T_SA_Sd_2017,MLD_T_SA_Sd_2018,MLD_T_SA_Sd_2019,MLD_T_SA_Sd_2020,MLD_T_SA_Sd_2021)
MLD_T_SA_Med_Vec = c(MLD_T_SA_Med_2017,MLD_T_SA_Med_2018,MLD_T_SA_Med_2019,MLD_T_SA_Med_2020,MLD_T_SA_Med_2021)
MLD_T_Pden_Avg_2021 = mean(BTM2021$MLD_T_Pden_Avg,na.rm = TRUE)
MLD_T_Pden_Avg_2020 = mean(BTM2020$MLD_T_Pden_Avg,na.rm = TRUE)
MLD_T_Pden_Avg_2019 = mean(BTM2019$MLD_T_Pden_Avg,na.rm = TRUE)
MLD_T_Pden_Avg_2018 = mean(BTM2018$MLD_T_Pden_Avg,na.rm = TRUE)
MLD_T_Pden_Avg_2017 = mean(BTM2017$MLD_T_Pden_Avg,na.rm = TRUE)
MLD_T_Pden_Sd_2021 = sd(BTM2021$MLD_T_Pden_Avg,na.rm = TRUE)
MLD_T_Pden_Sd_2020 = sd(BTM2020$MLD_T_Pden_Avg,na.rm = TRUE)
MLD_T_Pden_Sd_2019 = sd(BTM2019$MLD_T_Pden_Avg,na.rm = TRUE)
MLD_T_Pden_Sd_2018 = sd(BTM2018$MLD_T_Pden_Avg,na.rm = TRUE)
MLD_T_Pden_Sd_2017 = sd(BTM2017$MLD_T_Pden_Avg,na.rm = TRUE)
MLD_T_Pden_Med_2021 = median(BTM2021$MLD_T_Pden_Avg,na.rm = TRUE)
MLD_T_Pden_Med_2020 = median(BTM2020$MLD_T_Pden_Avg,na.rm = TRUE)
MLD_T_Pden_Med_2019 = median(BTM2019$MLD_T_Pden_Avg,na.rm = TRUE)
MLD_T_Pden_Med_2018 = median(BTM2018$MLD_T_Pden_Avg,na.rm = TRUE)
MLD_T_Pden_Med_2017 = median(BTM2017$MLD_T_Pden_Avg,na.rm = TRUE)
MLD_T_Pden_Avg_Vec = c(MLD_T_Pden_Avg_2017,MLD_T_Pden_Avg_2018,MLD_T_Pden_Avg_2019,MLD_T_Pden_Avg_2020,MLD_T_Pden_Avg_2021)
MLD_T_Pden_Sd_Vec = c(MLD_T_Pden_Sd_2017,MLD_T_Pden_Sd_2018,MLD_T_Pden_Sd_2019,MLD_T_Pden_Sd_2020,MLD_T_Pden_Sd_2021)
MLD_T_Pden_Med_Vec = c(MLD_T_Pden_Med_2017,MLD_T_Pden_Med_2018,MLD_T_Pden_Med_2019,MLD_T_Pden_Med_2020,MLD_T_Pden_Med_2021)
ggplot(NULL) +
geom_point(data = BTMALL, aes(x = MLD_T, y = MLD_T_Ptemp_Avg)) +
ggtitle("Scatter Plot of MLD_T vs MLD_T_Ptemp_Avg") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Mixed Layer Depth (from Temperature)") +
ylab("Average Potential Temperature above the MLD [C]")
ggplot(NULL) +
geom_point(data = BTMALL, aes(x = MLD_T, y = MLD_T_SA_Avg)) +
ggtitle("Scatter Plot of MLD_T vs MLD_T_SA_Avg") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Mixed Layer Depth (from Temperature)") +
ylab("Average Absolute Salinity above the MLD [g/kg]")
ggplot(NULL) +
geom_point(data = BTMALL, aes(x = MLD_T, y = MLD_T_Pden_Avg)) +
ggtitle("Scatter Plot of MLD_T vs MLD_T_Pden_Avg") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Mixed Layer Depth (from Temperature)") +
ylab("Average Potential Density above the MLD [kg/m^3]")
ggplot(NULL) +
geom_point(data = BTMALL, aes(x = MLD_D, y = MLD_D_Ptemp_Avg)) +
ggtitle("Scatter Plot of MLD_D vs MLD_D_Ptemp_Avg") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Mixed Layer Depth (from Density)") +
ylab("Average Potential Temperature above the MLD [C]")
ggplot(NULL) +
geom_point(data = BTMALL, aes(x = MLD_D, y = MLD_D_SA_Avg)) +
ggtitle("Scatter Plot of MLD_D vs MLD_D_SA_Avg") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Mixed Layer Depth (from Density)") +
ylab("Average Absolute Salinity above the MLD [g/kg]")
ggplot(NULL) +
geom_point(data = BTMALL, aes(x = MLD_D, y = MLD_D_Pden_Avg)) +
ggtitle("Scatter Plot of MLD_D vs MLD_D_Pden_Avg") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Mixed Layer Depth (from Density)") +
ylab("Average Potential Density above the MLD [kg/m^3]")
ggplot(NULL, aes(x = date_vec, y = MLD_T_Avg_Vec)) +
geom_path(data = NULL,aes(x = date_vec, y = MLD_T_Avg_Vec), color = "red") +
geom_point(data = NULL,aes(x = date_vec, y = MLD_T_Avg_Vec), color = "red", size = 4) +
geom_path(data = NULL,aes(x = date_vec, y = MLD_T_Med_Vec), color = "green") +
geom_point(data = NULL,aes(x = date_vec, y = MLD_T_Med_Vec), color = "green", size = 4) +
geom_point(data = BTMALL,aes(x=BTMALL$DATE,y=BTMALL$MLD_T),alpha = 0.1) +
geom_errorbar(aes(ymin = MLD_T_Avg_Vec - MLD_T_Sd_Vec, ymax = MLD_T_Avg_Vec + MLD_T_Sd_Vec),width = 50) +
ggtitle("Mixed Layer Depth from Temperature in Region 3") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Year") +
ylab("MLD_T [m]") +
labs(caption = "Error bars are the mean +- SD. Red = Mean, Green = Median")
ggplot(NULL, aes(x = date_vec, y = MLD_T_Ptemp_Avg_Vec)) +
geom_path(data = NULL,aes(x = date_vec, y = MLD_T_Ptemp_Avg_Vec), color = "red") +
geom_point(data = NULL,aes(x = date_vec, y = MLD_T_Ptemp_Avg_Vec), color = "red", size = 4) +
geom_path(data = NULL,aes(x = date_vec, y = MLD_T_Ptemp_Med_Vec), color = "green") +
geom_point(data = NULL,aes(x = date_vec, y = MLD_T_Ptemp_Med_Vec), color = "green", size = 4) +
geom_point(data = BTMALL,aes(x=BTMALL$DATE,y=BTMALL$MLD_T_Ptemp_Avg),alpha = 0.1) +
geom_errorbar(aes(ymin = MLD_T_Ptemp_Avg_Vec - MLD_T_Ptemp_Sd_Vec, ymax = MLD_T_Ptemp_Avg_Vec + MLD_T_Ptemp_Sd_Vec),width = 50) +
ggtitle("Potential Temperature above MLD in Region 3") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Year") +
ylab("Potential Temperature above the MLD [C]") +
labs(caption = "Error bars are the mean +- SD. Red = Mean, Green = Median")
ggplot(NULL, aes(x = date_vec, y = MLD_T_SA_Avg_Vec)) +
geom_path(data = NULL,aes(x = date_vec, y = MLD_T_SA_Avg_Vec), color = "red") +
geom_point(data = NULL,aes(x = date_vec, y = MLD_T_SA_Avg_Vec), color = "red", size = 4) +
geom_path(data = NULL,aes(x = date_vec, y = MLD_T_SA_Med_Vec), color = "green") +
geom_point(data = NULL,aes(x = date_vec, y = MLD_T_SA_Med_Vec), color = "green", size = 4) +
geom_point(data = BTMALL,aes(x=BTMALL$DATE,y=BTMALL$MLD_T_SA_Avg),alpha = 0.1) +
geom_errorbar(aes(ymin = MLD_T_SA_Avg_Vec - MLD_T_SA_Sd_Vec, ymax = MLD_T_SA_Avg_Vec + MLD_T_SA_Sd_Vec),width = 50) +
ggtitle("Absolute Salinity above the MLD in Region 3") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Year") +
ylab("Absolute Salinity above the MLD [g/kg]") +
labs(caption = "Error bars are the mean +- SD. Red = Mean, Green = Median")
ggplot(NULL, aes(x = date_vec, y = MLD_T_Pden_Avg_Vec)) +
geom_path(data = NULL,aes(x = date_vec, y = MLD_T_Pden_Avg_Vec), color = "red") +
geom_point(data = NULL,aes(x = date_vec, y = MLD_T_Pden_Avg_Vec), color = "red", size = 4) +
geom_path(data = NULL,aes(x = date_vec, y = MLD_T_Pden_Med_Vec), color = "green") +
geom_point(data = NULL,aes(x = date_vec, y = MLD_T_Pden_Med_Vec), color = "green", size = 4) +
geom_point(data = BTMALL,aes(x=BTMALL$DATE,y=BTMALL$MLD_T_Pden_Avg),alpha = 0.1) +
geom_errorbar(aes(ymin = MLD_T_Pden_Avg_Vec - MLD_T_Pden_Sd_Vec, ymax = MLD_T_Pden_Avg_Vec + MLD_T_Pden_Sd_Vec),width = 50) +
ggtitle("Potential Density above the MLD in Region 3") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Year") +
ylab("Potential Density above the MLD [kg/m^3]") +
labs(caption = "Error bars are the mean +- SD. Red = Mean, Green = Median")
NO_MLD_D = BTMALL[which(is.na(MLD_D[BTMALL_idx])),]
basemap(limits = c(-71,-51,71,77),rotate = TRUE, bathymetry = FALSE, glaciers = FALSE) +
ggtitle("CTDs with no MLD_D Value", "Test inconclusive.") +
theme(plot.title = element_text(hjust = 0.5)) +
geom_spatial_point(data = NO_MLD_D, aes(x = LON, y = LAT, color = D))
## Using lon and lat as longitude and latitude columns, respectively.
## projection transformed from +init=epsg:4326 to +proj=stere +lat_0=90 +lat_ts=71 +lon_0=-61 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
## Assuming `crs = 4326` in stat_spatial_identity()
ggplot(NULL) +
geom_histogram(bins = 40, color="darkblue", fill="lightblue") +
aes(NO_MLD_D$D) +
xlab("Depth of CTDs with no MLD_D Value") +
ylab("Number of Casts") +
ggtitle("Max Depth Distribution of Casts with no MLD_D Value","Test inconclusive") +
theme(plot.title = element_text(hjust = 0.5))
# Fresh Water Analysis.
# Calculates the mean and the median of freshwater content by year, and then plots it on a graph containing every data point.
# TO DO: Make a proper legend.
fresh2021avg = mean(BTM2021$FRESH)
fresh2020avg = mean(BTM2020$FRESH)
fresh2019avg = mean(BTM2019$FRESH)
fresh2018avg = mean(BTM2018$FRESH)
fresh2017avg = mean(BTM2017$FRESH)
fresh2021sd = sd(BTM2021$FRESH)
fresh2020sd = sd(BTM2020$FRESH)
fresh2019sd = sd(BTM2019$FRESH)
fresh2018sd = sd(BTM2018$FRESH)
fresh2017sd = sd(BTM2017$FRESH)
fresh2021med = median(BTM2021$FRESH)
fresh2020med = median(BTM2020$FRESH)
fresh2019med = median(BTM2019$FRESH)
fresh2018med = median(BTM2018$FRESH)
fresh2017med = median(BTM2017$FRESH)
fresh_vec_avg = c(fresh2017avg,fresh2018avg,fresh2019avg,fresh2020avg,fresh2021avg)
fresh_vec_med = c(fresh2017med,fresh2018med,fresh2019med,fresh2020med,fresh2021med)
fresh_vec_sd = c(fresh2017sd,fresh2018sd,fresh2019sd,fresh2020sd,fresh2021sd)
freshData = data.frame(date_vec,fresh_vec_avg,fresh_vec_med,fresh_vec_sd)
colnames(freshData) = c("Year","FreshAvg","FreshMed","SD")
ggplot(NULL, aes(x = date_vec, y = fresh_vec_avg)) +
geom_path(data = freshData,aes(x = Year, y = FreshAvg), color = "red") +
geom_point(data = freshData,aes(x = Year, y = FreshAvg), color = "red", size = 4) +
geom_path(data = freshData,aes(x = Year, y = FreshMed), color = "green") +
geom_point(data = freshData,aes(x = Year, y = FreshMed), color = "green", size = 4) +
geom_point(data = BTMALL,aes(x=BTMALL$DATE,y=BTMALL$FRESH),alpha = 0.1) +
geom_errorbar(aes(ymin = fresh_vec_avg - fresh_vec_sd, ymax = fresh_vec_avg + fresh_vec_sd),width = 50) +
ggtitle("Annual Average Freshwater Content in Region 3","Background contains all datapoints in chronological order") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Year") +
ylab("Freshwater Content [m]") +
labs(caption = "Freshwater Content from Spall(2013). Error bars are the mean +- SD. Red = Mean, Green = Median")
## Warning: Use of `BTMALL$DATE` is discouraged. Use `DATE` instead.
## Warning: Use of `BTMALL$FRESH` is discouraged. Use `FRESH` instead.
# TMAX and TMIN Analysis.
# The profiles in my region included a warm fresh layer over a cold fresh layer over a warm salty core. Because of this, the tmax calculations are actually immediately below the MLD and not indicative of deep temperature maximum. Because of this, I calculated a second deep tmax characterized by being below the tmin found in most of the profiles of my region.
tmax2021avg = mean(BTM2021$TMAX)
tmax2020avg = mean(BTM2020$TMAX)
tmax2019avg = mean(BTM2019$TMAX)
tmax2018avg = mean(BTM2018$TMAX)
tmax2017avg = mean(BTM2017$TMAX)
tmax2021sd = sd(BTM2021$TMAX)
tmax2020sd = sd(BTM2020$TMAX)
tmax2019sd = sd(BTM2019$TMAX)
tmax2018sd = sd(BTM2018$TMAX)
tmax2017sd = sd(BTM2017$TMAX)
tmax2021med = median(BTM2021$TMAX)
tmax2020med = median(BTM2020$TMAX)
tmax2019med = median(BTM2019$TMAX)
tmax2018med = median(BTM2018$TMAX)
tmax2017med = median(BTM2017$TMAX)
tmax_vec_avg = c(tmax2017avg,tmax2018avg,tmax2019avg,tmax2020avg,tmax2021avg)
tmax_vec_med = c(tmax2017med,tmax2018med,tmax2019med,tmax2020med,tmax2021med)
tmax_vec_sd = c(tmax2017sd,tmax2018sd,tmax2019sd,tmax2020sd,tmax2021sd)
tmaxData = data.frame(date_vec,tmax_vec_avg,tmax_vec_med,tmax_vec_sd)
colnames(tmaxData) = c("Year","Avg","Med","SD")
ggplot(NULL, aes(x = date_vec, y = tmax_vec_avg)) +
geom_path(data = tmaxData,aes(x = Year, y = Avg), color = "red") +
geom_point(data = tmaxData,aes(x = Year, y = Avg), color = "red", size = 4) +
geom_path(data = tmaxData,aes(x = Year, y = Med), color = "green") +
geom_point(data = tmaxData,aes(x = Year, y = Med), color = "green", size = 4) +
geom_point(data = BTMALL,aes(x=BTMALL$DATE,y=BTMALL$TMAX),alpha = 0.1) +
geom_errorbar(aes(ymin = tmax_vec_avg - tmax_vec_sd, ymax = tmax_vec_avg + tmax_vec_sd),width = 50) +
ggtitle("Annual Average Temperature Maximum in Region 3","Background contains all datapoints in chronological order") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Year") +
ylab("Temperature Maximum [C]") +
labs(caption = "Error bars are the mean +- SD. Red = Mean, Green = Median. CAUTION: TMAX VALUES ARE UNRELIABLE. REFER TO DTMAX VALUES.")
## Warning: Use of `BTMALL$DATE` is discouraged. Use `DATE` instead.
## Warning: Use of `BTMALL$TMAX` is discouraged. Use `TMAX` instead.
tmaxD2021avg = mean(BTM2021$TMAXD)
tmaxD2020avg = mean(BTM2020$TMAXD)
tmaxD2019avg = mean(BTM2019$TMAXD)
tmaxD2018avg = mean(BTM2018$TMAXD)
tmaxD2017avg = mean(BTM2017$TMAXD)
tmaxD2021sd = sd(BTM2021$TMAXD)
tmaxD2020sd = sd(BTM2020$TMAXD)
tmaxD2019sd = sd(BTM2019$TMAXD)
tmaxD2018sd = sd(BTM2018$TMAXD)
tmaxD2017sd = sd(BTM2017$TMAXD)
tmaxD2021med = median(BTM2021$TMAXD)
tmaxD2020med = median(BTM2020$TMAXD)
tmaxD2019med = median(BTM2019$TMAXD)
tmaxD2018med = median(BTM2018$TMAXD)
tmaxD2017med = median(BTM2017$TMAXD)
tmaxD_vec_avg = c(tmaxD2017avg,tmaxD2018avg,tmaxD2019avg,tmaxD2020avg,tmaxD2021avg)
tmaxD_vec_med = c(tmaxD2017med,tmaxD2018med,tmaxD2019med,tmaxD2020med,tmaxD2021med)
tmaxD_vec_sd = c(tmaxD2017sd,tmaxD2018sd,tmaxD2019sd,tmaxD2020sd,tmaxD2021sd)
tmaxDData = data.frame(date_vec,tmaxD_vec_avg,tmaxD_vec_med,tmaxD_vec_sd)
colnames(tmaxDData) = c("Year","Avg","Med","SD")
ggplot(NULL, aes(x = date_vec, y = tmaxD_vec_avg)) +
geom_path(data = tmaxDData,aes(x = Year, y = -Avg), color = "red") +
geom_point(data = tmaxDData,aes(x = Year, y = -Avg), color = "red", size = 4) +
geom_path(data = tmaxDData,aes(x = Year, y = -Med), color = "green") +
geom_point(data = tmaxDData,aes(x = Year, y = -Med), color = "green", size = 4) +
geom_point(data = BTMALL,aes(x=BTMALL$DATE,y=-BTMALL$TMAXD),alpha = 0.1) +
geom_errorbar(aes(ymin = -tmaxD_vec_avg - tmaxD_vec_sd, ymax = -tmaxD_vec_avg + tmaxD_vec_sd),width = 50) +
ggtitle("Annual Average Temperature Maximum Depth in Region 3","Background contains all datapoints in chronological order") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Year") +
ylab("Depth [m]") +
labs(caption = "Error bars are the mean +- SD. Red = Mean, Green = Median. CAUTION: TMAX VALUES ARE UNRELIABLE. REFER TO DTMAX VALUES.")
## Warning: Use of `BTMALL$DATE` is discouraged. Use `DATE` instead.
## Warning: Use of `BTMALL$TMAXD` is discouraged. Use `TMAXD` instead.
tmaxSA2021avg = mean(BTM2021$TMAXSA)
tmaxSA2020avg = mean(BTM2020$TMAXSA)
tmaxSA2019avg = mean(BTM2019$TMAXSA)
tmaxSA2018avg = mean(BTM2018$TMAXSA)
tmaxSA2017avg = mean(BTM2017$TMAXSA,na.rm = TRUE)
tmaxSA2021sd = sd(BTM2021$TMAXSA)
tmaxSA2020sd = sd(BTM2020$TMAXSA)
tmaxSA2019sd = sd(BTM2019$TMAXSA)
tmaxSA2018sd = sd(BTM2018$TMAXSA)
tmaxSA2017sd = sd(BTM2017$TMAXSA,na.rm = TRUE)
tmaxSA2021med = median(BTM2021$TMAXSA)
tmaxSA2020med = median(BTM2020$TMAXSA)
tmaxSA2019med = median(BTM2019$TMAXSA)
tmaxSA2018med = median(BTM2018$TMAXSA)
tmaxSA2017med = median(BTM2017$TMAXSA,na.rm = TRUE)
tmaxSA_vec_avg = c(tmaxSA2017avg,tmaxSA2018avg,tmaxSA2019avg,tmaxSA2020avg,tmaxSA2021avg)
tmaxSA_vec_med = c(tmaxSA2017med,tmaxSA2018med,tmaxSA2019med,tmaxSA2020med,tmaxSA2021med)
tmaxSA_vec_sd = c(tmaxSA2017sd,tmaxSA2018sd,tmaxSA2019sd,tmaxSA2020sd,tmaxSA2021sd)
tmaxSAData = data.frame(date_vec,tmaxSA_vec_avg,tmaxSA_vec_med,tmaxSA_vec_sd)
colnames(tmaxSAData) = c("Year","Avg","Med","SD")
ggplot(NULL, aes(x = date_vec, y = tmaxSA_vec_avg)) +
geom_path(data = tmaxSAData,aes(x = Year, y = Avg), color = "red") +
geom_point(data = tmaxSAData,aes(x = Year, y = Avg), color = "red", size = 4) +
geom_path(data = tmaxSAData,aes(x = Year, y = Med), color = "green") +
geom_point(data = tmaxSAData,aes(x = Year, y = Med), color = "green", size = 4) +
geom_point(data = BTMALL,aes(x=BTMALL$DATE,y=BTMALL$TMAXSA),alpha = 0.1) +
geom_errorbar(aes(ymin = tmaxSA_vec_avg - tmaxSA_vec_sd, ymax = tmaxSA_vec_avg + tmaxSA_vec_sd),width = 50) +
ggtitle("Annual Average Temperature Maximum Salinity Anomoly in Region 3","Background contains all datapoints in chronological order") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Year") +
ylab("Absolute Salinity [g/kg]") +
labs(caption = "Error bars are the mean +- SD. Red = Mean, Green = Median. CAUTION: TMAX VALUES ARE UNRELIABLE. REFER TO DTMAX VALUES.")
## Warning: Use of `BTMALL$DATE` is discouraged. Use `DATE` instead.
## Warning: Use of `BTMALL$TMAXSA` is discouraged. Use `TMAXSA` instead.
## Warning: Removed 2 rows containing missing values (geom_point).
# -------------------------------------------------------------------------------------------------------
tmin2021avg = mean(BTM2021$TMIN)
tmin2020avg = mean(BTM2020$TMIN)
tmin2019avg = mean(BTM2019$TMIN)
tmin2018avg = mean(BTM2018$TMIN)
tmin2017avg = mean(BTM2017$TMIN)
tmin2021sd = sd(BTM2021$TMIN)
tmin2020sd = sd(BTM2020$TMIN)
tmin2019sd = sd(BTM2019$TMIN)
tmin2018sd = sd(BTM2018$TMIN)
tmin2017sd = sd(BTM2017$TMIN)
tmin2021med = median(BTM2021$TMIN)
tmin2020med = median(BTM2020$TMIN)
tmin2019med = median(BTM2019$TMIN)
tmin2018med = median(BTM2018$TMIN)
tmin2017med = median(BTM2017$TMIN)
tmin_vec_avg = c(tmin2017avg,tmin2018avg,tmin2019avg,tmin2020avg,tmin2021avg)
tmin_vec_med = c(tmin2017med,tmin2018med,tmin2019med,tmin2020med,tmin2021med)
tmin_vec_sd = c(tmin2017sd,tmin2018sd,tmin2019sd,tmin2020sd,tmin2021sd)
tminData = data.frame(date_vec,tmin_vec_avg,tmin_vec_med,tmin_vec_sd)
colnames(tminData) = c("Year","Avg","Med","SD")
ggplot(NULL, aes(x = date_vec, y = tmin_vec_avg)) +
geom_path(data = tminData,aes(x = Year, y = Avg), color = "red") +
geom_point(data = tminData,aes(x = Year, y = Avg), color = "red", size = 4) +
geom_path(data = tminData,aes(x = Year, y = Med), color = "green") +
geom_point(data = tminData,aes(x = Year, y = Med), color = "green", size = 4) +
geom_point(data = BTMALL,aes(x=BTMALL$DATE,y=BTMALL$TMIN),alpha = 0.1) +
geom_errorbar(aes(ymin = tmin_vec_avg - tmin_vec_sd, ymax = tmin_vec_avg + tmin_vec_sd),width = 50) +
ggtitle("Annual Average Temperature Minimum in Region 3","Background contains all datapoints in chronological order") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Year") +
ylab("Temperature Minimum [C]") +
labs(caption = "Error bars are the mean +- SD. Red = Mean, Green = Median.")
## Warning: Use of `BTMALL$DATE` is discouraged. Use `DATE` instead.
## Warning: Use of `BTMALL$TMIN` is discouraged. Use `TMIN` instead.
tminD2021avg = mean(BTM2021$TMIND)
tminD2020avg = mean(BTM2020$TMIND)
tminD2019avg = mean(BTM2019$TMIND)
tminD2018avg = mean(BTM2018$TMIND)
tminD2017avg = mean(BTM2017$TMIND)
tminD2021sd = sd(BTM2021$TMIND)
tminD2020sd = sd(BTM2020$TMIND)
tminD2019sd = sd(BTM2019$TMIND)
tminD2018sd = sd(BTM2018$TMIND)
tminD2017sd = sd(BTM2017$TMIND)
tminD2021med = median(BTM2021$TMIND)
tminD2020med = median(BTM2020$TMIND)
tminD2019med = median(BTM2019$TMIND)
tminD2018med = median(BTM2018$TMIND)
tminD2017med = median(BTM2017$TMIND)
tminD_vec_avg = c(tminD2017avg,tminD2018avg,tminD2019avg,tminD2020avg,tminD2021avg)
tminD_vec_med = c(tminD2017med,tminD2018med,tminD2019med,tminD2020med,tminD2021med)
tminD_vec_sd = c(tminD2017sd,tminD2018sd,tminD2019sd,tminD2020sd,tminD2021sd)
tminDData = data.frame(date_vec,tminD_vec_avg,tminD_vec_med,tminD_vec_sd)
colnames(tminDData) = c("Year","Avg","Med","SD")
ggplot(NULL, aes(x = date_vec, y = tminD_vec_avg)) +
geom_path(data = tminDData,aes(x = Year, y = -Avg), color = "red") +
geom_point(data = tminDData,aes(x = Year, y = -Avg), color = "red", size = 4) +
geom_path(data = tminDData,aes(x = Year, y = -Med), color = "green") +
geom_point(data = tminDData,aes(x = Year, y = -Med), color = "green", size = 4) +
geom_point(data = BTMALL,aes(x=BTMALL$DATE,y=-BTMALL$TMIND),alpha = 0.1) +
geom_errorbar(aes(ymin = -tminD_vec_avg - tminD_vec_sd, ymax = -tminD_vec_avg + tminD_vec_sd),width = 50) +
ggtitle("Annual Average Temperature Minimum Depth in Region 3","Background contains all datapoints in chronological order") +
theme(plot.title = element_text(hjust = 0.5)) +
ylim(-200,0)+
xlab("Year") +
ylab("Depth [m]") +
labs(caption = "Error bars are the mean +- SD. Red = Mean, Green = Median.")
## Warning: Use of `BTMALL$DATE` is discouraged. Use `DATE` instead.
## Warning: Use of `BTMALL$TMIND` is discouraged. Use `TMIND` instead.
tminSA2021avg = mean(BTM2021$TMINSA)
tminSA2020avg = mean(BTM2020$TMINSA)
tminSA2019avg = mean(BTM2019$TMINSA)
tminSA2018avg = mean(BTM2018$TMINSA,na.rm=TRUE)
tminSA2017avg = mean(BTM2017$TMINSA)
tminSA2021sd = sd(BTM2021$TMINSA)
tminSA2020sd = sd(BTM2020$TMINSA)
tminSA2019sd = sd(BTM2019$TMINSA)
tminSA2018sd = sd(BTM2018$TMINSA,na.rm=TRUE)
tminSA2017sd = sd(BTM2017$TMINSA)
tminSA2021med = median(BTM2021$TMINSA)
tminSA2020med = median(BTM2020$TMINSA)
tminSA2019med = median(BTM2019$TMINSA)
tminSA2018med = median(BTM2018$TMINSA,na.rm=TRUE)
tminSA2017med = median(BTM2017$TMINSA)
tminSA_vec_avg = c(tminSA2017avg,tminSA2018avg,tminSA2019avg,tminSA2020avg,tminSA2021avg)
tminSA_vec_med = c(tminSA2017med,tminSA2018med,tminSA2019med,tminSA2020med,tminSA2021med)
tminSA_vec_sd = c(tminSA2017sd,tminSA2018sd,tminSA2019sd,tminSA2020sd,tminSA2021sd)
tminSAData = data.frame(date_vec,tminSA_vec_avg,tminSA_vec_med,tminSA_vec_sd)
colnames(tminSAData) = c("Year","Avg","Med","SD")
ggplot(NULL, aes(x = date_vec, y = tminSA_vec_avg)) +
geom_path(data = tminSAData,aes(x = Year, y = Avg), color = "red") +
geom_point(data = tminSAData,aes(x = Year, y = Avg), color = "red", size = 4) +
geom_path(data = tminSAData,aes(x = Year, y = Med), color = "green") +
geom_point(data = tminSAData,aes(x = Year, y = Med), color = "green", size = 4) +
geom_point(data = BTMALL,aes(x=BTMALL$DATE,y=BTMALL$TMINSA),alpha = 0.1) +
geom_errorbar(aes(ymin = tminSA_vec_avg - tminSA_vec_sd, ymax = tminSA_vec_avg + tminSA_vec_sd),width = 50) +
ggtitle("Annual Average Temperature Minimum Salinity Anomoly in Region 3","Background contains all datapoints in chronological order") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Year") +
ylab("Absolute Salinity [g/kg]") +
labs(caption = "Error bars are the mean +- SD. Red = Mean, Green = Median.")
## Warning: Use of `BTMALL$DATE` is discouraged. Use `DATE` instead.
## Warning: Use of `BTMALL$TMINSA` is discouraged. Use `TMINSA` instead.
## Warning: Removed 2 rows containing missing values (geom_point).
# -------------------------------------------------------------------------------------------------------
dtmax2021avg = mean(BTM2021$DTMAX)
dtmax2020avg = mean(BTM2020$DTMAX)
dtmax2019avg = mean(BTM2019$DTMAX)
dtmax2018avg = mean(BTM2018$DTMAX)
dtmax2017avg = mean(BTM2017$DTMAX)
dtmax2021sd = sd(BTM2021$DTMAX)
dtmax2020sd = sd(BTM2020$DTMAX)
dtmax2019sd = sd(BTM2019$DTMAX)
dtmax2018sd = sd(BTM2018$DTMAX)
dtmax2017sd = sd(BTM2017$DTMAX)
dtmax2021med = median(BTM2021$DTMAX)
dtmax2020med = median(BTM2020$DTMAX)
dtmax2019med = median(BTM2019$DTMAX)
dtmax2018med = median(BTM2018$DTMAX)
dtmax2017med = median(BTM2017$DTMAX)
dtmax_vec_avg = c(dtmax2017avg,dtmax2018avg,dtmax2019avg,dtmax2020avg,dtmax2021avg)
dtmax_vec_med = c(dtmax2017med,dtmax2018med,dtmax2019med,dtmax2020med,dtmax2021med)
dtmax_vec_sd = c(dtmax2017sd,dtmax2018sd,dtmax2019sd,dtmax2020sd,dtmax2021sd)
dtmaxData = data.frame(date_vec,dtmax_vec_avg,dtmax_vec_med,dtmax_vec_sd)
colnames(dtmaxData) = c("Year","Avg","Med","SD")
ggplot(NULL, aes(x = date_vec, y = dtmax_vec_avg)) +
geom_path(data = dtmaxData,aes(x = Year, y = Avg), color = "red") +
geom_point(data = dtmaxData,aes(x = Year, y = Avg), color = "red", size = 4) +
geom_path(data = dtmaxData,aes(x = Year, y = Med), color = "green") +
geom_point(data = dtmaxData,aes(x = Year, y = Med), color = "green", size = 4) +
geom_point(data = BTMALL,aes(x=BTMALL$DATE,y=BTMALL$DTMAX),alpha = 0.1) +
geom_errorbar(aes(ymin = dtmax_vec_avg - dtmax_vec_sd, ymax = dtmax_vec_avg + dtmax_vec_sd),width = 50) +
ggtitle("Annual Average Deep Temperature Maximum in Region 3","Background contains all datapoints in chronological order") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Year") +
ylab("Temperature Maximum [C]") +
labs(caption = "Error bars are the mean +- SD. Red = Mean, Green = Median.")
## Warning: Use of `BTMALL$DATE` is discouraged. Use `DATE` instead.
## Warning: Use of `BTMALL$DTMAX` is discouraged. Use `DTMAX` instead.
dtmaxd2021avg = mean(BTM2021$DTMAXD)
dtmaxd2020avg = mean(BTM2020$DTMAXD)
dtmaxd2019avg = mean(BTM2019$DTMAXD)
dtmaxd2018avg = mean(BTM2018$DTMAXD)
dtmaxd2017avg = mean(BTM2017$DTMAXD)
dtmaxd2021sd = sd(BTM2021$DTMAXD)
dtmaxd2020sd = sd(BTM2020$DTMAXD)
dtmaxd2019sd = sd(BTM2019$DTMAXD)
dtmaxd2018sd = sd(BTM2018$DTMAXD)
dtmaxd2017sd = sd(BTM2017$DTMAXD)
dtmaxd2021med = median(BTM2021$DTMAXD)
dtmaxd2020med = median(BTM2020$DTMAXD)
dtmaxd2019med = median(BTM2019$DTMAXD)
dtmaxd2018med = median(BTM2018$DTMAXD)
dtmaxd2017med = median(BTM2017$DTMAXD)
dtmaxd_vec_avg = c(dtmaxd2017avg,dtmaxd2018avg,dtmaxd2019avg,dtmaxd2020avg,dtmaxd2021avg)
dtmaxd_vec_med = c(dtmaxd2017med,dtmaxd2018med,dtmaxd2019med,dtmaxd2020med,dtmaxd2021med)
dtmaxd_vec_sd = c(dtmaxd2017sd,dtmaxd2018sd,dtmaxd2019sd,dtmaxd2020sd,dtmaxd2021sd)
dtmaxdData = data.frame(date_vec,dtmaxd_vec_avg,dtmaxd_vec_med,dtmaxd_vec_sd)
colnames(dtmaxdData) = c("Year","Avg","Med","SD")
ggplot(NULL, aes(x = date_vec, y = dtmaxd_vec_avg)) +
geom_path(data = dtmaxdData,aes(x = Year, y = Avg), color = "red") +
geom_point(data = dtmaxdData,aes(x = Year, y = Avg), color = "red", size = 4) +
geom_path(data = dtmaxdData,aes(x = Year, y = Med), color = "green") +
geom_point(data = dtmaxdData,aes(x = Year, y = Med), color = "green", size = 4) +
geom_point(data = BTMALL,aes(x=BTMALL$DATE,y=BTMALL$DTMAXD),alpha = 0.1) +
geom_errorbar(aes(ymin = dtmaxd_vec_avg - dtmaxd_vec_sd, ymax = dtmaxd_vec_avg + dtmaxd_vec_sd),width = 50) +
ggtitle("Annual Average Deep Temperature Maximum in Region 3","Background contains all datapoints in chronological order") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Year") +
ylab("Temperature Maximum Depth [m]") +
labs(caption = "Error bars are the mean +- SD. Red = Mean, Green = Median.")
## Warning: Use of `BTMALL$DATE` is discouraged. Use `DATE` instead.
## Warning: Use of `BTMALL$DTMAXD` is discouraged. Use `DTMAXD` instead.
dtmaxSA2021avg = mean(BTM2021$DTMAXSA)
dtmaxSA2020avg = mean(BTM2020$DTMAXSA,na.rm=TRUE)
dtmaxSA2019avg = mean(BTM2019$DTMAXSA)
dtmaxSA2018avg = mean(BTM2018$DTMAXSA)
dtmaxSA2017avg = mean(BTM2017$DTMAXSA,na.rm=TRUE)
dtmaxSA2021sd = sd(BTM2021$DTMAXSA)
dtmaxSA2020sd = sd(BTM2020$DTMAXSA,na.rm=TRUE)
dtmaxSA2019sd = sd(BTM2019$DTMAXSA)
dtmaxSA2018sd = sd(BTM2018$DTMAXSA)
dtmaxSA2017sd = sd(BTM2017$DTMAXSA,na.rm=TRUE)
dtmaxSA2021med = median(BTM2021$DTMAXSA)
dtmaxSA2020med = median(BTM2020$DTMAXSA,na.rm=TRUE)
dtmaxSA2019med = median(BTM2019$DTMAXSA)
dtmaxSA2018med = median(BTM2018$DTMAXSA)
dtmaxSA2017med = median(BTM2017$DTMAXSA,na.rm=TRUE)
dtmaxSA_vec_avg = c(dtmaxSA2017avg,dtmaxSA2018avg,dtmaxSA2019avg,dtmaxSA2020avg,dtmaxSA2021avg)
dtmaxSA_vec_med = c(dtmaxSA2017med,dtmaxSA2018med,dtmaxSA2019med,dtmaxSA2020med,dtmaxSA2021med)
dtmaxSA_vec_sd = c(dtmaxSA2017sd,dtmaxSA2018sd,dtmaxSA2019sd,dtmaxSA2020sd,dtmaxSA2021sd)
dtmaxSAData = data.frame(date_vec,dtmaxSA_vec_avg,dtmaxSA_vec_med,dtmaxSA_vec_sd)
colnames(dtmaxSAData) = c("Year","Avg","Med","SD")
ggplot(NULL, aes(x = date_vec, y = tmaxSA_vec_avg)) +
geom_path(data = dtmaxSAData,aes(x = Year, y = Avg), color = "red") +
geom_point(data = dtmaxSAData,aes(x = Year, y = Avg), color = "red", size = 4) +
geom_path(data = dtmaxSAData,aes(x = Year, y = Med), color = "green") +
geom_point(data = dtmaxSAData,aes(x = Year, y = Med), color = "green", size = 4) +
geom_point(data = BTMALL,aes(x=BTMALL$DATE,y=BTMALL$DTMAXSA),alpha = 0.1) +
geom_errorbar(aes(ymin = dtmaxSA_vec_avg - dtmaxSA_vec_sd, ymax = dtmaxSA_vec_avg + dtmaxSA_vec_sd),width = 50) +
ggtitle("Annual Average Deep Temperature Maximum Salinity Anomoly in Region 3","Background contains all datapoints in chronological order") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Year") +
ylab("Absolute Salinity [g/kg]") +
labs(caption = "Error bars are the mean +- SD. Red = Mean, Green = Median.")
## Warning: Use of `BTMALL$DATE` is discouraged. Use `DATE` instead.
## Warning: Use of `BTMALL$DTMAXSA` is discouraged. Use `DTMAXSA` instead.
## Warning: Removed 2 rows containing missing values (geom_point).
export = data.frame(BTMALL$LAT,BTMALL$LON,BTMALL$DATE,BTMALL$TMAX,BTMALL$TMAXD,BTMALL$TMAXSA)
colnames(export) = c("LAT","LON","DATE","TMAX","TMAXD","TMAXSA")
write.table(export, file = "mcappola_tmaxdata.dat", sep = "\t")
export2 = data.frame(BTMALL$LAT,BTMALL$LON,BTMALL$DATE,BTMALL$DTMAX,BTMALL$DTMAXD,BTMALL$DTMAXSA)
colnames(export) = c("LAT","LON","DATE","DTMAX","DTMAXD","DTMAXSA")
write.table(export, file = "mcappola_dtmaxdata.dat", sep = "\t")
# Multi-Plot +- 300m
# Updated (10/28) from 2021 data, to all data (2017-2021).
# TO DO: Establish a way to average all profiles into a single "Generic Region 3" profile.
a1 = BTMALL
b1 = which(a1$D <= 300)
c1 = a1$CAST[b1]
SHALLOW_idx = which(OMGALL$CAST == c1)
## Warning in OMGALL$CAST == c1: longer object length is not a multiple of shorter
## object length
SHALLOW = OMGALL[SHALLOW_idx,]
ggplot(SHALLOW,aes(x=T,y=-D)) +
xlim(-2,6) +
geom_path(group = SHALLOW$CAST) +
ggtitle("2017 - 2021 Shallow Temperature Profiles (<300m)") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Temp, C") +
ylab("Depth, m")
## Warning: Removed 6 row(s) containing missing values (geom_path).
ggplot(SHALLOW,aes(x=S,y=-D)) +
xlim(26,35) +
geom_path(group = SHALLOW$CAST) +
ggtitle("2017 - 2021 Shallow Salinity Profiles (<300m)") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Salinity, PSU") +
ylab("Depth, m")
a1 = BTM2021
b2 = which(a1$D >= 300)
c2 = a1$CAST[b2]
DEEP_idx = which(OMGALL$CAST == c2)
## Warning in OMGALL$CAST == c2: longer object length is not a multiple of shorter
## object length
DEEP = OMGALL[DEEP_idx,]
ggplot(DEEP,aes(x=T,y=-D)) +
#xlim(-2,6) +
geom_path(group = DEEP$CAST) +
ggtitle("2017 - 2021 Deep Temperature Profiles (>300m)") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Temp, C") +
ylab("Depth, m")
## Warning: Removed 14 row(s) containing missing values (geom_path).
ggplot(DEEP,aes(x=S,y=-D)) +
xlim(26,35) +
geom_path(group = DEEP$CAST) +
ggtitle("2017 - 2021 Deep Salinity Profiles (>300m)") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Salinity, PSU") +
ylab("Depth, m")
## Warning: Removed 14 row(s) containing missing values (geom_path).
ggplot(SHALLOW,mapping=aes(x=S,y=T,color=D)) +
xlim(26,35) +
ylim(-2,6) +
geom_point() +
ggtitle("2017 - 2021 Shallow T(S) Profiles (<300m)") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Salinity, psu") +
ylab("Temperature, C") +
scale_color_viridis_c(guide=guide_colorbar(barheight=15))
## Warning: Removed 22 rows containing missing values (geom_point).
ggplot(DEEP,mapping=aes(x=S,y=T,color=D)) +
xlim(26,35) +
ylim(-2,6) +
geom_point() +
ggtitle("2017 - 2021 Deep T(S) Profiles (>300m)") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Salinity, psu") +
ylab("Temperature, C") +
scale_color_viridis_c(guide=guide_colorbar(barheight=15))
## Warning: Removed 119 rows containing missing values (geom_point).
rm(a1,b1,c1,b2,c2)