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:
#install.packages("dplyr")
#install.packages("tidyr")
#install.packages("ggplot2")
#install.packages("leaflet")
#install.packages("gsw")
#install.packages("matrixStats")
library(matrixStats)
library(tidyr)
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:matrixStats':
##
## count## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)
library(gsw)
library(leaflet)#cd /Users/w5/Mast467/Code
#./Hw3code
#ls -la output.dat
cp /Users/w5/Mast467/Data2021/output2021.dat output.datyear = 2020
omg <- read.csv(file="/Users/w5/Mast467/Code/output.dat",sep=" ")depth <- omg[,3]
Sal <- omg[,5]
temp <- omg[,4]
lon <- omg[,1]
lat <- omg[,2]
cast <- omg[,6]z=depth
s=Sal
t=temp
#for (i in 27:28)
#{
p=gsw_p_from_z(-omg[,3],omg[,2]) #pressure
#gsw must have depth as a negative to produce positive pressure
#print(omg[i,3])
#print(omg[i,2])
#print(p)
#} Debugging for loop
SA=gsw_SA_from_SP(s,p,omg[,1],omg[,2]) #Absolute Salinity
pt=gsw_pt_from_t(SA,t,p,0) #Potential Temperature
CT=gsw_CT_from_t(SA,t,p) #Conservative Temperature
r0=gsw_sigma0(SA,CT) #Potential Density to p(atm)
r1=gsw_sigma1(SA,CT) #Potential Density to p(1km)
tf=gsw_CT_freezing(SA,p,1) #Freezing Point
HC = abs(pt/tf) #Heat Content
#n2 <- gsw_Nsquared(SA,CT,p,lat)
#n2 <- as.data.frame(n2)
#n_2<- as.integer(tail(n2,1))
#n= n_2+1
omg1 <- data.frame (pt,r0,r1,s,SA,t,tf,z,CT,depth,temp,Sal,p,lon,lat,cast)
#data frame with generic and calculated variablesy = NULL
#Dmax = 0
#fresh = vector()
for (i in 0:500)
{ omgcast<-dplyr::filter(omg1, cast==i)
last=length(omgcast[,1])
#Named value last as number of rows to sift out empty casts
if (last >= 1)
{
depth<-as.numeric(omgcast[,10])
temp<-as.numeric(omgcast[,11])
Sal<-as.numeric(omgcast[,12])
p<- as.numeric(omgcast[,13])
SA<- as.numeric(omgcast[,5])
pt<- as.numeric(omgcast[,1])
CT<- as.numeric(omgcast[,9])
r0<- as.numeric(omgcast[,2])
r1<- as.numeric(omgcast[,3])
tf<- as.numeric(omgcast[,7])
lat<- as.numeric(omgcast[,15])
lon<- as.numeric(omgcast[,14])
#extract measurements/ calculations into vector in data frame omgcast
maxdepth<-as.integer(tail(depth,1))
#print(maxdepth)
bindepth<-seq(2,maxdepth,1)
binSal<-binMeans(Sal, 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)
binr0<-binMeans(r0, x=depth, bx=bindepth, na.rm=TRUE)
binr1<-binMeans(r1, x=depth, bx=bindepth, na.rm=TRUE)
bintf<-binMeans(tf, x=depth, bx=bindepth, na.rm=TRUE)
binlat <- binMeans(lat, x=depth, bx=bindepth, na.rm=TRUE)
binlon <- binMeans(lon, x=depth, bx=bindepth, na.rm=TRUE)
bindepth = head(bindepth,-1)
bindata = data.frame (bindepth,binSal,bintemp,binSA,binpt,binp,binCT,binr0,binr1,bintf,binlat,binlon)
#bindata <- na.omit(bindata)
# Fresh Water Content
Sr = 34.8
Sr <- (rep(Sr,length(bindepth)))
fresh = round(sum(((Sr-binSal)/Sr), na.rm = TRUE), digits = 4)
# Mixed Layer Depth
#SST = bindata[4,3]
SST <- as.numeric(bindata[4,3])
Tsum = 0
Tavg = 0
Ssum = 0
n = 0
u = 4
while(bindata[u,3] <= SST+.5 && bindata[u,3] >= SST-.5)
{
n = n+1
Tsum = bindata[u,3]+Tsum
Ssum = bindata[u,2]+Ssum
if(is.na(bindata[u+1,3]))
{bindata[u+1,3]=FALSE}
u = u+1
}
Zsml <- u-1
Tavg = Tsum/n
Savg = Ssum/n #finding the temperature and salinity mean above the mixed layer depth
#print(Tavg)
#print(Zsml)
#print(Savg)
Tsqsum = 0
Ssqsum = 0
for(l in 2:Zsml)
{ dT = bindata[l,3]-Tavg
dS = bindata[l,2]-Savg
Tsqsum = Tsqsum + dT*dT
Ssqsum = Ssqsum + dS*dS
}
Ssdev = sqrt(Ssqsum/n)
Tsdev = sqrt(Tsqsum/n)
# attempting to put values for each cast into data frame
#Castvalues = data.frame(fresh, Ssdev, Tsdev)
bindepth_vec<-as.numeric(bindata[,1])
binSal_vec<-as.numeric(bindata[,2])
bintemp_vec<-as.numeric(bindata[,3])
binSA_vec<-as.numeric(bindata[,4])
binpt_vec<-as.numeric(bindata[,5])
binp_vec<-as.numeric(bindata[,6])
binCT_vec<-as.numeric(bindata[,7])
binr0_vec<-as.numeric(bindata[,8])
binr1_vec<-as.numeric(bindata[,9])
bintf_vec<-as.numeric(bindata[,10])
binlat_vec<-as.numeric(bindata[,11])
binlon_vec<-as.numeric(bindata[,12])
#convert bin values to vectors and put them in a data frame
binvectors <- data.frame(bindepth_vec,binSal_vec,bintemp_vec,binSA_vec,binpt_vec,binp_vec,binCT_vec,binr0_vec,binr1_vec,bintf_vec,binlat_vec,binlon_vec)
#binvectors <- na.omit(binvectors)
#finding tmax and the depth at which tmax occurs
tmax=-99
for (ii in 50:(tail(bindepth_vec, n=1)-20))
{
if(is.na(binvectors[ii,3]))
{binvectors[ii,3] = FALSE}
else{ if (binvectors[ii,3] >= tmax)
{ tmax = binvectors[ii,3]
depth_tmax=ii+2
Sal_tmax = binvectors[ii,2]
lat_tmax = binvectors[ii,11]
lon_tmax = binvectors[ii,12]
Tmax_values <- data.frame(tmax,Sal_tmax,lat_tmax,lon_tmax,depth_tmax,ii)
}
}
}
y<- rbind.data.frame(y, Tmax_values)
}
}
#write.csv(y,file = "/Users/w5/Mast467/Data/TmaxData")
#median(y[,1], na.rm = FALSE)#y16
#y17
#y18
#y19
#y20
#y21
#write.table(y,file = "/Users/w5/Mast467/Data/TmaxData")
#Tmaxdat2 <- read.csv(file="/Users/w5/Mast467/Data/TmaxData",sep=",")
#Tmaxmedian2021 = median(y[,1], na.rm = FALSE)
#alltmax<- rbind(y16,y17,y19,y20,y21)
#Tmaxmedian2016 = median(y16[,1])
#Tmaxmedian2017 = median(y17[,1])
#Tmaxmedian2018 = median(y18[,1])
#Tmaxmedian2019 = median(y19[,1])
#Tmaxmedian2020 = median(y20[,1])
#Tmaxmedian2021 = median(y21[,1])#Tmaxplotdata <- data.frame(year = c(2016,2017,2018,2019,2020,2021),Tmaxmedians = c(Tmaxmedian2016,Tmaxmedian2017,Tmaxmedian2018,Tmaxmedian2019,Tmaxmedian2020,Tmaxmedian2021))#```{r} ggplot(Tmaxplotdata,aes(x=year,y=Tmaxmedians)) + geom_path(color = âredâ) + geom_point(size = 4)+ ylab(âTmax,Câ) + xlab(âYearâ)
```r
Tmaxmedian2016 = 1.818
Tmaxmedian2017 = 1.52
Tmaxmedian2018 = 1.56
Tmaxmedian2019 = 1.46
Tmaxmedian2020 = 1.197
Tmaxmedian2021 = 1.207#```{r} ggplot(bindata,aes(x=bintemp,y=bindepth)) + geom_path() + geom_point(aes(x= bindata[Zsml-1,3], y=Zsml), colour = âredâ) + ylab(âAveraged Depth, mâ) + xlab(âTemperature, Câ) + ylim(1.5*(Zsml),0) + xlim(-1,.5)
#```{r}
ggplot(bindata,aes(x=bintemp,y=bindepth)) +
geom_path() +
geom_point(aes(x= bindata[Zsml-1,3], y=Zsml), colour = "red") +
ylab("Averaged Depth, m") +
xlab("Temperature, C") +
ylim(1.5*(Zsml),0) +
xlim(-1,.5)
#{r} ggplot(Bindata,aes(x=binSal,y=bindepth)) + geom_path() + ylab("Averaged Depth, m") + xlab("Averaged Salinity, PPT") #```{r} ggplot(Bindata,aes(x=bintemp,y=bindepth)) + geom_path() + ylab(âAveraged Depth, mâ) + xlab(âTemperature, Câ) + ylim(150,0)
#```{r}
ggplot(omg,aes(x=depth,y=temp)) +
geom_path() +
ylab("Temperature, C") +
xlab("Depth, m") +
xlim(750,0) +
coord_flip(expand=FALSE)#{r} ggplot(omg,aes(x=depth,y=Sal)) + geom_path() + ylab("Salinity, PPT") + xlab("Depth, m") + coord_flip(expand=FALSE) #{r} ggplot(omg,aes(x=Sal,y=temp)) + geom_path() + ylab("Temperature, C") + xlab("Salinity, PPT") + coord_flip(expand=FALSE) #{r} ggplot(omg,aes(x=Sal,y=temp,color=depth)) + geom_path() + ylab("Temperature, C") + xlab("Salinity, PPT") + scale_color_viridis_c(guide=guide_colorbar(barheight=16.5)) + coord_flip(expand=FALSE) #{r} ggplot(omg,aes(x=depth,y=temp,color=Sal)) + geom_path() + ylab("Temperature, C") + xlab("Depth, m") + scale_color_viridis_c(guide=guide_colorbar(barheight=16.5)) + coord_flip(expand=FALSE) #{r} ggplot(omg,aes(x=depth,y=Sal,color=temp)) + geom_path() + ylab("Salinity, PPT") + xlab("Depth, m") + scale_color_viridis_c(guide=guide_colorbar(barheight=16.5)) + coord_flip(expand=FALSE) #{r} ggplot(omg,aes(x=depth,y=Sal,color=temp)) + geom_path() + ylab("Salinity, PPT") + xlab("Depth, m") + ylim(32,35) + xlim(750,0) + scale_color_viridis_c(guide=guide_colorbar(barheight=16.5)) + coord_flip(expand=FALSE)
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.