OMG00

#install.packages("dplyr")
#install.packages("tidyr")
#install.packages("leaflet")
#install.packages("knitr")
#install.packages("rmarkdown")
#install.packages("gsw")
#install.packages("matrixStats")
#install.packages("ggspatial")

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

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)

## Warning: package 'tidyr' was built under R version 4.1.2

library (ggplot2)
library (leaflet)

## Warning: package 'leaflet' was built under R version 4.1.2

library (gsw)

## Warning: package 'gsw' was built under R version 4.1.2

library (matrixStats)

## Warning: package 'matrixStats' was built under R version 4.1.2

## 
## Attaching package: 'matrixStats'

## The following object is masked from 'package:dplyr':
## 
##     count

##shell("cd C:/Users/ingri/MAST467/Data/2016/output16.dat")
shell("C:/Users/ingri/MAST467/Data/2016/N16.cmd", translate=TRUE)
Omg16<-read.csv("C:/Users/ingri/MAST467/Data/2016/output16.dat",sep=" ")

##shell("cd C:/Users/ingri/MAST467/Data/2017/output17.dat")
shell("C:/Users/ingri/MAST467/Data/2017/N17.cmd", translate=TRUE)
Omg17<-read.csv("C:/Users/ingri/MAST467/Data/2017/output17.dat",sep=" ")

##shell("cd C:/Users/ingri/MAST467/Data/2018/output18.dat")
shell("C:/Users/ingri/MAST467/Data/2018/N18.cmd", translate=TRUE)
Omg18<-read.csv("C:/Users/ingri/MAST467/Data/2018/output18.dat",sep=" ")

##shell("cd C:/Users/ingri/MAST467/Data/2019/output19.dat")
shell("C:/Users/ingri/MAST467/Data/2019/N19.cmd", translate=TRUE)
Omg19<-read.csv("C:/Users/ingri/MAST467/Data/2019/output19.dat",sep=" ")

##shell("cd C:/Users/ingri/MAST467/Data/2020/output20.dat")
shell("C:/Users/ingri/MAST467/Data/2020/N20.cmd", translate=TRUE)
Omg20<-read.csv("C:/Users/ingri/MAST467/Data/2020/output20.dat",sep=" ")

##shell("cd C:/Users/ingri/MAST467/Data/2021/output21.dat")
shell("C:/Users/ingri/MAST467/Data/2021/N21.cmd", translate=TRUE)
Omg21<-read.csv("C:/Users/ingri/MAST467/Data/2021/output21.dat",sep=" ")

p16= gsw_p_from_z(-Omg16$dep, Omg16$lat)
SA16 = gsw_SA_from_SP(Omg16$sal, p16, Omg16$lat, Omg16$lon)
pt16 = gsw_pt_from_t(SA16, Omg16$tem, p16, 0)
CT16= gsw_CT_from_t(SA16, Omg16$tem, p16)


dat16.1<-data.frame(Omg16$cast, Omg16$lon, Omg16$lat, Omg16$dep, Omg16$sal, Omg16$tem, p16, SA16, pt16, CT16)

colnames(dat16.1) = c("cast16", "lon16", "lat16", "dep16", "sal16", "tem16", "p16", "SA16", "pt16", "CT16")

###Many steps are this chunk is borrowed from Sam's code
###Defining empty vectors/variables/etc. BEFORE for loop begins...
Dmax16 = 0
fresh16 = 0
dat_all = NULL
y16= NULL #y is an empty dataframe that we will store the output of each iteration for

