library(readxl)
library(readr)
library(ggplot2)
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✓ tibble 3.1.2 ✓ dplyr 1.0.7
## ✓ tidyr 1.1.3 ✓ stringr 1.4.0
## ✓ purrr 0.3.4 ✓ forcats 0.5.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(janitor)##
## Attaching package: 'janitor'## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.testlibrary(lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(seacarb)## Loading required package: oce## Loading required package: gsw## Loading required package: sf## Linking to GEOS 3.8.1, GDAL 3.1.4, PROJ 6.3.1## Loading required package: testthat##
## Attaching package: 'testthat'## The following object is masked from 'package:dplyr':
##
## matches## The following object is masked from 'package:purrr':
##
## is_null## The following object is masked from 'package:tidyr':
##
## matcheslibrary(dplyr)
library(marelac)## Loading required package: shape##
## Attaching package: 'marelac'## The following objects are masked from 'package:oce':
##
## coriolis, gravitylibrary(gsubfn)## Loading required package: protolibrary(stringr)
### CO2 data: https://www.ncei.noaa.gov/access/ocean-carbon-data-system/oceans/Moorings/ndp097.html# Add your own working directory
setwd("~/Documents/Projects/CoastalOA/Code")
path <- "https://www.ncei.noaa.gov/access/ocean-carbon-data-system/oceans/Moorings/ndp097.html"
# Find names of files at NCEI portal
files <- readLines(path)[grep('.txt',readLines(path))]
files <- paste(str_extract_all(files,"(?<=href=\").+(?=\">timeseries)"))rm(list=setdiff(ls(), c("path","files")))
for(i in 1:length(files)){
#Read in Data
buoy <-read_tsv(file = files[i],comment="#", col_names = TRUE)
#Read in Metadata
metadata <- read_tsv(file = files[i],n_max=107)
#Parse metadata
site <-str_match(metadata[6,1],": (.*)")[2]
location <-str_match(metadata[7,1]," (.*)")[2]
lat <- str_match(location, ": (.*?),")[2]
long <- str_match(location, ", (.*)")[2]
buoy$lat <- lat
buoy$long <- long
buoy$site <- site
#Progress bar
print(paste("Working on: ",site," | Site ",i," out of ",length(files),sep=""))
pb <- txtProgressBar(min = 0, max = length(files), style = 3)
####### Speciate carbonate system + error #########
#Create dataset with complete cases of SST, SSS, pCO2_sw,pCO2_air, and pH
carb <- buoy[with(buoy,!is.na(buoy$SST)) & !is.na(buoy$SSS) & !is.na(buoy$pCO2_sw) & !is.na(buoy$pCO2_air) & !is.na(buoy$pH_sw),]
if(nrow(carb)>0){ #some sites have no pH data. Only run when full dataset is available
#Speciate carbonate system
t_carb <- carb(21,carb$pCO2_sw, carb$pH_sw, S=carb$SSS, T=carb$SST, Patm=1, P=0, Pt=0, Sit=0,k1k2="x", kf="x", ks="d", pHscale="T", b="u74", gas="insitu")
t_carb$datetime_utc <- carb$datetime_utc #add datetime
#Calculate buffer factors
buffers <- buffesm(21,carb$pCO2_sw, carb$pH_sw, S=carb$SSS, T=carb$SST,Patm = 1, P = 0, Pt = 0, Sit = 0,k1k2 = "x", kf = "x", ks = "d", pHscale = "T", b = "u74", warn = "n")
#combine carb system and buffers
t_carb <- bind_cols(t_carb, buffers)
#Calculate errors assuming 1uatm error pCO2 and 0.01 error pH
error <- errors(warn="n",21,evar1=1, evar2=0.01, carb$pCO2_sw, carb$pH_sw, S=carb$SSS, T=carb$SST, Patm=1, P=0, Pt=0, Sit=0,k1k2="x", kf="x", ks="d", pHscale="T", b="u74", gas="insitu")
names(error) <- paste(names(error),"_error",sep="")
#add to rest of carb data
t_carb <- bind_cols(t_carb,error)
#### convert pH to [H+] for hourly, daily, etc averaging
t_carb$H_plus <- 10^(-t_carb$pH)
# Remove redundant columns (pCO2, pH)
t_carb <- subset(t_carb, select=c(datetime_utc,H_plus,CO2,HCO3,CO3,DIC,ALK,OmegaAragonite, OmegaCalcite,gammaDIC,betaDIC,omegaDIC,gammaALK,betaALK,omegaALK,R,CO2_error,HCO3_error,CO3_error,DIC_error,ALK_error,OmegaAragonite_error,OmegaCalcite_error))
# add carb data to buoy data
buoy <- left_join(buoy, t_carb,by = "datetime_utc")
############### calculate pCO2_thermal and pCO2_non-thermal ######################
## using period of record mean SST and pCO2 because most years have data gaps, meaning calc would be biased by seasons with more/less data
# "pCO2 at Tobs" in Takahashi et al 2002
buoy$pCO2_therm <- mean(buoy$pCO2_sw,na.rm = TRUE) * exp(0.0423 * (buoy$SST - mean(buoy$SST,na.rm = TRUE)))
# "pCO2 at Tmean" in Takahashi et al 2002
buoy$pCO2_nontherm <- buoy$pCO2_sw * exp(0.0423 * (mean(buoy$SST,na.rm = TRUE) - buoy$SST))
}
# Merge site with previous sites
if (exists("buoy") == TRUE & i == 1) {
t <- buoy
} else {
if (exists("buoy") == TRUE & i > 1 & nrow(carb)>0)
t <- bind_rows(t, buoy)
}
setTxtProgressBar(pb, i)
rm(list=setdiff(ls(), c("path","files","i","t")))
} ## [1] "Working on: CCE1 | Site 1 out of 38"
##
|
| | 0%
|
|== | 3%[1] "Working on: Papa | Site 2 out of 38"
##
|
| | 0%
|
|==== | 5%[1] "Working on: KEO | Site 3 out of 38"
##
|
| | 0%
|
|====== | 8%[1] "Working on: JKEO | Site 4 out of 38"
##
|
| | 0%
|
|======= | 11%[1] "Working on: MOSEAN/WHOTS | Site 5 out of 38"
##
|
| | 0%
|
|========= | 13%[1] "Working on: TAO110W | Site 6 out of 38"
##
|
| | 0%
|
|=========== | 16%[1] "Working on: TAO125W | Site 7 out of 38"
##
|
| | 0%
|
|============= | 18%[1] "Working on: TAO140W | Site 8 out of 38"
##
|
| | 0%
|
|=============== | 21%[1] "Working on: TAO155W | Site 9 out of 38"
##
|
| | 0%
|
|================= | 24%[1] "Working on: TAO170W | Site 10 out of 38"
##
|
| | 0%
|
|================== | 26%[1] "Working on: TAO165E | Site 11 out of 38"
##
|
| | 0%
|
|==================== | 29%[1] "Working on: TAO8S165E | Site 12 out of 38"
##
|
| | 0%
|
|====================== | 32%[1] "Working on: Stratus | Site 13 out of 38"
##
|
| | 0%
|
|======================== | 34%[1] "Working on: BTM | Site 14 out of 38"
##
|
| | 0%
|
|========================== | 37%[1] "Working on: Iceland | Site 15 out of 38"
##
|
| | 0%
|
|============================ | 39%[1] "Working on: BOBOA | Site 16 out of 38"
##
|
| | 0%
|
|============================= | 42%[1] "Working on: SOFS | Site 17 out of 38"
##
|
| | 0%
|
|=============================== | 45%[1] "Working on: GAKOA | Site 18 out of 38"
##
|
| | 0%
|
|================================= | 47%[1] "Working on: Kodiak | Site 19 out of 38"
##
|
| | 0%
|
|=================================== | 50%[1] "Working on: SEAK | Site 20 out of 38"
##
|
| | 0%
|
|===================================== | 53%[1] "Working on: M2 | Site 21 out of 38"
##
|
| | 0%
|
|======================================= | 55%[1] "Working on: Cape Elizabeth | Site 22 out of 38"
##
|
| | 0%
|
|========================================= | 58%[1] "Working on: Chá bă | Site 23 out of 38"
##
|
| | 0%
|
|========================================== | 61%[1] "Working on: CCE2 | Site 24 out of 38"
##
|
| | 0%
|
|============================================ | 63%[1] "Working on: NH-10 | Site 25 out of 38"
##
|
| | 0%
|
|============================================== | 66%[1] "Working on: Ala Wai | Site 26 out of 38"
##
|
| | 0%
|
|================================================ | 68%[1] "Working on: Chuuk | Site 27 out of 38"
##
|
| | 0%
|
|================================================== | 71%[1] "Working on: CRIMP1 | Site 28 out of 38"
##
|
| | 0%
|
|==================================================== | 74%[1] "Working on: CRIMP2 | Site 29 out of 38"
##
|
| | 0%
|
|===================================================== | 76%[1] "Working on: Kaneohe | Site 30 out of 38"
##
|
| | 0%
|
|======================================================= | 79%[1] "Working on: Kilo Nalu | Site 31 out of 38"
##
|
| | 0%
|
|========================================================= | 82%[1] "Working on: Gray's Reef | Site 32 out of 38"
##
|
| | 0%
|
|=========================================================== | 84%[1] "Working on: Gulf of Maine | Site 33 out of 38"
##
|
| | 0%
|
|============================================================= | 87%[1] "Working on: Crescent Reef | Site 34 out of 38"
##
|
| | 0%
|
|=============================================================== | 89%[1] "Working on: Hog Reef | Site 35 out of 38"
##
|
| | 0%
|
|================================================================ | 92%[1] "Working on: Coastal MS | Site 36 out of 38"
##
|
| | 0%
|
|================================================================== | 95%[1] "Working on: Cheeca Rocks | Site 37 out of 38"
##
|
| | 0%
|
|==================================================================== | 97%[1] "Working on: La Parguera | Site 38 out of 38"
##
|
| | 0%
|
|======================================================================| 100%# add your own path to save data
write.csv(t,file=paste("~/Documents/Projects/CoastalOA/Code/Data/Summary","t.csv",sep="/"),row.names=FALSE)#Plotting
#ggplot(t)+ geom_point(aes(datetime_utc,pCO2_sw))+geom_point(aes(datetime_utc,DOXY),color="red")+facet_wrap(~site)