OMG00

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:

echo = F
library(devtools)

## Warning: package 'devtools' was built under R version 4.0.5

## Loading required package: usethis

## Warning: package 'usethis' was built under R version 4.0.5

library(dplyr)

## Warning: package 'dplyr' was built under R version 4.0.5

## 
## 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.0.5

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.0.5

library(leaflet)

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

library(devtools)
library(gsw)

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

#install_github("MikkoVihtakari/PlotSvalbard")
library(PlotSvalbard)
library(gsw)

echo = F
#shell("cd c:/Users/Saman/Onedrive/Desktop/Mast467/Data/", translate=TRUE)
#shell("dir")
#shell("c:/Users/Saman/OneDrive/Desktop/Mast467/Data/nasa.cmd",translate=TRUE)
Omg<- read.csv("c:/Users/Saman/Onedrive/Desktop/Mast467/Data/output.dat",sep=" ")

#Calculate Values using gsw package
#library(gsw)

# Pressure from depth
P_1 = gsw_p_from_z(-Omg$dep,Omg$lat)
#absolute salinity from practical salinity
SA_1 = gsw_SA_from_SP(Omg$sal,P_1,Omg$lon,Omg$lat)
#potential temperature from in situ t
PT_1 = gsw_pt_from_t(SA_1,Omg$temp,P_1,0)
#conservative temp from pt
CT_1 = gsw_CT_from_t(SA_1,Omg$temp,P_1)
#potential density at p=0
#r0_1 = gsw_sigma0(SA_1,CT_1)
#potential density with p=1000
#r1_1 = gsw_sigma1(SA_1,CT_1)
# Freezing Point
#This code gets an error
#tf = gsw_CT_freezing_poly(SA,p,1)

#Put OMG data into data frame
#dat1 = data.frame(Omg$lon, Omg$lat, Omg$dep, Omg$temp, Omg$sal, Omg$cast)

dat1 <- data.frame(Omg$cast, Omg$lon, Omg$lat, Omg$dep, Omg$sal, Omg$temp, P_1, SA_1, PT_1, CT_1)

#names(dat1) = NULL
#rename columns
#colnames(omgdf) = c("Cast", "Lon", "Lat", "Dep", "Sal", "Temp")

colnames(dat1) = c("cast", "Lon", "Lat", "Dep", "Sal", "Temp", "P", "SA", "PT", "CT")

#Isolate a single cast (156)
#C156 = dat1[dat1$cast==156,]

#write.csv(dat1,file = "c:/Users/Saman/Onedrive/Desktop/Mast467/Data/dat1.csv")

#Install package matrixStats
#Load it
library(matrixStats)

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

## 
## Attaching package: 'matrixStats'

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