###Looping through each cast to define bin averaged variables
for (i in 1:17) #i refers to cast
{data16<-dplyr::filter(dat16.1, cast16==i) #data stores dat16.1, but separated by cast
#print(i)
last=length(data16[,1])


Depth16<-as.numeric(data16[,4])
Temp16<-as.numeric(data16[,6])
Salt16<-as.numeric(data16[,5])
Long16<-as.numeric(data16[,2])
Lat16<-as.numeric(data16[,3])
p16<-as.numeric(data16[,7])
SA16<-as.numeric(data16[,8])
pt16<-as.numeric(data16[,9])
CT16<-as.numeric(data16[,10])


#tail function goes to last value above. integer means it has no decimals vs. numeric
maxDepth16<-as.integer(tail(Depth16,1)) 
Dmax16 = maxDepth16

#Beginning the conditional statement for loop
if (last>=1)
{
  
  #Calculating and defining bin averages
  binDepth16<-seq(2,maxDepth16,1)#starts at 2m, goes to maxdepth, sorts by 1m
  
  binSalt16<-binMeans(Salt16, x=Depth16, bx=binDepth16, na.rm=TRUE)
  binTemp16<-binMeans(Temp16, x=Depth16, bx=binDepth16, na.rm=TRUE)
  binp16<-binMeans(p16, x=Depth16, bx=binDepth16, na.rm=TRUE)
  binSA16<-binMeans(SA16, x=Depth16, bx=binDepth16, na.rm=TRUE)
  binpt16<-binMeans(pt16, x=Depth16, bx=binDepth16, na.rm=TRUE)
  binCT16<-binMeans(CT16, x=Depth16, bx=binDepth16, na.rm=TRUE)
  
  ###makes binDepth equal length to binSalt, binTemp, etc. borrowed from      michael's code
  binDepth16 = head(binDepth16,-1)
  
  
  ###calculating the amount of freshwater at the surface. borrowed from sam's code.
  s_ref = 34.8
  s_ref_vec16 = c(rep(s_ref, length(binDepth16)))
  fresh16<-as.integer(sum(((s_ref_vec16-binSalt16)/s_ref_vec16),na.rm=TRUE),digits=4)
  #print(fresh16)
  
  ###creating dataframe for binvalues, adjusting their lengths so they all match
  n=maxDepth16-2
  
  ###calculating first stability max in terms of bin averages values
  N2 = gsw_Nsquared(binSA16, binCT16, binp16, Lat16)
  #print(N2)
  
  ###creating new dataframe with binned data and Nsquared
  dat16.2<- data.frame(binDepth16[1:n], binSalt16[1:n], binTemp16[1:n], Long16[1:n], Lat16[1:n], binp16[1:n], binSA16[1:n], binpt16[1:n], binCT16[1:n], N2$N2[1:n])
  
  ###some binDepths don't exist and produce NAs in date 2. returns dat2 without rows with NA
  dat16.2 = na.omit(dat16.2)
  
  ###define last as the number of rows in dat2 
  last = length(dat16.2[,1]) 
  
  ###use new letter for index (we already used "i"). we are ignoring the first 50m of the profile; we are looking for the (deep) atlantic layer
  j=40 
  
  ###starting a new conditional within the previous one
  ###(dat16.2[,10]) is N2
  
  while(dat16.2[j,10] >= dat16.2[j-1,10]) 
    
  {
    
    ###sorts through each row (binDepth). looking for max N2 (below 50m)
    j=j+1 
    
  }###end of "while" statement
  
  
  ###setting Tmax to an arbitrary low number
  Tmax= -99
  
  ###beginning another for loop (still within the for loop and if-statement)
  ###need new index "k"
  for(k in j:last)
  {
    ###looking at everything below Tmax
    if(dat16.2[k,8] >= Tmax)
    {Tmax = dat16.2[k,8]
    D_Tmax = dat16.2[k,1]
    SA_Tmax = dat16.2[k,7]
    Lat_Tmax16 = dat16.2[k,5]
    Long_Tmax16 = dat16.2[k,4]
    
    ###creating a new dataframe for this output
    dat16.3 <- data.frame(Tmax, D_Tmax, SA_Tmax, Lat_Tmax16, Long_Tmax16, i, 2016)
    
    }###end if-statment
  }###end (second) for-loop
  
  ###combining output for each for-loop iteration. makes dataframe with one row per cast
  y16<- rbind.data.frame(y16,dat16.3)
  
}###end if-statement
}###end (first) for-loop

y16<- y16[,c(6, 7, 1, 2, 3, 4, 5)]      #rearrange columns so cast is first

colnames(y16)<- c("cast", "year", "tmax", "d_tmax", "sa_tmax", "lat_tmax", "lon_tmax")


#create datafile of just "y" dataframe

#write.csv(y16,file = C:\Users\ingri\MAST467\Data\2016\"dat3.csv")

p17= gsw_p_from_z(-Omg17$dep, Omg17$lat)
SA17 = gsw_SA_from_SP(Omg17$sal, p17, Omg17$lat, Omg17$lon)
pt17 = gsw_pt_from_t(SA17, Omg17$tem, p17, 0)
CT17= gsw_CT_from_t(SA17, Omg17$tem, p17)

dat17.1<-data.frame(Omg17$cast, Omg17$lon, Omg17$lat, Omg17$dep, Omg17$sal, Omg17$tem, p17, SA17, pt17, CT17)

colnames(dat17.1) = c("cast17", "lon17", "lat17", "dep17", "sal17", "tem17", "p17", "SA17", "pt17", "CT17")

###Many steps are this chunk is borrowed from Sam's code
###Defining empty vectors/variables/etc. BEFORE for loop begins...
Dmax17 = 0
fresh17 = 0
dat_all = NULL
y17= NULL #y is an empty dataframe that we will store the output of each iteration for

