library(dplyr)
##
## 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
library(tidyr)
library(ggplot2)
library(leaflet)
library(ggspatial)
library(ggOceanMaps)
## Registered S3 methods overwritten by 'stars':
## method from
## st_bbox.SpatRaster sf
## st_crs.SpatRaster sf
library(ggOceanMapsData)
library(gsw)
library(matrixStats)
##
## 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.
# 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. Splits this frame into yearly data sets.
a = diff(OMGALL$D)
b = which(a < 0)
BTMALL_idx = append(b,length(a))
BTMALL = OMGALL[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",]
rm(a,b)
# Freshwater Content.
# For loop that records the freshwater content for each profile. Chunk than adds
# this vector to the BTMALL data set, then splits this frame into yearly data sets.
# TO DO: Binned values might be helpful for other functions. Perhaps make a new
# data.frame of binned values only, and record/append them in this loop. Will
# need to learn how to do this.
fresh = vector()
for (i in BTMALL_idx) {
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)
s_ref = 34.8
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)
rm(data,Depth,Temp,Salt,maxDepth,binDepth,binSalt,binTemp,s_ref_vec,s_ref)
}
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",]
fresh2021avg = mean(BTM2021$FRESH)
fresh2020avg = mean(BTM2020$FRESH)
fresh2019avg = mean(BTM2019$FRESH)
fresh2018avg = mean(BTM2018$FRESH)
fresh2017avg = mean(BTM2017$FRESH)
fresh_vec_avg = c(fresh2017avg,fresh2018avg,fresh2019avg,fresh2020avg,fresh2021avg)
freshAVG = data.frame(fresh_vec_avg,c(mean(BTM2017$DATE),mean(BTM2018$DATE),mean(BTM2019$DATE),mean(BTM2020$DATE),mean(BTM2021$DATE)))
colnames(freshAVG) = c("Fresh","Year")
ggplot(NULL) +
geom_path(data = freshAVG,aes(x = Year, y = Fresh), color = "red") +
geom_point(data = freshAVG,aes(x = Year, y = Fresh), color = "red", size = 4) +
geom_point(data = BTMALL,aes(x=BTMALL$DATE,y=BTMALL$FRESH),alpha = 0.2) +
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). Caution: Data has not been checked for even spatial distribution yet.")
## Warning: Use of `BTMALL$DATE` is discouraged. Use `DATE` instead.
## Warning: Use of `BTMALL$FRESH` is discouraged. Use `FRESH` instead.
# Mixed Layer Depth.
# For loop that records the MLD from a Temperature threshold method and a
# Potential Density threshold method. Then updates the results to the BTMALL
# data frame.
# TO DO: Visually check 15 profiles from different areas in the region to
# verify accuracy and compare methods. Search for profiles that have NA for both
# methods. If any, why? Search for profiles where the two methods are extremely
# different. If so, which is better? Make logic code that chooses one for graphing
# and make code that can display a graph with MLD reference lines. Make a graph
# showing if avg MLD has changed in the data set. Calculate the MLD properties,
# see/show how they have changed.
# BIGBRAIN STUFF: Combine this chunk with the freshwater content loop.
temp_crit = 0.2 # Thompson (1976)
pden_crit = 0.125 # Miller (1976)
MLD_T = vector()
MLD_T_Ptemp_Avg = vector()
MLD_T_Ptemp_SD = vector()
MLD_T_SA_Avg = vector()
MLD_T_SA_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()
# N2_min = vector()
for (i in BTMALL_idx) {
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))]
# N2_min[i] = data_N2_min
SST = binTemp[1]
a = binDepth[which(binTemp <= SST - temp_crit)]
MLD_T[i] = a[1]
if (is.na(MLD_T[i]) == FALSE) {
MLD_T_Ptemp_Avg[i] = mean(binTemp[1:a[1]])
MLD_T_Ptemp_SD[i] = sd((binTemp[1:a[1]]))
MLD_T_SA_Avg[i] = mean(binSA[1:a[1]])
MLD_T_SA_SD[i] = sd(binSA[1:a[1]])
}
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
}
a = binDepth[which(binPDen >= binPDen[1]*pden_crit + binPDen[1])]
MLD_D[i] = a[1]
if (is.na(MLD_D[i]) == FALSE) {
MLD_D_Ptemp_Avg[i] = mean(binTemp[1:a[1]])
MLD_D_Ptemp_SD[i] = sd((binTemp[1:a[1]]))
MLD_D_SA_Avg[i] = mean(binSA[1:a[1]])
MLD_D_SA_SD[i] = sd(binSA[1:a[1]])
}
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_T = MLD_T[BTMALL_idx]
MLD_T_Ptemp_Avg = MLD_T_Ptemp_Avg[BTMALL_idx]
MLD_T_Ptemp_SD = MLD_T_Ptemp_SD[BTMALL_idx]
MLD_T_SA_Avg = MLD_T_SA_Avg[BTMALL_idx]
MLD_T_SA_SD = MLD_T_SA_SD[BTMALL_idx]
MLD_D = MLD_D[BTMALL_idx]
MLD_D_Ptemp_Avg = MLD_D_Ptemp_Avg[BTMALL_idx]
MLD_D_Ptemp_SD = MLD_D_Ptemp_SD[BTMALL_idx]
MLD_D_SA_Avg = MLD_D_SA_Avg[BTMALL_idx]
MLD_D_SA_SD = MLD_D_SA_SD[BTMALL_idx]
#N2_min = N2_min[BTMALL_idx]
BTMALL$MLD_T = MLD_T
BTMALL$MLD_D = MLD_D
BTMALL$MLD_T_Ptemp_Avg = MLD_T_Ptemp_Avg
BTMALL$MLD_T_Ptemp_SD = MLD_T_Ptemp_SD
BTMALL$MLD_T_SA_Avg = MLD_T_SA_Avg
BTMALL$MLD_T_SA_SD = MLD_T_SA_SD
BTMALL$MLD_D_Ptemp_Avg = MLD_D_Ptemp_Avg
BTMALL$MLD_D_Ptemp_SD = MLD_D_Ptemp_SD
BTMALL$MLD_D_SA_Avg = MLD_D_SA_Avg
BTMALL$MLD_D_SA_SD = MLD_D_SA_SD
# 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
# Single Profile Plots.
# TEST is the cast number from OMGALL for the desired profiles.
TEST = 186.1
temp_crit = 0.2 # Thompson (1976)
pden_crit = 0.125 # Miller (1976)
data = filter(OMGALL,CAST == TEST)
mld = which(BTMALL$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))]
SST = binTemp[1]
a = binDepth[which(binTemp <= SST - temp_crit)]
MLD_T_TEST = a[1]
a = binDepth[which(binPDen >= binPDen[1]*pden_crit + binPDen[1])]
MLD_D_TEST = a[1]
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")
## Warning: Use of `data$D` is discouraged. Use `D` instead.
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 = data, aes(x=T,y=-D)) +
ggtitle("Temperature Profile with MLD_T",data$CAST) +
geom_hline(BTMALL,yintercept = -BTMALL[mld,]$MLD_T, color = "red") +
geom_vline(BTMALL,xintercept = BTMALL[mld,]$MLD_T_Ptemp_Avg, color = "blue") +
geom_hline(NULL,yintercept = -data_N2_min, color = "darkgreen") +
ylim(-(BTMALL[mld,]$MLD_T + 50), 0)+
xlim(-2, 6) +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Temperature [C]") +
ylab("Depth, [m]") +
annotate("text", x = -1.5, y = -(data_N2_min-5), label = "N^2 Min", color = "darkgreen")+
annotate("text", x = -1.5, y = -(BTMALL[mld,]$MLD_T -5), label = "MLD_T", color = "red")+
annotate("text", x = BTMALL[mld,]$MLD_T_Ptemp_Avg + 0.7, y = 0, label = "PTemp Avg", color = "blue")+
annotate("text", x = -0.7, y = -(data_N2_min-5), label = round(data_N2_min, digits = 2), color = "darkgreen")+
annotate("text", x = -0.8, y = -(BTMALL[mld,]$MLD_T -5), label = BTMALL[mld,]$MLD_T, color = "red")+
annotate("text", x = BTMALL[mld,]$MLD_T_Ptemp_Avg + 1.6, y = 0, label = round(BTMALL[mld,]$MLD_T_Ptemp_Avg, digits = 2), color = "blue")
## Warning: geom_hline(): Ignoring `mapping` because `yintercept` was provided.
## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.
## Warning: Removed 6140 row(s) containing missing values (geom_path).
ggplot(NULL) +
geom_path(data = data,aes(x=S,y=-D)) +
ggtitle("Salinity Profile with MLD_T",data$CAST) +
geom_hline(BTMALL,yintercept = -BTMALL[mld,]$MLD_T, color = "red") +
geom_vline(BTMALL,xintercept = BTMALL[mld,]$MLD_T_SA_Avg, color = "blue") +
geom_hline(NULL,yintercept = -data_N2_min, color = "green")+
ylim(-(BTMALL[mld,]$MLD_T + 50), 0)+
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Salinity [g/kg]") +
ylab("Depth [m]") +
annotate("text", x = data$S[1]+0.5, y = 0, label = "N^2 Min", color = "darkgreen")+
annotate("text", x = data$S[1]+0.5, y = -(BTMALL[mld,]$MLD_T -5), label = "MLD_T", color = "red") +
annotate("text", x = BTMALL[mld,]$MLD_T_SA_Avg + 0.7, y = 0, label = "SA Avg", color = "blue")+
annotate("text", x = data$S[1]+1, y = 0, label = round(data_N2_min, digits = 2), color = "darkgreen")+
annotate("text", x = data$S[1]+1, y = -(BTMALL[mld,]$MLD_T -5), label = BTMALL[mld,]$MLD_T, color = "red")+
annotate("text", x = BTMALL[mld,]$MLD_T_SA_Avg + 1.6, y = 0, label = round(BTMALL[mld,]$MLD_T_SA_Avg, digits = 2), color = "blue")
## Warning: geom_hline(): Ignoring `mapping` because `yintercept` was provided.
## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.
## Warning: Removed 6140 row(s) containing missing values (geom_path).
ggplot(NULL) +
geom_path(data = data, aes(x=T,y=-D)) +
ggtitle("Temperature Profile with MLD_D",data$CAST) +
geom_hline(BTMALL,yintercept = -BTMALL[mld,]$MLD_D, color = "red") +
geom_vline(BTMALL,xintercept = BTMALL[mld,]$MLD_D_Ptemp_Avg, color = "blue") +
geom_hline(NULL,yintercept = -data_N2_min, color = "green")+
ylim(-(BTMALL[mld,]$MLD_D + 50), 0)+
xlim(-2, 6) +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Temperature [C]") +
ylab("Depth, [m]") +
annotate("text", x = -1.5, y = 0, label = "N^2 Min", color = "darkgreen")+
annotate("text", x = -1.5, y = -(BTMALL[mld,]$MLD_D -5), label = "MLD_D", color = "red")+
annotate("text", x = BTMALL[mld,]$MLD_D_Ptemp_Avg + 0.7, y = 0, label = "PTemp Avg", color = "blue")+
annotate("text", x = -0.7, y = 0, label = round(data_N2_min, digits = 2), color = "darkgreen")+
annotate("text", x = -0.8, y = -(BTMALL[mld,]$MLD_D -5), label = BTMALL[mld,]$MLD_D, color = "red")+
annotate("text", x = BTMALL[mld,]$MLD_D_Ptemp_Avg + 1.6, y = 0, label = round(BTMALL[mld,]$MLD_D_Ptemp_Avg, digits = 2), color = "blue")
## Warning: geom_hline(): Ignoring `mapping` because `yintercept` was provided.
## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.
## Warning: Removed 4934 row(s) containing missing values (geom_path).
ggplot(NULL) +
geom_path(data = data,aes(x=S,y=-D)) +
ggtitle("Salinity Profile with MLD_D",data$CAST) +
geom_hline(BTMALL,yintercept = -BTMALL[mld,]$MLD_D, color = "red") +
geom_vline(BTMALL,xintercept = BTMALL[mld,]$MLD_D_SA_Avg, color = "blue") +
geom_hline(NULL,yintercept = -(data_N2_min), color = "green")+
ylim(-(BTMALL[mld,]$MLD_D + 50), 0)+
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Salinity [g/kg]") +
ylab("Depth [m]") +
annotate("text", x = data$S[1]+0.5, y = 0, label = "N^2 Min", color = "darkgreen")+
annotate("text", x = data$S[1]+0.5, y = -(BTMALL[mld,]$MLD_D -5), label = "MLD_D", color = "red") +
annotate("text", x = BTMALL[mld,]$MLD_D_SA_Avg + 0.7, y = 0, label = "SA Avg", color = "blue")+
annotate("text", x = data$S[1]+1, y = 0, label = round(data_N2_min, digits = 2), color = "darkgreen")+
annotate("text", x = data$S[1]+1, y = -(BTMALL[mld,]$MLD_D -5), label = BTMALL[mld,]$MLD_D, color = "red")+
annotate("text", x = BTMALL[mld,]$MLD_D_SA_Avg + 1.6, y = 0, label = round(BTMALL[mld,]$MLD_D_SA_Avg, digits = 2), color = "blue")
## Warning: geom_hline(): Ignoring `mapping` because `yintercept` was provided.
## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.
## Warning: Removed 4934 row(s) containing missing values (geom_path).
# 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)
#MAPS. I am using the ggOceanMaps package, which requires the ggOceanMapsData package to run.
# I recommend making Bathymetry and Glaciers = "FALSE" when experimenting as they take up
# a lot of computing resources.
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")
## Warning in sp::proj4string(shapefiles$land): CRS object has comment, which is
## lost in output
## 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()
OMGALL_SINGULAR = OMGALL[BTMALL_idx,]
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 = OMGALL_SINGULAR, aes(x = LON, y = LAT, color = D))
## Warning in sp::proj4string(shapefiles$land): CRS object has comment, which is
## lost in output
## 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()
# Binned Value Graphs.
# Retained code for future use. Will need to make a new data.frame containing only
# binned values before this can be re-adapted.
#ggplot(NULL) +
# geom_path(aes(binTemp,-binDepth)) +
# ggtitle("Binned Temperature Profile") +
# theme(plot.title = element_text(hjust = 0.5)) +
# xlab("binTemp, C") +
# ylab("Depth, m")
#ggplot(NULL) +
# geom_path(aes(binSalt,-binDepth)) +
# ggtitle("Binned Salinity Profile w/ Freshwater Content Reference Line") +
# theme(plot.title = element_text(hjust = 0.5)) +
# xlab("binSalt, g/kg") +
# ylab("Depth, m") +
# geom_hline(aes(yintercept = -fresh,color = fresh))+
# labs(colour = "Freshwater Content")