library(dplyr)
#Create forloop to print maximum depth for each cast
Dmax=0
frsh=0
#fresh = vector()
y= NULL
for(i in 127:157)
#for (i in 144:152)
{data<-dplyr::filter(dat1,cast==i)  
  print(i)                         #prints cast number
last=length(data[,1])              

 Depth<- as.numeric(data[,4])
 Temp<- as.numeric(data[,6])
 Salt<- as.numeric(data[,5])
 Lat<- as.numeric(data[,3])
 Lon<- as.numeric(data[,2])
 P<- as.numeric(data[,7])
 SA<- as.numeric(data[,8])
 PT<- as.numeric(data[,9])
 CT<- as.numeric(data[,10])
 
 maxDepth<- as.integer(tail(Depth,1)) #tail function goes to last value
 print(maxDepth)              #prints max depth
 
 Dmax=maxDepth                
 
 if(last >=1)           #only proceeds if cast has data

    {
   
 binDepth<-seq(2,maxDepth,1)  #starts at 2m, goes to maxdepth, sorts by 1m
 binSalt<-binMeans(Salt,x=Depth,bx=binDepth,na.rm=TRUE)
 binTemp<-binMeans(Temp,x=Depth,bx=binDepth,na.rm=TRUE)
 binSA<-binMeans(SA,x=Depth,bx=binDepth,na.rm=TRUE)
 binPT<-binMeans(PT,x=Depth,bx=binDepth,na.rm=TRUE) 
 binP<-binMeans(P, x=Depth, bx=binDepth, na.rm=TRUE)
 binCT<-binMeans(CT, x=Depth, bx=binDepth, na.rm=TRUE)
 binLat<-binMeans(Lat, x=Depth, bx=binDepth, na.rm=TRUE)
 
# n = maxDepth-2
# binDepth<-binDepth[1:n]
 
 binDepth<-head(binDepth,-1)  #makes binDepth equal length to binSalt,etc.
      #I borrowed this from Michael's code
 
  n = maxDepth -2

N2<- gsw_Nsquared(binSA, binCT, binP, binLat)
#if (N2$N2[1:n] <= 0) {N2$N2[1:n] = 0}
 
dat2<-data.frame(binDepth[1:n],binSalt[1:n],binTemp[1:n],Lon[1:n],Lat[1:n], binP[1:n], binSA[1:n],binPT[1:n],binCT[1:n], N2$N2[1:n])

dat2 = na.omit(dat2) # Remove rows with NA's

colnames(dat2) = c("binD", "binS", "binT", "Lon", "Lat", "binP", "binSA", "binPT", "binCT", "N2")



  last = length(dat2[,1])
  j=50
  while(dat2[j,10] >= dat2[j-1,10])
  {#print(dat2[j,10])
    
    j=j+1
  }
  
  Tmax = -99
  
  for(k in j:last)
  {
    if( dat2[k,8] >= Tmax)
      {Tmax = dat2[k,8]
      D_Tmax = dat2[k,1]
      S_Tmax = dat2[k,7]
      Lat_Tmax = dat2[k,5]
      Lon_Tmax = dat2[k,4]
      
      dat3 <- data.frame(Tmax, D_Tmax, S_Tmax, Lat_Tmax, Lon_Tmax, i)
      #colnames(dat3) = c("Tmax", "D", "S", "Lat", "Lon")
    
    }}
  y<- rbind.data.frame(y, dat3)
 } 
 
}    # END LOOP

## [1] 127
## [1] 450
## [1] 128
## integer(0)
## [1] 129
## [1] 480

## Warning in gsw_Nsquared(binSA, binCT, binP, binLat): gsw_Nsquared() has some NA inputs ... these are not handled well

## [1] 130
## [1] 338
## [1] 131
## [1] 386
## [1] 132
## integer(0)
## [1] 133
## integer(0)
## [1] 134
## integer(0)
## [1] 135
## integer(0)
## [1] 136
## [1] 337

## Warning in gsw_Nsquared(binSA, binCT, binP, binLat): gsw_Nsquared() has some NA inputs ... these are not handled well

## [1] 137
## [1] 362
## [1] 138
## [1] 502
## [1] 139
## integer(0)
## [1] 140
## [1] 718

## Warning in gsw_Nsquared(binSA, binCT, binP, binLat): gsw_Nsquared() has some NA inputs ... these are not handled well

## [1] 141
## [1] 641
## [1] 142
## [1] 531
## [1] 143
## [1] 445

## Warning in gsw_Nsquared(binSA, binCT, binP, binLat): gsw_Nsquared() has some NA inputs ... these are not handled well

## [1] 144
## [1] 453
## [1] 145
## [1] 342
## [1] 146
## [1] 467
## [1] 147
## [1] 461
## [1] 148
## [1] 417
## [1] 149
## [1] 447
## [1] 150
## [1] 427
## [1] 151
## integer(0)
## [1] 152
## [1] 315
## [1] 153
## integer(0)
## [1] 154
## [1] 278
## [1] 155
## integer(0)
## [1] 156
## [1] 173
## [1] 157
## [1] 275

y<- y[,c(6,1,2,3,4,5)]
colnames(y)<- c("Cast", "Tmax", "D_Tmax", "S_Tmax", "Lat", "Long")

write.csv(y,file = "c:/Users/Saman/Onedrive/Desktop/Mast467/Data/dat3.csv")