###Looping through each cast to define bin averaged variables
for (i in 185:217) #i refers to cast
{data17<-dplyr::filter(dat17.1, cast17==i)
  #print(i)
  last=length(data17[,1])
  
  
  Depth17<-as.numeric(data17[,4])
  Temp17<-as.numeric(data17[,6])
  Salt17<-as.numeric(data17[,5])
  Long17<-as.numeric(data17[,2])
  Lat17<-as.numeric(data17[,3])
  p17<-as.numeric(data17[,7])
  SA17<-as.numeric(data17[,8])
  pt17<-as.numeric(data17[,9])
  CT17<-as.numeric(data17[,10])
  
  
#tail function goes to last value above. integer means it has no decimals vs. numeric
  maxDepth17<-as.integer(tail(Depth17,1)) 
  Dmax17 = maxDepth17
  
#Beginning the conditional statement for loop
  if (last>=1)
  {
    
#Calculating and defining bin averages
  binDepth17<-seq(2,maxDepth17,1)#starts at 2m, goes to maxdepth, sorts by 1m
  
  binSalt17<-binMeans(Salt17, x=Depth17, bx=binDepth17, na.rm=TRUE)
  binTemp17<-binMeans(Temp17, x=Depth17, bx=binDepth17, na.rm=TRUE)
  binp17<-binMeans(p17, x=Depth17, bx=binDepth17, na.rm=TRUE)
  binSA17<-binMeans(SA17, x=Depth17, bx=binDepth17, na.rm=TRUE)
  binpt17<-binMeans(pt17, x=Depth17, bx=binDepth17, na.rm=TRUE)
  binCT17<-binMeans(CT17, x=Depth17, bx=binDepth17, na.rm=TRUE)
  
###makes binDepth equal length to binSalt, binTemp, etc. borrowed from      michael's code
  binDepth17 = head(binDepth17,-1)
  
  
###calculating the amount of freshwater at the surface. borrowed from sam's code.
  s_ref = 34.8
  s_ref_vec17 = c(rep(s_ref, length(binDepth17)))
  fresh17<-as.integer(sum(((s_ref_vec17-binSalt17)/s_ref_vec17),na.rm=TRUE),digits=4)
  #print(fresh17)
  
###creating dataframe for binvalues, adjusting their lengths so they all match
  n=maxDepth17-2
  
###calculating first stability max in terms of bin averages values
N2 = gsw_Nsquared(binSA17, binCT17, binp17, Lat17)
#print(N2)

###creating new dataframe with binned data and Nsquared
dat17.2<- data.frame(binDepth17[1:n], binSalt17[1:n], binTemp17[1:n], Long17[1:n], Lat17[1:n], binp17[1:n], binSA17[1:n], binpt17[1:n], binCT17[1:n], N2$N2[1:n])

###some binDepths don't exist and produce NAs in date 2. returns dat2 without rows with NA
dat17.2 = na.omit(dat17.2)
  
###define last as the number of rows in dat2 
last = length(dat17.2[,1]) 

###use new letter for index (we already used "i"). we are ignoring the first 50m of the profile; we are looking for the (deep) atlantic layer
j=40 

###starting a new conditional within the previous one
###(dat17.2[,10]) is N2

while(dat17.2[j,10] >= dat17.2[j-1,10]) 
  
{
  
###sorts through each row (binDepth). looking for max N2 (below 50m)
  j=j+1 
  
}###end of "while" statement


###setting Tmax to an arbitrary low number
Tmax= -99

###beginning another for loop (still within the for loop and if-statement)
###need new index "k"
for(k in j:last)
{
  ###looking at everything below Tmax
  if(dat17.2[k,8] >= Tmax)
  {Tmax = dat17.2[k,8]
  D_Tmax = dat17.2[k,1]
  SA_Tmax = dat17.2[k,7]
  Lat_Tmax17 = dat17.2[k,5]
  Long_Tmax17 = dat17.2[k,4]
  
  ###creating a new dataframe for this output
  dat17.3 <- data.frame(Tmax, D_Tmax, SA_Tmax, Lat_Tmax17, Long_Tmax17, i, 2017)
  
  }###end if-statment
}###end (second) for-loop

###combining output for each for-loop iteration. makes dataframe with one row per cast
y17<- rbind.data.frame(y17,dat17.3)
  
  }###end if-statement
}###end (first) for-loop

## Warning in gsw_Nsquared(binSA17, binCT17, binp17, Lat17): gsw_Nsquared() has some NA inputs ... these are not handled well

## Warning in gsw_Nsquared(binSA17, binCT17, binp17, Lat17): gsw_Nsquared() has some NA inputs ... these are not handled well

## Warning in gsw_Nsquared(binSA17, binCT17, binp17, Lat17): gsw_Nsquared() has some NA inputs ... these are not handled well

## Warning in gsw_Nsquared(binSA17, binCT17, binp17, Lat17): gsw_Nsquared() has some NA inputs ... these are not handled well

