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=" ")
#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)

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

colnames(dat1) = c("Cast", "Lon", "Lat", "Dep", "Sal", "Temp")

#Isolate a single cast (156)
#C156 = dat1[dat1$cast==156,]
#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()
for(i in 136:136)
{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])
 #SA<- as.numeric(data[,7])
 #PT<-as.numeric(data[,9])
 Lat<- as.numeric(data[,3])
 Lon<- as.numeric(data[,2])
 
 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) 
 
# 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
 
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

 }}      #output should be cast number, max depth, freshwater
## [1] 136
## [1] 337
## [1] 9
#create dataframe for binvalues, adjust their lengths so they all match

n=maxDepth-2
dat2<-data.frame(binDepth[1:n],binSalt[1:n],binTemp[1:n],Lon[1:n],Lat[1:n])

#dat2<-data.frame(binDepth,binSalt,binSA,binTemp,binPT)
colnames(dat2) = c("binD", "binS", "binT", "Lon", "Lat")
#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))

#ggplot(dat2,mapping=aes(x=binS,y=binPT))+
# geom_path() +
#  ylab("PT, C")+
#  xlab("Salinity")
#ggplot(dat2, mapping=aes(y=binD, x=binS))+
#  geom_path()+
#  ylim(300,0)
sst<-as.numeric(dat2[1,3])
up<- sst + 0.2
low<- sst - 0.2
print(sst)
## [1] 2.794444
#last=length(dat2)
  i=2
  while(dat2[i,3] <=up && dat2[i,3]>=low)  #&& i<=last )
  {
    print(dat2[i,3])
    
    SA2<- gsw_SA_from_SP(dat2[,2], dat2[,1], dat2[,4], dat2[,5])
    PT2<- gsw_pt_from_t(SA2, dat2[,3], dat2[,1], 0)
    
    i=i+1
    #print(mean(dat2[(2:i),5]))
    #print(sd(dat2[(2:i),5]))
    
    #print(mean(dat2[(2:i),3]))
    #print(sd(dat2[(2:i),3]))

  }
## [1] 2.784444
## [1] 2.716667
## [1] 2.688889
## [1] 2.68
## [1] 2.677778
## [1] 2.675556
## [1] 2.67
## [1] 2.668889
## [1] 2.662222
## [1] 2.66
## [1] 2.66
## [1] 2.66
## [1] 2.658889
## [1] 2.658889
## [1] 2.658889
## [1] 2.66
## [1] 2.66
## [1] 2.66
## [1] 2.662222
## [1] 2.612222
  print(i)
## [1] 22
z_SML<-binDepth[20]
print(z_SML)
## [1] 21
mean(PT2[1:z_SML])
## [1] 2.677033
sd(PT2[1:z_SML])
## [1] 0.04177571
mean(SA2[1:z_SML])
## [1] 30.0355
sd(SA2[1:z_SML])
## [1] 0.1015537
plot(PT2, binDepth, ylim = rev(range(binDepth)), xlab = "potential temperature, C", ylab = "Depth, m", col="black", lwd = 2, type="l")
points(PT2[z_SML], binDepth[z_SML], col = "red", lwd=2, pch = 20)

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

#Calculate Values using gsw package
#library(gsw)

# Pressure from depth
P = gsw_p_from_z(-Omg$dep,Omg$lat)
#absolute salinity from practical salinity
SA = gsw_SA_from_SP(Omg$sal,P,Omg$lon,Omg$lat)
#potential temperature from in situ t
PT = gsw_pt_from_t(SA,Omg$temp,P,0)
#conservative temp from pt
CT = gsw_CT_from_t(SA,Omg$temp,P)
#potential density at p=0
r0 = gsw_sigma0(SA,CT)
#potential density with p=1000
r1 = gsw_sigma1(SA,CT)
# Freezing Point
#This code gets an error
#tf = gsw_CT_freezing_poly(SA,p,1)
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)