#write.csv(dat2,"c:/Users/Saman/Onedrive/Desktop/Mast467/Data/dat2.csv") 

 
#ref<-34.8                              #set reference salinity
#s_ref<- c(rep(ref,length(binDepth)))   #make ref sal a vector with the                                                           #samelength as bin depth
  
#fresh[i] = round(sum(((s_ref - binSalt)/s_ref),na.rm = TRUE),digits=4)

#frsh<-as.integer(sum(((s_ref-binSalt)/s_ref),na.rm=TRUE),digits=4) #frsh is the sum of all                                                   #freshwater in the bin
#print(frsh)      #print freshwater


#write.csv(dat2,"c:/Users/Saman/Onedrive/Desktop/Mast467/Data/dat2.csv") 

#colnames(dat2) = c("binD", "binS", "binT", "Lon", "Lat", "binP", "binSA", "binPT", "binCT", "N2")



#output should be cast number, max depth, freshwater

#create dataframe for binvalues, adjust their lengths so they all match

#n=maxDepth-2

#N2<- gsw_Nsquared(binSA, binCT, binP, Lat)

#dat2<-data.frame(binDepth[1:n],binSalt[1:n],binTemp[1:n],Lon[1:n],Lat[1:n], binP[1:n], binSA[1:n],binPT[1:n],binCT[1:n], N2$N2[1:n])

#colnames(dat2) = c("binD", "binS", "binT", "Lon", "Lat", "binP", "binSA", "binPT", "binCT", "N2")

#graphing 

require(fields)

## Loading required package: fields

## Warning: package 'fields' was built under R version 4.0.5

## Loading required package: spam

## Warning: package 'spam' was built under R version 4.0.5

## Loading required package: dotCall64

## Warning: package 'dotCall64' was built under R version 4.0.5

## Loading required package: grid

## Spam version 2.7-0 (2021-06-25) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction 
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.

## 
## Attaching package: 'spam'

## The following objects are masked from 'package:base':
## 
##     backsolve, forwardsolve

## Loading required package: viridis

## Warning: package 'viridis' was built under R version 4.0.5

## Loading required package: viridisLite

## Warning: package 'viridisLite' was built under R version 4.0.5

## 
## Try help(fields) to get started.

## 
## Attaching package: 'fields'

## The following object is masked from 'package:leaflet':
## 
##     addLegend

plot(binSalt, binDepth, ylim = rev(range(binDepth)), xlab = "Salinity (ppt)", ylab = "Depth (m)", type = "l")

plot(data$Sal, data$Dep, ylim = rev(range(data$Dep)),xlim = range(data$Sal), xlab = "salinity (ppt)", ylab = "depth (m)", col="darkred", lwd = 2, type="l")
points(binSalt, binDepth, col="steelblue", lwd = 2, type = "l")
legend(x=33, y=400, legend=c("bin average", "original"), fill = c(4,2))

plot(data$Temp, data$Dep, ylim = rev(range(data$Dep)), xlab = "temperature (C)", ylab = "depth (m)", col="darkred", lwd = 4, type="l")
points(binTemp, binDepth, col="steelblue", lwd = 2, type = "l")
legend(x=5, y=500, legend=c("bin average", "original"), fill = c(4,2))

library(gsw)

sst<-as.numeric(dat2[1,3])
up<- sst + 0.2
low<- sst - 0.2
print(sst)

## [1] 1.348889

#last=length(dat2)
  i=2
  while(dat2[i,3] <=up && dat2[i,3]>=low) 
  {
    print(dat2[i,3])
    
 
    i=i+1

  }

## [1] 1.353333
## [1] 1.258889
## [1] 1.237778

z_SML<-binDepth[i]
print(z_SML)

## [1] 6

mean(PT[1:z_SML])

## [1] 1.376661

sd(PT[1:z_SML])

## [1] 0.035592

mean(SA[1:z_SML])

## [1] 13.05941

sd(SA[1:z_SML])

## [1] 0.4206699

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

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

echo = F
## Plot sal vs. depth for range
library(ggplot2)