## Warning in gsw_Nsquared(binSA17, binCT17, binp17, Lat17): gsw_Nsquared() has some NA inputs ... these are not handled well

###rearranging columns so cast is first
y17 <- y17[,c(6, 7, 1, 2, 3, 4, 5)]

###renaming columns
colnames(y17)<-c("cast", "year", "tmax", "d_tmax", "sa_tmax", "lat_tmax", "lon_tmax")

###creating datafile of just "y"dataframe
###write.csv(y17, file = "C:/Users/ingri/MAST467/Data/dat17.3.csv")

p19 = gsw_p_from_z(-Omg19$dep, Omg19$lat)
SA19 = gsw_SA_from_SP(Omg19$sal, p19, Omg19$lat, Omg19$lon)
pt19 = gsw_pt_from_t(SA19, Omg19$tem, p19, 0)
CT19= gsw_CT_from_t(SA19, Omg19$tem, p19)

dat19.1<-data.frame(Omg19$cast, Omg19$lon, Omg19$lat, Omg19$dep, Omg19$sal, Omg19$tem, p19, SA19, pt19, CT19)

colnames(dat19.1) = c("cast19", "lon19", "lat19", "dep19", "sal19", "tem19", "p19", "SA19", "pt19", "CT19")

###Many steps are this chunk is borrowed from Sam's code
###Defining empty vectors/variables/etc. BEFORE for loop begins...
Dmax19 = 0
fresh19 = 0
dat_all = NULL
y19= NULL #y is an empty dataframe that we will store the output of each iteration for

###Looping through each cast to define bin averaged variables
for (i in 106:147) #i refers to cast
{data19<-dplyr::filter(dat19.1, cast19==i)

  print(i)
  last=length(data19[,1])
  
  
  Depth19<-as.numeric(data19[,4])
  Temp19<-as.numeric(data19[,6])
  Salt19<-as.numeric(data19[,5])
  Long19<-as.numeric(data19[,2])
  Lat19<-as.numeric(data19[,3])
  p19<-as.numeric(data19[,7])
  SA19<-as.numeric(data19[,8])
  pt19<-as.numeric(data19[,9])
  CT19<-as.numeric(data19[,10])
  
  
#tail function goes to last value above. integer means it has no decimals vs. numeric
  maxDepth19<-as.integer(tail(Depth19,1)) 
  Dmax19 = maxDepth19
  
#Beginning the conditional statement for loop
  if (last>=1)
  {
    
#Calculating and defining bin averages
  binDepth19<-seq(2,maxDepth19,1)#starts at 2m, goes to maxdepth, sorts by 1m
  
  binSalt19<-binMeans(Salt19, x=Depth19, bx=binDepth19, na.rm=TRUE)
  binTemp19<-binMeans(Temp19, x=Depth19, bx=binDepth19, na.rm=TRUE)
  binp19<-binMeans(p19, x=Depth19, bx=binDepth19, na.rm=TRUE)
  binSA19<-binMeans(SA19, x=Depth19, bx=binDepth19, na.rm=TRUE)
  binpt19<-binMeans(pt19, x=Depth19, bx=binDepth19, na.rm=TRUE)
  binCT19<-binMeans(CT19, x=Depth19, bx=binDepth19, na.rm=TRUE)
  
###makes binDepth equal length to binSalt, binTemp, etc. borrowed from      michael's code
  binDepth19 = head(binDepth19,-1)
  
  
###calculating the amount of freshwater at the surface. borrowed from sam's code.
  s_ref = 34.8
  s_ref_vec19 = c(rep(s_ref, length(binDepth19)))
  fresh19<-as.integer(sum(((s_ref_vec19-binSalt19)/s_ref_vec19),na.rm=TRUE),digits=4)
  print(fresh19)
  
###creating dataframe for binvalues, adjusting their lengths so they all match
  n=maxDepth19-2
  
###calculating first stability max in terms of bin averages values
N2 = gsw_Nsquared(binSA19, binCT19, binp19, Lat19)
#print(N2)

###creating new dataframe with binned data and Nsquared
dat19.2<- data.frame(binDepth19[1:n], binSalt19[1:n], binTemp19[1:n], Long19[1:n], Lat19[1:n], binp19[1:n], binSA19[1:n], binpt19[1:n], binCT19[1:n], N2$N2[1:n])

###some binDepths don't exist and produce NAs in date 2. returns dat2 without rows with NA
dat19.2 = na.omit(dat19.2)
  
###define last as the number of rows in dat2 
last = length(dat19.2[,1]) 

###use new letter for index (we already used "i"). we are ignoring the first 50m of the profile; we are looking for the (deep) atlantic layer
j=50 

###starting a new conditional within the previous one
###(dat19.2[,10]) is N2
while(dat19.2[j,10] >= dat19.2[j-1,10]) 
{
###sorts through each row (binDepth). looking for max N2 (below 50m)
  j=j+1 
  
}###end of "while" statement


###setting Tmax to an arbitrary low number
Tmax= -99

###beginning another for loop (still within the for loop and if-statement)
###need new index "k"
for(k in j:last)
{
  ###looking at everything below Tmax
  if(dat19.2[k,8] >= Tmax)
  {Tmax = dat19.2[k,8]
  D_Tmax = dat19.2[k,1]
  SA_Tmax = dat19.2[k,7]
  Lat_Tmax19 = dat19.2[k,5]
  Long_Tmax19 = dat19.2[k,4]
  
  ###creating a new dataframe for this output
  dat19.3 <- data.frame(Tmax, D_Tmax, SA_Tmax, Lat_Tmax19, Long_Tmax19, i, 2019)
  
  }###end if-statment
}###end (second) for-loop

###combining output for each for-loop iteration. makes dataframe with one row per cast
y19<- rbind.data.frame(y19,dat19.3)
  
  }###end if-statement
}###end (first) for-loop

## [1] 106
## [1] 5
## [1] 107
## [1] 10
## [1] 108
## [1] 10
## [1] 109
## [1] 5

## Warning in gsw_Nsquared(binSA19, binCT19, binp19, Lat19): gsw_Nsquared() has some NA inputs ... these are not handled well

## [1] 110
## [1] 12
## [1] 111
## [1] 112
## [1] 113
## [1] 114
## [1] 115
## [1] 116
## [1] 117
## [1] 118
## [1] 119
## [1] 120
## [1] 121
## [1] 122
## [1] 123
## [1] 124
## [1] 125
## [1] 126
## [1] 127
## [1] 128
## [1] 129
## [1] 130
## [1] 131
## [1] 132
## [1] 133
## [1] 134
## [1] 135
## [1] 136
## [1] 137
## [1] 138
## [1] 139
## [1] 140
## [1] 141
## [1] 142
## [1] 143
## [1] 144
## [1] 145
## [1] 146
## [1] 8

## Warning in gsw_Nsquared(binSA19, binCT19, binp19, Lat19): gsw_Nsquared() has some NA inputs ... these are not handled well

## [1] 147
## [1] 9

###rearranging columns so cast is first
y19 <- y19[,c(6, 7, 1, 2, 3, 4, 5)]

###renaming columns
colnames(y19)<-c("cast", "year", "tmax", "d_tmax", "sa_tmax", "lat_tmax", "lon_tmax")

###creating datafile of just "y"dataframe
###write.csv(y19, file = "C:/Users/ingri/MAST467/Data/dat19.3.csv")

p20 = gsw_p_from_z(-Omg20$dep, Omg20$lat)
SA20 = gsw_SA_from_SP(Omg20$sal, p20, Omg20$lat, Omg20$lon)
pt20 = gsw_pt_from_t(SA20, Omg20$tem, p20, 0)
CT20= gsw_CT_from_t(SA20, Omg20$tem, p20)

dat20.1<-data.frame(Omg20$cast, Omg20$lon, Omg20$lat, Omg20$dep, Omg20$sal, Omg20$tem, p20, SA20, pt20, CT20)

colnames(dat20.1) = c("cast20", "lon20", "lat20", "dep20", "sal20", "tem20", "p20", "SA20", "pt20", "CT20")

###Many steps are this chunk is borrowed from Sam's code
###Defining empty vectors/variables/etc. BEFORE for loop begins...
Dmax20 = 0
fresh20 = 0
dat_all = NULL
y20= NULL #y is an empty dataframe that we will store the output of each iteration for