#Individual profiles
ggplot(Omg, mapping=aes (x=dep, y=sal, group=cast)) +
  geom_path() +
  ylim(c(30,35)) +
  xlim(c(720,0)) +
  ylab("Salinity, psu") +
  xlab("Depth, m") +
  coord_flip(expand=FALSE) + 
  facet_wrap(~cast)

## Warning: Removed 1898 row(s) containing missing values (geom_path).

#Profiles on one plot
ggplot(Omg, mapping=aes (x=dep, y=sal, group=cast)) +
  geom_path() +
  xlim(c(720,0)) +
  ylab("Salinity, psu") +
  xlab("Depth, m") +
  coord_flip(expand=FALSE) 

## Plot temp vs. depth for individual casts
library(ggplot2)
ggplot(Omg, mapping=aes (x=dep, y=temp, group=cast)) +
  geom_path() +
  #ylim(c(-3.0,6.5)) +
  xlim(c(720,0)) +
  ylab("Temperature, C") +
  xlab("Depth, m") +
  coord_flip(expand=FALSE) +
  facet_wrap(~cast)

#Plot temp vs. depth for all casts
ggplot(Omg, mapping=aes (x=dep, y=temp, group=cast)) +
  geom_path() +
  #ylim(c(-3.0,6.5)) +
  xlim(c(720,0)) +
  ylab("Temperature, C") +
  xlab("Depth, m") +
  coord_flip(expand=FALSE) 

#Salinity and Temperature (all casts)
ggplot(Omg,mapping=aes(x=sal,y=temp,group=cast,color=dep)) +
  geom_point() +
  xlim(c(30,35)) +
  ylim(c(-2,4)) +
  xlab("Salinity, psu") +
  ylab("Temperature, C") +
  scale_color_viridis_c(guide=guide_colorbar(barheight=15,reverse=TRUE)) #+

## Warning: Removed 1911 rows containing missing values (geom_point).

  #facet_wrap(~cast)

#Isolate a single cast (156)
#C156 = dat1[dat1$cast==156,]

#Create temperature profile for this single cast
#ggplot(C156, mapping=aes (x=z, y=t)) +
 # geom_path() +
 # xlim(c(200,0)) +
#  ylab("Temperature, C") +
#  xlab("Depth, m") +
#  coord_flip(expand=FALSE) +
#  ggtitle("Temperature Profile Cast 156")

#Create Salinity profile for single cast
#ggplot(C156, mapping=aes (x=z, y=s)) +
#  geom_path() +
#  xlim(c(200,0)) +
#  ylab("Salinity, psu") +
#  xlab("Depth, m") +
#  coord_flip(expand=FALSE) +
#  ggtitle("Salinity Profile Cast 156")

#New salinity profile for depths between 150 and 180 m
#ggplot(C156, mapping=aes (x=z, y=s)) +
#  geom_path() +
#  xlim(c(180,150)) +
#  ylim(c(34.2,34.6)) +
#  ylab("Salinity, psu") +
#  xlab("Depth, m") +
#  coord_flip(expand=FALSE) +
#  ggtitle("Salinity Profile Cast 156")

#Graph another T profile for this same depth range
#ggplot(C156, mapping=aes (x=z, y=t)) +
#  geom_path() +
#  xlim(c(180,150)) +
#  ylim(c(-.55,-.12)) +
#  ylab("Temperature, C") +
#  xlab("Depth, m") +
#  coord_flip(expand=FALSE) +
#  ggtitle("Temperature Profile Cast 156")

#Download package that can make T-S Diagrams
#library(devtools)
#library(PlotSvalbard)

#ts_plot(dt=C156,temp_col= "T",sal_col="S")

#I found this package that makes T-S diagrams using ggplot
#I'm not entirely sure what's going on with this function. The graph lists potential temperature on the y axis, but I am not sure if it is just taking my in situ T values or if the conversion is built within the function

#Trying a different packege for T-S Diagram

#install_github("MikkoVihtakari/ggOceanMaps")
#install_github("MikkoVihtakari/ggOceanPlots")
#library(ggOceanPlots)
#library(ggOceanMaps)