###Looping through each cast to define bin averaged variables
for (i in 1:297) #i refers to cast
{data20<-dplyr::filter(dat20.1, cast20==i)
  #print(i)
  last=length(data20[,1])
  
  
  Depth20<-as.numeric(data20[,4])
  Temp20<-as.numeric(data20[,6])
  Salt20<-as.numeric(data20[,5])
  Long20<-as.numeric(data20[,2])
  Lat20<-as.numeric(data20[,3])
  p20<-as.numeric(data20[,7])
  SA20<-as.numeric(data20[,8])
  pt20<-as.numeric(data20[,9])
  CT20<-as.numeric(data20[,10])
  
  
#tail function goes to last value above. integer means it has no decimals vs. numeric
  maxDepth20<-as.integer(tail(Depth20,1)) 
  Dmax20 = maxDepth20
  
#Beginning the conditional statement for loop
  if (last>=1)
  {
    
#Calculating and defining bin averages
  binDepth20<-seq(2,maxDepth20,1)#starts at 2m, goes to maxdepth, sorts by 1m
  
  binSalt20<-binMeans(Salt20, x=Depth20, bx=binDepth20, na.rm=TRUE)
  binTemp20<-binMeans(Temp20, x=Depth20, bx=binDepth20, na.rm=TRUE)
  binp20<-binMeans(p20, x=Depth20, bx=binDepth20, na.rm=TRUE)
  binSA20<-binMeans(SA20, x=Depth20, bx=binDepth20, na.rm=TRUE)
  binpt20<-binMeans(pt20, x=Depth20, bx=binDepth20, na.rm=TRUE)
  binCT20<-binMeans(CT20, x=Depth20, bx=binDepth20, na.rm=TRUE)
  
###makes binDepth equal length to binSalt, binTemp, etc. borrowed from      michael's code
  binDepth20 = head(binDepth20,-1)
  
  
###calculating the amount of freshwater at the surface. borrowed from sam's code.
  s_ref = 34.8
  s_ref_vec20 = c(rep(s_ref, length(binDepth20)))
  fresh20<-as.integer(sum(((s_ref_vec20-binSalt20)/s_ref_vec20),na.rm=TRUE),digits=4)
  #print(fresh20)
  
###creating dataframe for binvalues, adjusting their lengths so they all match
  n=maxDepth20-2
  
###calculating first stability max in terms of bin averages values
N2 = gsw_Nsquared(binSA20, binCT20, binp20, Lat20)
#print(N2)

###creating new dataframe with binned data and Nsquared
dat20.2<- data.frame(binDepth20[1:n], binSalt20[1:n], binTemp20[1:n], Long20[1:n], Lat20[1:n], binp20[1:n], binSA20[1:n], binpt20[1:n], binCT20[1:n], N2$N2[1:n])

###some binDepths don't exist and produce NAs in date 2. returns dat2 without rows with NA
dat20.2 = na.omit(dat20.2)
  
###define last as the number of rows in dat2 
last = length(dat20.2[,1]) 

###use new letter for index (we already used "i"). we are ignoring the first 50m of the profile; we are looking for the (deep) atlantic layer
j=40 

###starting a new conditional within the previous one
###(dat20.2[,10]) is N2
while(dat20.2[j,10] >= dat20.2[j-1,10]) 
{
###sorts through each row (binDepth). looking for max N2 (below 50m)
  j=j+1 
  
}###end of "while" statement


###setting Tmax to an arbitrary low number
Tmax= -99

###beginning another for loop (still within the for loop and if-statement)
###need new index "k"
for(k in j:last)
{
  ###looking at everything below Tmax
  if(dat20.2[k,8] >= Tmax)
  {Tmax = dat20.2[k,8]
  D_Tmax = dat20.2[k,1]
  SA_Tmax = dat20.2[k,7]
  Lat_Tmax20 = dat20.2[k,5]
  Long_Tmax20 = dat20.2[k,4]
  
  ###creating a new dataframe for this output
  dat20.3 <- data.frame(Tmax, D_Tmax, SA_Tmax, Lat_Tmax20, Long_Tmax20, i, 2020)
  
  }###end if-statment
}###end (second) for-loop

###combining output for each for-loop iteration. makes dataframe with one row per cast
y20<- rbind.data.frame(y20,dat20.3)
  
  }###end if-statement
}###end (first) for-loop

## Warning in gsw_Nsquared(binSA20, binCT20, binp20, Lat20): gsw_Nsquared() has some NA inputs ... these are not handled well

## Warning in gsw_Nsquared(binSA20, binCT20, binp20, Lat20): gsw_Nsquared() has some NA inputs ... these are not handled well

## Warning in gsw_Nsquared(binSA20, binCT20, binp20, Lat20): gsw_Nsquared() has some NA inputs ... these are not handled well

## Warning in gsw_Nsquared(binSA20, binCT20, binp20, Lat20): gsw_Nsquared() has some NA inputs ... these are not handled well

###rearranging columns so cast is first
y20 <- y20[,c(6, 7, 1, 2, 3, 4, 5)]

###renaming columns
colnames(y20)<-c("cast", "year", "tmax", "d_tmax", "sa_tmax", "lat_tmax", "lon_tmax")

###creating datafile of just "y"dataframe
###write.csv(y20, file = "C:/Users/ingri/MAST467/Data/dat20.3.csv")

p21 = gsw_p_from_z(-Omg21$dep, Omg21$lat)
SA21 = gsw_SA_from_SP(Omg21$sal, p21, Omg21$lat, Omg21$lon)
pt21 = gsw_pt_from_t(SA21, Omg21$tem, p21, 0)
CT21= gsw_CT_from_t(SA21, Omg21$tem, p21)


dat21.1<-data.frame(Omg21$cast, Omg21$lon, Omg21$lat, Omg21$dep, Omg21$sal, Omg21$tem, p21, SA21, pt21, CT21)

colnames(dat21.1) = c("cast21", "lon21", "lat21", "dep21", "sal21", "tem21", "p21", "SA21", "pt21", "CT21")

###Many steps are this chunk is borrowed from Sam's code
###Defining empty vectors/variables/etc. BEFORE for loop begins...
Dmax21 = 0
fresh21 = 0
dat_all = NULL
y21= NULL #y is an empty dataframe that we will store the output of each iteration for

###Looping through each cast to define bin averaged variables
for (i in 1:293) #i refers to cast
{data21<-dplyr::filter(dat21.1, cast21==i)
  #print(i)
  last=length(data21[,1])
  
  
  Depth21<-as.numeric(data21[,4])
  Temp21<-as.numeric(data21[,6])
  Salt21<-as.numeric(data21[,5])
  Long21<-as.numeric(data21[,2])
  Lat21<-as.numeric(data21[,3])
  p21<-as.numeric(data21[,7])
  SA21<-as.numeric(data21[,8])
  pt21<-as.numeric(data21[,9])
  CT21<-as.numeric(data21[,10])
  
  
#tail function goes to last value above. integer means it has no decimals vs. numeric
  maxDepth21<-as.integer(tail(Depth21,1)) 
  Dmax21 = maxDepth21
  
#Beginning the conditional statement for loop
  if (last>=1)
  {
    
#Calculating and defining bin averages
  binDepth21<-seq(2,maxDepth21,1)#starts at 2m, goes to maxdepth, sorts by 1m
  
  binSalt21<-binMeans(Salt21, x=Depth21, bx=binDepth21, na.rm=TRUE)
  binTemp21<-binMeans(Temp21, x=Depth21, bx=binDepth21, na.rm=TRUE)
  binp21<-binMeans(p21, x=Depth21, bx=binDepth21, na.rm=TRUE)
  binSA21<-binMeans(SA21, x=Depth21, bx=binDepth21, na.rm=TRUE)
  binpt21<-binMeans(pt21, x=Depth21, bx=binDepth21, na.rm=TRUE)
  binCT21<-binMeans(CT21, x=Depth21, bx=binDepth21, na.rm=TRUE)
  
###makes binDepth equal length to binSalt, binTemp, etc. borrowed from      michael's code
  binDepth21 = head(binDepth21,-1)
  
  
###calculating the amount of freshwater at the surface. borrowed from sam's code.
  s_ref = 34.8
  s_ref_vec21 = c(rep(s_ref, length(binDepth21)))
  fresh21<-as.integer(sum(((s_ref_vec21-binSalt21)/s_ref_vec21),na.rm=TRUE),digits=4)
  #print(fresh21)
  
###creating dataframe for binvalues, adjusting their lengths so they all match
  n=maxDepth21-2
  
###calculating first stability max in terms of bin averages values
N2 = gsw_Nsquared(binSA21, binCT21, binp21, Lat21)
#print(N2)

###creating new dataframe with binned data and Nsquared
dat21.2<- data.frame(binDepth21[1:n], binSalt21[1:n], binTemp21[1:n], Long21[1:n], Lat21[1:n], binp21[1:n], binSA21[1:n], binpt21[1:n], binCT21[1:n], N2$N2[1:n])

###some binDepths don't exist and produce NAs in date 2. returns dat2 without rows with NA
dat21.2 = na.omit(dat21.2)
  
###define last as the number of rows in dat2 
last = length(dat21.2[,1]) 

###use new letter for index (we already used "i"). we are ignoring the first 50m of the profile; we are looking for the (deep) atlantic layer
j=40 

###starting a new conditional within the previous one
###(dat21.2[,10]) is N2
while(dat21.2[j,10] >= dat21.2[j-1,10]) 
{
###sorts through each row (binDepth). looking for max N2 (below 50m)
  j=j+1 
  
}###end of "while" statement


###setting Tmax to an arbitrary low number
Tmax= -99

###beginning another for loop (still within the for loop and if-statement)
###need new index "k"
for(k in j:last)
{
  ###looking at everything below Tmax
  if(dat21.2[k,8] >= Tmax)
  {Tmax = dat21.2[k,8]
  D_Tmax = dat21.2[k,1]
  SA_Tmax = dat21.2[k,7]
  Lat_Tmax21 = dat21.2[k,5]
  Long_Tmax21 = dat21.2[k,4]
  
  ###creating a new dataframe for this output
  dat21.3 <- data.frame(Tmax, D_Tmax, SA_Tmax, Lat_Tmax21, Long_Tmax21, i, 2021)
  
  }###end if-statment
}###end (second) for-loop

###combining output for each for-loop iteration. makes dataframe with one row per cast
y21<- rbind.data.frame(y21,dat21.3)
  
  }###end if-statement
}###end (first) for-loop

## Warning in gsw_Nsquared(binSA21, binCT21, binp21, Lat21): gsw_Nsquared() has some NA inputs ... these are not handled well

###rearranging columns so cast is first
y21 <- y21[,c(6, 7, 1, 2, 3, 4, 5)]

###renaming columns
colnames(y21)<-c("cast", "year", "tmax", "d_tmax", "sa_tmax", "lat_tmax", "lon_tmax")

###creating datafile of just "y"dataframe
###write.csv(y21, file = "C:/Users/ingri/MAST467/Data/dat21.3.csv")

###Using while loop to identify the surface mixed layer (SML)
sst <-as.numeric(dat19.2[1,3])
print(sst)

## [1] 7.382222

  i = 2
  while(dat19.2[i,3]  <= (sst+.2) && dat19.2[i,3] >= (sst-.2))
{print(dat19.2[i,3])
    
   
    i = i+1
    
    
}###end loop

## [1] 7.371111
## [1] 7.37
## [1] 7.37
## [1] 7.37
## [1] 7.351111
## [1] 7.338889
## [1] 7.326
## [1] 7.314444
## [1] 7.285556
## [1] 7.27
## [1] 7.264444
## [1] 7.25
## [1] 7.237778
## [1] 7.211111
## [1] 7.186667

  print(i)

## [1] 17

  z_SML <- binDepth19[16]
  print(z_SML)

## [1] 17

#r plotting-binDepth-binSalt-and-binTemp-with-SML
###calculate the mean value and standard deviation of potential temperature and absolute salinity 
  mean(pt19[1:z_SML])

## [1] 7.376787

  sd(pt19[1:z_SML])

## [1] 0.007761541

  mean(SA19[1:z_SML])

## [1] 33.00729

  sd(SA19[1:z_SML])

## [1] 0.03899222

plot(binpt19, binDepth19, ylim = rev(range(binDepth19)), xlab = "potential temperature (C)", ylab = "depth (m)", col="black", lwd = 2, type="l")
points(binpt19[z_SML], binDepth19[z_SML], col = "red", lwd=2, pch = 20)

plot(binSA19, binDepth19, ylim = rev(range(binDepth19)), xlab = "Absolute Salinity (ppt)", ylab = "depth (m)", col="black", lwd = 2, type="l")
points(binSA19[z_SML], binDepth19[z_SML], col = "red", lwd=2, pch = 20)

###Borrrowed from Sam's code
#Combine all the "y[yr]" dataframes into one large one

d_all<- rbind(y16,y17,y19,y20,y21)
#write.csv(d_all,file = "C:\Users\ingri\MAST467\Data"\all_dat.csv")

a = mean(y16[,3])
b = mean(y17[,3])
c = mean(y19[,3])
d = mean(y20[,3])
e = mean(y21[,3])

T_avg = rbind(a,b,c,d,e)
rm(a,b,c,d,e)

Year = c(2016,2017,2019,2020,2021)

a = mean(y16[,4])
b = mean(y17[,4])
c = mean(y19[,4])
d = mean(y20[,4])
e = mean(y21[,4])

D_Tavg = rbind(a,b,c,d,e)
rm(a,b,c,d,e)

a = mean(y16[,5])
b = mean(y17[,5])
c = mean(y19[,5])
d = mean(y20[,5])
e = mean(y21[,5])

SA_Tavg = rbind(a,b,c,d,e)
rm(a,b,c,d,e)

avg_df = data.frame(Year, T_avg, D_Tavg, SA_Tavg)

###PLOTTING
ggplot(d_all) + 
  geom_point(aes(x=year, y=tmax))+
  stat_summary(aes(x = year,y = tmax,group=1), fun=median, color="red", geom="line",group=1)+
  ggtitle("Median T Max for Region 5 (2016, 2017, 2019-2021)")+
  stat_summary(aes(x = year,y = tmax,group=1), fun=median, color="red", geom="point",group=1, size = 3)+
  theme(legend.position = "none")

ggplot(d_all) + 
  geom_point(aes(x=year, y=d_tmax))+
  stat_summary(aes(x = year,y = d_tmax,group=1), fun=median, color="red", geom="line",group=1)+
  ggtitle("Median Depth of TMax for Region 5 (2016, 2017, 2019-2021)")+
  stat_summary(aes(x = year,y = d_tmax,group=1), fun=median, color="red", geom="point",group=1, size = 3)+
  theme(legend.position = "none")

ggplot(d_all) + 
  geom_point(aes(x=year, y=sa_tmax))+
  stat_summary(aes(x = year,y = sa_tmax,group=1), fun=median, color="red", geom="line",group=1)+
  ggtitle("Median Salinity of TMax for Region 5 (2016, 2017, 2019-2021)")+
  stat_summary(aes(x = year,y = sa_tmax,group=1), fun=median, color="red", geom="point",group=1, size = 3)+
  theme(legend.position = "none")