R/correct.merge.IOP.profile.R
In belasi01/Riops: IOPs data processing
Defines functions
correct.merge.IOP.profile
Documented in
correct.merge.IOP.profile
#' Apply corrections and merge IOP data from an optical package
#'
#' Apply corrections and merge data coming from various instruments for a
#' given vertical profile in the water column.
#'
#' @param instrument is a list or a data frame of instruments available. For example,
#' instrument = list(ASPH=1, HS6=1, CTD.UQAR=1)
#' @param parameters is a list or a data frame with the processing
#' parameters provided by the user. When called from
#' \code{\link{process.IOPs}}, the parameters are red in cast.info.dat
#' and cal.info.dat. A detailed description of each parameter is provided
#' below (Details section)
#'
#' @return
#' It retruns a large list with all the raw and corrected data of each instrument
#' The list is also save as IOP.RData format. One or tow additional RData are save
#' to store the interpolated data for the down and up cast respectively
#' (IOP.fitted.down.RData and IOP.fitted.up.RData). The upcast
#' is optional and depends on the user input (parameter maxx).
#'
#' @details
#' This is the heart of the \pkg{Riops} package. All the data are red,
#' matched using a unique time frame, corrected and interpolated on a common
#' grid of equally spaced depth. The detailed description of the correction
#' is not described here (To be done in a future version). Here is a summary
#' of the minimum require to run the code succesfully.
#'
#' The processing parameters are stored in two ascii files named cast.info.dat
#' and cal.info.dat. The former file includes the following parameters
#' \itemize{
#' \item{\strong{lon}} { is the longitude of the station. It is not used in the processing
#' it will be in the PDF reporting (see \code{\link{create.report}})}
#' \item{\strong{lat}} { in the latitude of the station. It is not used in the processing
#' it will be in the PDF reporting (see \code{\link{create.report}})}
#' \item{\strong{cast}} { is a character string corresponding to either MINIDAS cast number
#' or file extension generated by the WAP from DH4 archive file (e.g. "001")}
#' \item{\strong{Time0.CTD}}{ is a character string for CTD start time. It will be transformed
#' into a time stamp using \code{\link{as.POSIXct}} and should respect this format "2015-05-06 12:18:01".
#' However, if NA the program it check in other instruments (e.g. ASPH or BB3) to find the starting
#' time. If no absolute time is found, the user will be asked to add a time stamp manually in the
#' cast.info.dat file)}
#' \item{\strong{Time0.LISST}}{ is a character string for LISST starting time. It will be transformed
#' into a time stamp using \code{\link{as.POSIXct}} and should respect this format "2015-05-06 12:18:01".
#' However, if NA the program it check in other instruments (e.g. ASPH or BB3) to find the starting
#' time. If no absolute time is found, the user will be asked to add a time stamp manually in the
#' cast.info.dat file)}
#' \item{\strong{minx}}{ is the index of the CTD timer corresponding to the start of the profile.
#' If NA, the user will be prompted to select click on the start of the profile on a plot
#' showing the depth versus time. This a call to the function \code{\link{identify}}. Once
#' this is done, the user can adjust the begining of the cast by changing the values in
#' the cast.info.dat file afterward (won't be prompted if minx!=NA).}
#' \item{\strong{maxx}}{ is the index of the CTD corresponding to the End of the profile.
#' Similar to \strong{minx}.}
#' \item{\strong{Zint}}{ is the depth interval for profile smooting
#' (e.g. 1 will result in a fitted profile every 1 meter)}
#' \item{\strong{depth.interval.for.smoothing}}{ is the depth interval, in meter, that will
#' be used to smooth data with the loess function. It will be converted into a "span" value.
#' (see \code{\link{loess}})}
#' \item{\strong{asph.skip}}{ is the number or records to skip at the begining of the a-sphere file.
#' This is necessary only when the a-sphere file contains 2 or more profiles. This
#' happens some time when the MiniDAS reset the cast number. Otherwise 0 will fit
#' most situations.}
#' \item{\strong{maxbb}}{ is a parameter for the backscattering plots generated in the PDF report.}
#' \item{\strong{Ndepth.to.plot}}{ is the number of spectra along the depth profile to
#' plot in the report. 8 to 12 usually is enough}
#'}
#'
#' The calibration information are stored in the cal.info.dat.
#' The file includes the following information
#' (NOTE: any field can be ommited with out problem):
#'
#' \itemize{
#' \item{\strong{Tref.ASPH}} { is the temperature of pure water used by
#' Hobilabs for the ASPH calibration.
#' In 2010 and 2014 it was 13.2,
#' in 2013 it was 19,
#' in 2016 it was 14.4.}
#' \item{\strong{HS6.CALYEAR}}{ is the year of the HS6 calibration.}
#' \item{\strong{Tref.ACS}}{ is the temperature of pure water used by
#' Wetlabs for the ACS calibration. The IML ACs calibration of 2013 was 20.3}
#' \item{\strong{scat.correction}}{ is the method for the scattering correction of AC-S abssorption
#' (eg : "mckee","zaneveld", "baseline","none")}
#' \item{\strong{blank.ASPH}}{ is a string for the path of the blank for ASPH as created by
#' \code{\link{analyse.ASPH.blank}}.}
#' \item{\strong{blank.ACS}}{ is a string for the path of the blank for ACS as created by
#' \code{\link{analyse.ACs.blank}}.}
#' \item{\strong{blank.BB9}}{ is a string for the path of the blank for BB9.}
#' \item{\strong{blank.BB3}}{ is a string for the path of the blank for BB3.}
#' }
#'
#'
#' @author Simon Belanger
#'
#' @export
#'
#'
correct.merge.IOP.profile <- function(instrument, parameters){
####### Get prepared to process data #####
print("Get prepared to process data")
path = parameters$path
cast = parameters$cast
Time0.CTD = parameters$Time0.CTD
Time0.LISST = parameters$Time0.LISST
minx= parameters$minx
maxx= parameters$maxx
Zint= parameters$Zint
asph.skip = parameters$asph.skip
depth.span= parameters$depth.interval.for.smoothing
# Check which instrument is missing in the list
if (is.null(instrument$CTD.UQAR)) instrument$CTD.UQAR <- 0
if (is.null(instrument$HS6)) instrument$HS6 <- 0
if (is.null(instrument$FLECO)) instrument$FLECO <- 0
if (is.null(instrument$ASPH)) instrument$ASPH <- 0
if (is.null(instrument$CTD.DH4)) instrument$CTD.DH4 <- 0
if (is.null(instrument$ACS)) instrument$ACS <- 0
if (is.null(instrument$BB9)) instrument$BB9 <- 0
if (is.null(instrument$BB3)) instrument$BB3 <- 0
if (is.null(instrument$LISST)) instrument$LISST <- 0
if (is.null(instrument$FLBBCD)) instrument$FLBBCD <- 0
if (is.null(instrument$FLCHL)) instrument$FLCHL <- 0
# Check information from cal.info
if (!is.null(parameters$Tref.ASPH)) {
Tref.ASPH = parameters$Tref.ASPH
} else {
if (instrument$ASPH == 1) {
print("WARNING: A-sphere calibration temperature of PW must be provided!!!")
print("Add Tref.ASPH in the cal.info.dat file")
print("Stop processing")
retrun(0)
}
Tref.ASPH = NA
}
if (!is.null(parameters$Tref.ACS) & !is.null(parameters$scat.correction)) {
Tref.ACS = parameters$Tref.ACS
scat.correction = parameters$scat.correction
} else {
if (instrument$ACS == 1) {
print("WARNING: AC-s calibration temperature of PW and scat.correction method must be provided!!!")
print("Add Tref.ACS and scat.correction in the cal.info.dat file")
print("Stop processing")
retrun(0)
}
Tref.ACS = NA
scat.correction = NA
}
if (!is.null(parameters$HS6.CALYEAR) & !is.na(parameters$HS6.CALYEAR)) {
HS6.CALYEAR = parameters$HS6.CALYEAR
print(paste("Hydroscat-6 calibration year is", HS6.CALYEAR))
} else {
print(paste("No Hydroscat-6 calibration year provided"))
if (instrument$HS6 == 1) {
print("WARNING: HS-6 calibration year must be provided!!!")
print("Add HS6.CALYEAR in the cal.info.dat file")
print("Stop processing")
return(0)
}
HS6.CALYEAR = NA
}
if (!is.null(parameters$blank.BB9) & !is.na(parameters$blank.BB9)) {
index_ext = length(unlist(strsplit(parameters$blank.BB9, "[.]")))
ext = unlist(strsplit(parameters$blank.BB9, "[.]"))[index_ext]
if (ext == "dev") BB9.DEVICE.FILE = TRUE else BB9.DEVICE.FILE = FALSE
}
if (!is.null(parameters$blank.BB3) & !is.na(parameters$blank.BB3)) {
index_ext = length(unlist(strsplit(parameters$blank.BB3, "[.]")))
ext = unlist(strsplit(parameters$blank.BB3, "[.]"))[index_ext]
if (ext == "dev") BB3.DEVICE.FILE = TRUE else BB3.DEVICE.FILE = FALSE
}
setwd(path)
####### Read all data#########
print("Read all data")
# This was added in order to process previous instrument file
if (!is.null(instrument$CTD.IML)) instrument$CTD.DH4 = instrument$CTD.IML
if (instrument$CTD.UQAR == 1) {
filen = paste("CTD",cast, ".TXT", sep="")
# If the standard CTD file do not exist,
# check for a cnv file instead
if (file.exists(filen)) {
CTD = read.CTD(filen)
} else {
filen = list.files( pattern = "*\\.cnv")
if (file.exists(filen)) {
print("Reading CNV file using oce read function")
CTD = list()
CTD.OCE = read.ctd.sbe(filen)
# Extract data from the CTD object and store as list
if (is.na(Time0.CTD)) Time0.CTD = as.POSIXct(format(CTD.OCE[["startTime"]], tz="GMT"), tz="GMT")
if (is.na(Time0.CTD)) {
print("No CTD time available")
print("Please add a time string for Time0.CTD")
print("in cast.info.dat")
print("STOP PROCESSING")
return(0)
}
if (is.numeric(CTD.OCE[["time"]])) { # This condition was added because
#the "time" in the cnv files are now output
#in POSIXct by the oce package since 2018.
CTD$Time = Time0.CTD + CTD.OCE[["time"]]
} else CTD$Time = CTD.OCE[["time"]]
CTD$Temp = CTD.OCE[["temperature"]]
CTD$Sal = CTD.OCE[["salinity"]]
if (!is.null(CTD.OCE[["par"]])) CTD$PAR<- CTD.OCE[["par"]]
if (!is.null(CTD.OCE[["turbidity"]])) CTD$turbidity <- CTD.OCE[["turbidity"]]
if (!is.null(CTD.OCE[["fluorescence"]])) CTD$fluorescence <- CTD.OCE[["fluorescence"]]
if (!is.null(CTD.OCE[["depth"]])) CTD$Depth= swDepth(CTD.OCE[["depth"]]) else CTD$Depth= swDepth(CTD.OCE[["pressure"]])
if (!is.null(CTD.OCE[["sigmaTheta"]])) CTD$sigmaTheta=CTD.OCE[["sigmaTheta"]]
}
}
} else {
print("No CTD.UQAR to process")
if (instrument$CTD.DH4 == 1) {
print("Read Port Assossiation for DH4")
default.DH4.ports.file <- paste( Sys.getenv("R_IOPs_DATA_DIR"), "DH4.ports.dat", sep = "/")
DH4.ports.file <- paste(parameters$path, "DH4.ports.dat", sep = "/")
if(!file.exists(DH4.ports.file)) {
file.copy(from = default.DH4.ports.file, to = DH4.ports.file)
cat("EDIT file", DH4.ports.file, "and CUSTOMIZE IT\n")
stop()
}
DH4.ports <- read.table(DH4.ports.file, sep=";", header = T,
colClasses = c("character", "character"))
ix.inst = which(DH4.ports$Sensor == "CTD")
pattern = paste("_", DH4.ports$Port[ix.inst],"_", sep="")
filen = list.files( pattern = pattern)
if(!file.exists(filen)) {
print(paste("CTD file not found in", getwd()))
print("Check whether or not the port number is good")
print("Use capital letters for sensor name")
print(DH4.ports)
stop()
} else CTD = read.CTD.DH4(filen)
} else {
print("No CTD to process")
# print("Abort processing")
CTD = list()
# return(1)
}
}
if (instrument$ASPH == 1) {
filen = list.files( pattern = "CST.*\\.txt")
if (length(filen) == 0) filen = list.files( pattern = "ASPH.*\\.txt")
if (length(filen) > 1) {
ix.filen = which.min(str_length(filen))
filen = filen[ix.filen]
}
ASPH = read.ASPH(filen, skip=asph.skip) # modifiy by Simon Belanger
} else {
print("No ASPH to process")
ASPH = list()
# This condition is added because the CTD from IML has no absolute time information
# So the absolute time is given by the ASPH instrument.
# if (instrument$CTD.DH4 == 1) {
# print("Abort processing")
# return(0)
# }
}
if (instrument$HS6 == 1) {
filen = paste("./HS6",cast,".dat", sep="")
if (!file.exists(filen)) filen = list.files( pattern = ".*\\HS.*\\.dat")
HS6 = read.HS6(filen)
} else {
print("No Hydroscat-6 to process")
HS6 = list()
}
if (instrument$FLECO == 1) {
FLECO = read.FLECO(paste("./FL",cast,".TXT", sep=""))
} else {
print("No fluorescence ECOtriplet to process")
FLECO = list()
}
if (instrument$BB9 == 1) {
# BB9 attached to a DH4...
if (instrument$CTD.DH4 == 1) {
ix.inst = which(DH4.ports$Sensor == "BB9")
pattern = paste("_", DH4.ports$Port[ix.inst],"_", sep="")
filen = list.files( pattern = pattern)
if(!file.exists(filen)) {
print(paste("BB9 file not found in", getwd()))
print("Check whether or not the port number is good")
print("Use capital letters for sensor name")
print(DH4.ports)
stop()
} else BB9 = read.BB9(filen)
# extract the port number
BB9.port = as.numeric(DH4.ports$Port[ix.inst])
} else {
# BB9 in stand alone mode.... data could be provided in calibrated format,
# or in raw format which the prefered method.
pattern = "BB9"
filen = list.files( pattern = pattern)
if(!file.exists(filen)) {
print(paste("BB9 file not found in", getwd()))
print("BB9 should appear in the file name...")
stop()
} else {
# Check the extension. If raw, then apply calibration from the device file.
index_ext = length(unlist(strsplit(filen, "[.]"))) # for station names with periods, ex. G604.5
ext = unlist(strsplit(filen, "[.]"))[index_ext]
if (ext == "raw") {
print("BB9 deployed in Stand Alone mode and raw format")
BB9.raw <- read.BB9(filen, raw = TRUE)
if (BB9.DEVICE.FILE) {
BB9 <- apply.ECO.cal(BB9.raw, dev.file=parameters$blank.BB9, dark.file=NA,ECO.type="BB9")
} else {
print("Raw BB9 detected but no device file provided!!! ")
print("Please put the path and file name of the device file")
print("in the cal.info.dat for the >blank.BB9< parameter")
BB9.DEVICE.FILE == FALSE
}
} else {
print("BB9 deployed in Stand Alone mode and eng format")
BB9 <- read.BB9(filen)
}
}
}
} else {
print("No BB9 to process")
BB9 = list()
}
if (instrument$ACS == 1) {
# ACS attached to a DH4...
if (instrument$CTD.DH4 == 1) {
ix.inst = which(DH4.ports$Sensor == "ACS")
pattern = paste("_", DH4.ports$Port[ix.inst],"_", sep="")
filen = list.files( pattern = pattern)
if(!file.exists(filen)) {
print(paste("ACS file not found in", getwd()))
print("Check whether or not the port number is good")
print("Use capital letters for sensor name")
print(DH4.ports)
stop()
} else ACS = read.ACs(filen)
# extract the port number from name
ACS.port = as.numeric(DH4.ports$Port[ix.inst])
} else {
# ACS in stand alone mode.... data could be provided in calibrated format,
# or in raw format which the prefered method.
pattern = "ACS"
filen = list.files( pattern = pattern)
if(!file.exists(filen)) {
print(paste("ACS file not found in", getwd()))
print("ACS should appear in the file name...")
stop()
} else {
ACS <- read.ACs(filen)
}
}
} else {
print("No ACS to process")
ACS = list()
}
if (instrument$LISST == 1){
ix.inst = which(DH4.ports$Sensor == "LISST")
pattern = paste("_", DH4.ports$Port[ix.inst],"_", sep="")
filen = list.files( pattern = pattern)
if(!file.exists(filen)) {
print("Check for LISST in stand alone mode")
filen = list.files( pattern = "L.*\\.asc")
if (!file.exists(filen)) {
print(paste("LISST file not found in", getwd()))
print("Check whether or not the port number is good")
print("Use capital letters for sensor name")
print(DH4.ports)
stop()
} else LISST = read.LISST(filen)
} else LISST = read.LISST(filen)
} else {
print("No LISST to process")
LISST = list()
}
if (instrument$BB3 == 1){
ix.inst = which(DH4.ports$Sensor == "BB3")
pattern = paste("_", DH4.ports$Port[ix.inst],"_", sep="")
filen = list.files( pattern = pattern)
if(!file.exists(filen)) {
print(paste("BB3 file not found in", getwd()))
print("Check whether or not the port number is good")
print("Use capital letters for sensor name")
print(DH4.ports)
stop()
} else BB3 = read.BB3(filen)
BB3.port = as.numeric(DH4.ports$Port[ix.inst])
} else {
print("No BB3 to process")
BB3 = list()
}
if (instrument$FLBBCD == 1){
ix.inst = which(DH4.ports$Sensor == "FLBBCD")
pattern = paste("_", DH4.ports$Port[ix.inst],"_", sep="")
filen = list.files( pattern = pattern)
if(!file.exists(filen)) {
print(paste("FLBBCD file not found in", getwd()))
print("Check whether or not the port number is good")
print("Use capital letters for sensor name")
print(DH4.ports)
stop()
} else FLBBCD = read.FLBBCD(filen)
FLBBCD.port = as.numeric(DH4.ports$Port[ix.inst])
} else {
print("No FLBBCD to process")
FLBBCD = list()
}
if (instrument$FLCHL == 1){
ix.inst = which(DH4.ports$Sensor == "FLCHL")
pattern = paste("_", DH4.ports$Port[ix.inst],"_", sep="")
filen = list.files( pattern = pattern)
if(!file.exists(filen)) {
print(paste("FLCHL file not found in", getwd()))
print("Check whether or not the port number is good")
print("Use capital letters for sensor name")
print(DH4.ports)
stop()
} else FLCHL = read.FLCHL(filen)
FLCHL.port = as.numeric(DH4.ports$Port[ix.inst])
} else {
print("No FLCHL to process")
FLCHL = list()
}
####### Adjusting time stamp to each sensor; two cases are considered####
######## for 1) IML optical package and 2) UQAR optical package, or SBE CTD
if (instrument$CTD.DH4 == 1) {
# Reads Time OFFSET for each sensor
filen = list.files( pattern = "_TO.")
TO = read.table(filen, header=T)
if (is.na(Time0.CTD) & instrument$ASPH == 1) {
#Time0.CTD = find.Time0.CTD(path,cast)
# Find time0 from ASPH to apply time stamp to CTD, BB9 and ACS
print("Create Time Stamp from ASPHERE")
ix.maxZ.CTD = which.max(CTD$Depth)
ix.maxZ.ASPH = which.max(ASPH$depth)
Time0.CTD = ASPH$time[ix.maxZ.ASPH]- CTD$Timer[ix.maxZ.CTD]/1000
# Plot
plot(CTD$Time, CTD$Depth, pch=19,cex=0.5, main=" depth matching (adjust cast.info if needed)")
points(ASPH$time, ASPH$depth, col=2, pch=19, cex=0.2)
legend("topright", c("CTD", "ASPH"), col=c(1,2), pch=c(19,19))
}
if (is.na(Time0.CTD) & instrument$BB3 == 1) {
#Time0.CTD = find.Time0.CTD(path,cast)
# Find time0 from ASPH to apply time stamp to CTD, BB9 and ACS
print("Create Time Stamp from BB3")
Time0.CTD = BB3$Time[1]
}
if (is.na(Time0.CTD)) {
print("No instrument have absolute time string")
print("You need to provide a POSXCT time string")
print("for Time0.CTD and Time0.LISST in the cast.info.dat file")
print("(e.g.2015-05-06 12:18:01 ")
stop()
}
CTD$Time = Time0.CTD + CTD$Timer/1000
# NOTE the Timer, when using the DH4, are in milliseconds. This is why we divide by 1000
if (instrument$ACS == 1) ACS$Time = Time0.CTD + ACS$Timer/1000 - TO$offset[ACS.port]/1000
if (instrument$BB9 == 1) BB9$Time = Time0.CTD + BB9$Timer/1000 - TO$offset[BB9.port]/1000
if (instrument$BB3 == 1) BB3$Time = Time0.CTD + BB3$Timer/1000 - TO$offset[BB3.port]/1000
if (instrument$FLBBCD == 1) FLBBCD$Time = Time0.CTD + FLBBCD$Timer/1000 - TO$offset[FLBBCD.port]/1000
if (instrument$FLCHL == 1) FLCHL$Time = Time0.CTD + FLCHL$Timer/1000 - TO$offset[FLCHL.port]/1000
dt = difftime(max(CTD$Time), min(CTD$Time), units="secs")
# use a 10 seconds time window to define the span parameter for smoothing data
span.CTD = max(10 / as.numeric(dt), 0.1)
# Find time0 from ASPH to apply time stamp to LISST
if (instrument$LISST == 1) {
if (is.na(Time0.LISST) & instrument$ASPH == 1) {
# Check if LISST is in stand alone mode or attached to DH4
ix.inst = which(DH4.ports$Sensor == "LISST")
if(length(ix.inst) == 0) {
ix.maxZ.LISST = which.max(LISST$Depth)
ix.maxZ.ASPH = which.max(ASPH$depth)
# the 1.02932 is the asquisition rate in seconds of the LISST in stand alone mode.
Time0.LISST = ASPH$time[ix.maxZ.ASPH] - ix.maxZ.LISST*1.029432
LISST$Time = Time0.LISST + (1:length(LISST$time))*1.029432
} else {
ix.maxZ.LISST = which.max(LISST$Depth)
ix.maxZ.ASPH = which.max(ASPH$depth)
Time0.LISST = ASPH$time[ix.maxZ.ASPH] - ix.maxZ.LISST
LISST$Time = Time0.LISST + (1:length(LISST$time))
}
}
ix.inst = which(DH4.ports$Sensor == "LISST")
if(length(ix.inst) == 0) {
# the LISST time step in stand alone mode is 1.0294, while it is equal to 1 when attached to DH4
LISST$Time = Time0.LISST + (1:length(LISST$time))*1.029432
} else LISST$Time = Time0.LISST + (1:length(LISST$time))
plot(LISST$Time, LISST$Depth, pch=19,cex=0.5, main=" depth matching (adjust cast.info if needed)")
points(ASPH$time, ASPH$depth, col=2, pch=19, cex=0.2)
legend("topright", c("LISST", "ASPH"), col=c(1,2), pch=c(19,19))
}
}
if (instrument$CTD.UQAR == 1) {
# Assume CTD has the correct time
if (!is.na(Time0.CTD)) {
diffCTD = difftime(CTD$Time[1], Time0.CTD, units="secs")
CTD$Time = CTD$Time + diffCTD
} else {
Time0.CTD = CTD$Time[1]
diffCTD = difftime(CTD$Time[1], Time0.CTD, units="secs") # so difftime ==0
}
print(diffCTD)
CTD$Time <- CTD$Time + diffCTD
#### modification for processing time series
if (parameters$timeserie == 1) {
print("Correct time difference based on the time when each instruments hit the water")
plot(CTD$Time, CTD$Depth)
print("Click for the CTD time stamp")
ixCTD = identify(CTD$Time, CTD$Depth)
} else {
ix.maxZ.CTD = which.max(CTD$Depth)
}
if (instrument$ASPH == 1) {
ix.maxZ.ASPH = which.max(ASPH$depth)
if (parameters$timeserie == 1) {
plot(ASPH$time, ASPH$depth)
print("Click for the ASPH time stamp")
ixASPH = identify(ASPH$time, ASPH$depth)
diffASPH = difftime(CTD$Time[ixCTD], ASPH$time[ixASPH], units="secs")
} else {
diffASPH = difftime(CTD$Time[ix.maxZ.CTD], ASPH$time[ix.maxZ.ASPH], units="secs")
}
ASPH$time = ASPH$time + diffASPH
plot(ASPH$time, ASPH$depth, pch=19,cex=0.5, main="Check depth matching",
sub = "Adjust Time0.CTD in cast.info if needed")
points(CTD$Time, CTD$Depth, col=2, pch=19, cex=0.2)
legend("topright", c("ASPH", "CTD"), col=c(1,2), pch=c(19,19))
text(ASPH$time[ix.maxZ.ASPH], 0.8*ASPH$depth[ix.maxZ.ASPH],
paste('Time diff:', as.character(signif(diffASPH,3))))
png("CTDvsASPH_depth.png", height = 5, width = 7, units = "in", res=300)
plot(ASPH$time, ASPH$depth, pch=19,cex=0.5, main="Check depth matching",
sub = "Adjust Time0.CTD in cast.info if needed")
points(CTD$Time, CTD$Depth, col=2, pch=19, cex=0.2)
legend("topright", c("ASPH", "CTD"), col=c(1,2), pch=c(19,19))
text(ASPH$time[ix.maxZ.ASPH], 0.8*ASPH$depth[ix.maxZ.ASPH],
paste('Time diff:', as.character(signif(diffASPH,3))))
dev.off()
}
if (instrument$HS6 == 1) {
ix.maxZ.HS6 = which.max(HS6$depth)
if (parameters$timeserie == 1) {
print("Correct time difference based on the time when each instruments hit the water")
plot(HS6$Time, HS6$depth)
print("Click for the HS6 time stamp")
ixHS6 = identify(HS6$Time, HS6$depth)
plot(CTD$Time, CTD$Depth)
print("Click for the CTD time stamp")
ixCTD = identify(CTD$Time, CTD$Depth)
diffHS6 = difftime(CTD$Time[ixCTD], HS6$Time[ixHS6], units="secs")
HS6$Time = HS6$Time + diffHS6
} else {
diff.z = abs(HS6$depth[ix.maxZ.HS6] - CTD$Depth[ix.maxZ.CTD])
if (diff.z < 1) {
diffHS6 = difftime(CTD$Time[ix.maxZ.CTD], HS6$Time[ix.maxZ.HS6], units="secs")
HS6$Time = HS6$Time + diffHS6
} else {
print("Correct time difference based on the time when each instruments hit the water")
plot(HS6$Time, HS6$depth)
print("Click for the HS6 time stamp")
ixHS6 = identify(HS6$Time, HS6$depth)
plot(CTD$Time, CTD$Depth)
print("Click for the CTD time stamp")
ixCTD = identify(CTD$Time, CTD$Depth)
diffHS6 = difftime(CTD$Time[ixCTD], HS6$Time[ixHS6], units="secs")
HS6$Time = HS6$Time + diffHS6
}
}
plot(HS6$Time, HS6$depth, pch=19,cex=0.5, main="Check depth matching",
sub = "Adjust Time0.CTD in cast.info if needed")
points(CTD$Time, CTD$Depth, col=2, pch=19, cex=0.2)
legend("topright", c("HS6", "CTD"), col=c(1,2), pch=c(19,19))
text(HS6$Time[ix.maxZ.HS6], 0.8*HS6$depth[ix.maxZ.HS6],
paste('Time diff:', as.character(signif(diffHS6,3))))
png("CTDvsHS6_depth.png", height = 5, width = 7, units = "in", res=300)
plot(HS6$Time, HS6$depth, pch=19,cex=0.5, main="Check depth matching",
sub = "Adjust Time0.CTD in cast.info if needed")
points(CTD$Time, CTD$Depth, col=2, pch=19, cex=0.2)
legend("topright", c("HS6", "CTD"), col=c(1,2), pch=c(19,19))
text(HS6$Time[ix.maxZ.HS6], 0.8*HS6$depth[ix.maxZ.HS6],
paste('Time diff:', as.character(signif(diffHS6,3))))
dev.off()
}
if (instrument$FLECO ==1) {
print("FOR ECOFL, hit when the instrument leave or enter the water")
plot(FLECO$Time, FLECO$FL[,3])
print("Click for the FLECO time stamp")
ixFLECO = identify(FLECO$Time, FLECO$FL[,3])
plot(CTD$Time, CTD$Sal)
print("Click for the CTD time stamp")
ixCTD = identify(CTD$Time, CTD$Sal)
diffFLECO = difftime(CTD$Time[ixCTD], FLECO$Time[ixFLECO], units="secs")
FLECO$Time = FLECO$Time + diffFLECO
}
# New cases added when CTD and BB9 or ACS are deployed in Stand alone mode
if (instrument$BB9 == 1 & instrument$CTD.DH4 == 0) {
BB9$Time <- Time0.CTD + BB9$Timer # Here the timer is in second
plot(BB9$Time, BB9$Betau[,2])
print("Click for the BB9 time stamp when exiting the water")
ixBB9 = identify(BB9$Time, BB9$Betau[,2])
plot(CTD$Time, CTD$Sal)
print("Click for the CTD time stamp when exiting the water")
ixCTD = identify(CTD$Time, CTD$Sal)
diffBB9 = difftime(CTD$Time[ixCTD], BB9$Time[ixBB9], units="secs")
BB9$Time = BB9$Time + diffBB9
}
if (instrument$ACS == 1 & instrument$CTD.DH4 == 0) {
ACS$Time <- Time0.CTD + ACS$Timer/1000 # Here the timer is in millsecond
plot(ACS$Time, ACS$c[,10])
print("Click for the ACS time stamp when exiting the water")
ixACS = identify(ACS$Time, ACS$c[,10])
plot(CTD$Time, CTD$Sal)
print("Click for the CTD time stamp when exiting the water")
ixCTD = identify(CTD$Time, CTD$Sal)
diffACS = difftime(CTD$Time[ixCTD], ACS$Time[ixACS], units="secs")
ACS$Time = ACS$Time + diffACS
}
dt = difftime(max(CTD$Time), min(CTD$Time), units="secs")
# use a 10 seconds time window to define the span parameter
span.CTD = max(10 / as.numeric(dt), 0.1)
}
if (is.null(CTD$Time)) {
print("No CTD data; no time to adjust")
### Check if HS6 and ASPH are present and adjust the depth...
if (instrument$ASPH == 1 && instrument$HS6 == 1) {
ix.maxZ.ASPH = which.max(ASPH$depth)
ix.maxZ.HS6 = which.max(HS6$depth)
diffASPH = difftime(HS6$Time[ix.maxZ.HS6], ASPH$time[ix.maxZ.ASPH], units="secs")
ASPH$time = ASPH$time + diffASPH
plot(ASPH$time, ASPH$depth, pch=19,cex=0.5, main="Check depth matching",
sub = "Adjust Time0.CTD in cast.info if needed")
points(HS6$Time, HS6$depth, col=2, pch=19, cex=0.2)
legend("topright", c("ASPH", "HS6"), col=c(1,2), pch=c(19,19))
text(ASPH$time[ix.maxZ.ASPH], 0.8*ASPH$depth[ix.maxZ.ASPH],
paste('Time diff:', as.character(signif(diffASPH,3))))
}
}
####### Add Depth to instruments without pressure sensors#####
# Match BB9, BB3, FLBBCD, FLECO and ACs to obtain depth from the CTD using time stamp.
if (instrument$ACS == 1 & !is.null(CTD$Time)) {
print("Add depth to ACS")
df = as.data.frame(cbind(CTD$Time, CTD$Depth))
names(df) = c("t","z")
mod = loess(z ~ t, data=df, span=span.CTD)
ACS$Depth = predict(mod, ACS$Time)
}
if (instrument$BB9 == 1 & !is.null(CTD$Time)){
print("Add depth to BB9")
df = as.data.frame(cbind(CTD$Time, CTD$Depth))
names(df) = c("t","z")
mod = loess(z ~ t, data=df, span=span.CTD)
BB9$Depth = predict(mod, BB9$Time)
}
if (instrument$BB3 == 1 & !is.null(CTD$Time)){
print("Add depth to BB3")
df = as.data.frame(cbind(CTD$Time, CTD$Depth))
names(df) = c("t","z")
mod = loess(z ~ t, data=df, span=span.CTD)
BB3$Depth = predict(mod, BB3$Time)
}
if (instrument$FLECO == 1 & !is.null(CTD$Time)) {
print("Add depth to FLECO")
df = as.data.frame(cbind(CTD$Time, CTD$Depth))
names(df) = c("t","z")
mod = loess(z ~ t, data=df, span=span.CTD)
FLECO$Depth = predict(mod, FLECO$Time)
}
if (instrument$FLECO == 1 & is.null(CTD$Time)) {
print("NO CTD available to set the depth to FLECO")
print("Check for depth from another instrument")
if (instrument$HS6 == 1) {
df = as.data.frame(cbind(HS6$Time, HS6$depth))
names(df) = c("t","z")
dt = difftime(max(HS6$Time), min(HS6$Time), units="secs")
# use a 10 seconds time window to define the span parameter
span.HS6 = max(10 / as.numeric(dt), 0.1)
mod = loess(z ~ t, data=df, span=span.HS6)
FLECO$Depth = predict(mod, FLECO$Time)
}
if (instrument$HS6 == 0 && instrument$ASPH == 1) {
df = as.data.frame(cbind(ASPH$time, ASPH$depth))
names(df) = c("t","z")
dt = difftime(max(ASPH$time), min(ASPH$time), units="secs")
# use a 10 seconds time window to define the span parameter
span.ASPH = max(10 / as.numeric(dt), 0.1)
mod = loess(z ~ t, data=df, span=span.ASPH)
FLECO$Depth = predict(mod, FLECO$Time)
}
if (instrument$HS6 == 0 && instrument$ASPH == 0) {
print("No depth sensor available for FLECO")
print("Set instrument to 0")
instrument$FLECO <- 0
}
}
if (instrument$FLBBCD == 1 & !is.null(CTD$Time)) {
print("Add depth to FLBBCD")
df = as.data.frame(cbind(CTD$Time, CTD$Depth))
names(df) = c("t","z")
mod = loess(z ~ t, data=df, span=span.CTD)
FLBBCD$Depth = predict(mod, FLBBCD$Time)
}
if (instrument$FLCHL == 1 & !is.null(CTD$Time)) {
print("Add depth to FLCHL")
df = as.data.frame(cbind(CTD$Time, CTD$Depth))
names(df) = c("t","z")
mod = loess(z ~ t, data=df, span=span.CTD)
FLCHL$Depth = predict(mod, FLCHL$Time)
}
####### Add Temperature and/or Salinity to instruments for various correction #####
# Add Temperature and Salinity in the ACS file for TS correction
if (instrument$ACS == 1 & !is.null(CTD$Time)) {
df = as.data.frame(cbind(CTD$Time, CTD$Temp))
names(df) = c("t","T")
mod = loess(T ~ t, data=df, span=span.CTD)
ACS$T = predict(mod, ACS$Time)
df = as.data.frame(cbind(CTD$Time, CTD$Sal))
names(df) = c("t","S")
mod = loess(S ~ t, data=df, span=span.CTD)
ACS$S = predict(mod, ACS$Time)
}
# Add Temperature and Salinity in the ASPH file for TS correction
if (instrument$ASPH == 1 & !is.null(CTD$Time)) {
df = as.data.frame(cbind(CTD$Time, CTD$Temp))
names(df) = c("t","T")
mod = loess(T ~ t, data=df, span=span.CTD)
ASPH$T = predict(mod, ASPH$time)
df = as.data.frame(cbind(CTD$Time, CTD$Sal))
names(df) = c("t","S")
mod = loess(S ~ t, data=df, span=span.CTD)
ASPH$S = predict(mod, ASPH$time)
}
# Add Salinity in the BB9 file for betaP calculation
if (instrument$BB9 == 1 & !is.null(CTD$Time)) {
df = as.data.frame(cbind(CTD$Time, CTD$Sal))
names(df) = c("t","S")
mod = loess(S ~ t, data=df, span=span.CTD)
BB9$S = predict(mod, BB9$Time)
}
# Add Salinity in the BB3 file for betaP calculation
if (instrument$BB3 == 1 & !is.null(CTD$Time)) {
df = as.data.frame(cbind(CTD$Time, CTD$Sal))
names(df) = c("t","S")
mod = loess(S ~ t, data=df, span=span.CTD)
BB3$S = predict(mod, BB3$Time)
}
# Add Salinity in the FLBBCD file for betaP calculation
if (instrument$FLBBCD == 1 & !is.null(CTD$Time)) {
df = as.data.frame(cbind(CTD$Time, CTD$Sal))
names(df) = c("t","S")
mod = loess(S ~ t, data=df, span=span.CTD)
FLBBCD$S = predict(mod, FLBBCD$Time)
}
# Add Salinity in the HS6 file for betaP calculation
if (instrument$HS6 == 1 & !is.null(CTD$Time)) {
df = as.data.frame(cbind(CTD$Time, CTD$Sal))
names(df) = c("t","S")
mod = loess(S ~ t, data=df, span=span.CTD)
HS6$S = predict(mod, HS6$Time)
}
# This case exceptions have been writen to manage data with no CTD
if (instrument$ASPH == 1 & is.null(CTD$Time)) {
Salinity = as.numeric(readline("No CTD data; Enter a salinity value for beta correction of ASPH: "))
ASPH$S = rep(Salinity, length(ASPH$time))
Temperature = as.numeric(readline("No CTD data; Enter a temperature value for beta correction of ASPH: "))
ASPH$T = rep(Temperature, length(ASPH$time))
}
if (instrument$ACS == 1 & is.null(CTD$Time)) {
Salinity = as.numeric(readline("No CTD data; Enter a salinity value for beta correction of ACS: "))
ACS$S = rep(Salinity, length(ACS$Time))
Temperature = as.numeric(readline("No CTD data; Enter a Temperature value for beta correction of ACS: "))
ACS$T = rep(Temperature, length(ACS$Time))
}
if (instrument$HS6 == 1 & is.null(CTD$Time)) {
Salinity = as.numeric(readline("No CTD data; Enter a salinity value for beta correction of HS6: "))
HS6$S = rep(Salinity, length(HS6$Time))
}
if (instrument$BB9 == 1 & is.null(CTD$Time)) {
Salinity = as.numeric(readline("No CTD data; Enter a salinity value for beta correction of BB9: "))
BB9$S = rep(Salinity, length(BB9$Time))
}
if (instrument$BB3 == 1 & is.null(CTD$Time)) {
Salinity = as.numeric(readline("No CTD data; Enter a salinity value for beta correction of BB3: "))
BB3$S = rep(Salinity, length(BB3$Time))
}
if (instrument$FLBBCD == 1 & is.null(CTD$Time)) {
Salinity = as.numeric(readline("No CTD data; Enter a salinity value for beta correction of FLBBCD: "))
FLBBCD$S = rep(Salinity, length(FLBBCD$Time))
}
####### Apply TS correction to ASPH and/or ACS ##############################
if (instrument$ASPH == 1) {
# create a matrix of coefficients for ASPH
print("Apply TS correction to ASPH")
PsiT = spline(Tdf$wavelength, Tdf$PsiT, xout=ASPH$wl,method="natural")$y
PsiS = spline(Sdf$wavelength, Sdf$PsiS_ASW, xout=ASPH$wl,method="natural")$y
PsiTm = matrix(PsiT, nrow=length(ASPH$depth), ncol=length(ASPH$wl), byrow=T)
delta_T = (ASPH$T-Tref.ASPH)
delta_Tm = matrix(delta_T, nrow=length(ASPH$depth), ncol=length(ASPH$wl), byrow=F)
Tcor_factor = PsiTm * delta_Tm
PsiSm = matrix(PsiS, nrow=length(ASPH$depth), ncol=length(ASPH$wl), byrow=T)
Sm = matrix(ASPH$S, nrow=length(ASPH$depth), ncol=length(ASPH$wl), byrow=F)
Scor_factor = PsiSm * Sm
# Apply correction to ASPH
ASPH$a.TScor = ASPH$a - Tcor_factor - Scor_factor
}
if (instrument$ACS ==1) {
print("Apply TS correction to ACS a and c beams")
PsiT = spline(TS4.cor.df$wavelength, TS4.cor.df$PsiT, xout=ACS$a.wl,method="natural")$y
PsiS = spline(TS4.cor.df$wavelength, TS4.cor.df$PsiSa, xout=ACS$a.wl,method="natural")$y
PsiSc = spline(TS4.cor.df$wavelength, TS4.cor.df$PsiSc, xout=ACS$c.wl,method="natural")$y
PsiTm = matrix(PsiT, nrow=length(ACS$Depth), ncol=length(ACS$a.wl), byrow=T)
delta_T = (ACS$T-Tref.ACS)
delta_Tm = matrix(delta_T, nrow=length(ACS$Depth), ncol=length(ACS$a.wl), byrow=F)
Tcor_factor = PsiTm * delta_Tm
# create a matrix of coefficients for ACs
PsiSm = matrix(PsiS, nrow=length(ACS$Depth), ncol=length(ACS$a.wl), byrow=T)
PsiScm = matrix(PsiSc, nrow=length(ACS$Depth), ncol=length(ACS$c.wl), byrow=T)
Sm = matrix(ACS$S, nrow=length(ACS$Depth), ncol=length(ACS$a.wl), byrow=T)
Scor_factor.a = PsiSm * Sm
Scor_factor.c = PsiScm * Sm
# Apply correction to ACS
ACS$a.TScor = ACS$a - Tcor_factor - Scor_factor.a
ACS$c.TScor = ACS$c - Tcor_factor - Scor_factor.c
}
####### Apply blank correction to ASC and ASPH####
if (instrument$ASPH == 1) {
if (file.exists(parameters$blank.ASPH)){
print("Apply blank correction to ASPH")
load(file=parameters$blank.ASPH)
ASPH.blank.m = matrix(ASPH.blank$smooth, nrow=length(ASPH$depth), ncol=length(ASPH$wl), byrow=T)
ASPH$a.corrected = ASPH$a.TScor -ASPH.blank.m
} else {
print("No blank correction to ASPH")
ASPH$a.corrected = ASPH$a.TScor
}
}
if (instrument$ACS ==1) {
if (file.exists(parameters$blank.ACS)){
print("Apply blank correction to ACS")
load(file=parameters$blank.ACS)
ACS.blank.a.m = matrix(ACS.blank$a.smooth, nrow=length(ACS$Depth), ncol=length(ACS$a.wl), byrow=T)
ACS.blank.c.m = matrix(ACS.blank$c.smooth, nrow=length(ACS$Depth), ncol=length(ACS$c.wl), byrow=T)
ACS$a.TScor.offset = ACS$a.TScor - ACS.blank.a.m
ACS$c.corrected = ACS$c.TScor - ACS.blank.c.m
} else {
print("No blank correction to ACS")
ACS$a.TScor.offset = ACS$a.TScor
ACS$c.corrected = ACS$c.TScor
}
####### Apply scattering correction to ACS ####
print("Apply scattering correction to ACS absorption")
ix.NIR = which(ACS$c.wl > 720)
a.nir = apply(ACS$a.TScor.offset[,ix.NIR],1,mean, na.rm=T)
c.nir = apply(ACS$c.corrected[,ix.NIR],1,mean, na.rm=T)
#a.nir[a.nir < 0] = NA
#c.nir[c.nir < 0] = NA
#a.nir[a.nir > 5] = NA
#c.nir[c.nir > 50] = NA
# Mckee et al Opt. Express 2008 method
if (scat.correction == "mckee") {
if (instrument$BB9 ==1) {
# First correct BB9 data using ACS absorption corrected using zaneveld
b.nir = c.nir - a.nir
scat.factor = a.nir / b.nir
scat.factor[scat.factor < 0] = 0
scat.factor.m = matrix(scat.factor, nrow=length(ACS$Depth), ncol=length(ACS$a.wl), byrow=F)
b_star = ACS$c.corrected - ACS$a.TScor.offset
ACS$a.corrected = ACS$a.TScor.offset - (scat.factor.m * b_star)
# correct bb9 data
# Compute Pure Water beta
S.fac = (1+(0.3*BB9$S/37))*1e-4 * (1+ cos(117/180*pi)*cos(117/180*pi)*((1-0.09)/(1+0.09)))
wl.fac = 1.38*(BB9$waves/500)^-4.32
S.fac.m = matrix(S.fac, nrow=length(BB9$Depth), ncol=length(BB9$waves), byrow=F)
wl.fac.m = matrix(wl.fac, nrow=length(BB9$Depth), ncol=length(BB9$waves),byrow=T)
BetaW = S.fac.m * wl.fac.m
bbW = (0.0029308*(BB9$waves/500)^-4.24)/2
bbW.m = matrix(bbW, nrow=length(BB9$Depth), ncol=length(BB9$waves), byrow=T)
Beta2bb <- 2*pi*1.1
ix.BB9.wl = rep(NA,9)
for (j in 1:9){
ix.BB9.wl[j] = which.min(abs(ACS$a.wl-BB9$waves[j]))
}
dt = difftime(max(BB9$Time), min(BB9$Time), units="secs")
# use a 10 seconds time window to define the span parameter
span.BB9 = max(10 / as.numeric(dt), 0.1)
tmp = fit.with.loess(ACS$a.wl[ix.BB9.wl], as.numeric(ACS$Time),
ACS$a.corrected[,ix.BB9.wl], span.BB9, as.numeric(BB9$Time))
a.bb9 = tmp$aop.fitted
# New pathlenght estimated by Doxaran et al 2016
BB9$Beta.corrected <- BB9$Betau * exp(0.01635*a.bb9)
BB9$BetaP.corrected <- BB9$Beta.corrected - BetaW
BB9$bbP.corrected <- Beta2bb * BB9$BetaP.corrected
# Extrapolate bbp to match ACS data matrix.
# Only use the extremity to avoid extrapolation artefact
# A nls fit could be implemented eventually
good.ix <- which(!is.na(BB9$bbP.corrected[,1]))
bbP.acs = t(apply(BB9$bbP.corrected[good.ix,c(2,9)],1,
function(y){
res = spline(BB9$waves[c(2,9)],y,xout=ACS$c.wl, method="natural")$y
}))
# now interpolate for depth
dt = difftime(max(ACS$Time), min(ACS$Time), units="secs")
span.ACS = max(10 / as.numeric(dt), 0.1)
tmp = fit.with.loess(ACS$c.wl, as.numeric(BB9$Time[good.ix]),
bbP.acs, span.ACS, as.numeric(ACS$Time))
bbP.acs = tmp$aop.fitted # overwrite
bbp.tild.0 = bbP.acs / (1.5*b_star) # from Mckee et al. (Fig. 4)
fa0 = apply(bbp.tild.0, 2, function(x){
y = 2.699e-3 + 4.636*x - 3.746e1 *x^2 + 3.177e2*x^3 -1.166e3*x^4})
fc0 = apply(bbp.tild.0, 2, function(x){
y = 6.809e-3/(8.502e-3+x)-(1.918e-2)})
bp1 = b_star / (1 - fa0 - fc0) # from Mckee et al. (eq. 7)
bbp.tild.1 = bbP.acs / bp1 # get a first approximation of bbp tild
for (j in 1:10) { # Replace the condition by a fixed number of iterations.
# This is to speed up the processing
prev.bbp.tild = bbp.tild.1
fa1 = apply(bbp.tild.1, 2, function(x){
y = 2.699e-3 + 4.636*x - 3.746e1 *x^2 + 3.177e2*x^3 -1.166e3*x^4})
fc1 = apply(bbp.tild.1, 2, function(x){
y = 6.809e-3/(8.502e-3+x)-(1.918e-2)})
an = ACS$a.TScor.offset - (fa1*b_star)/(1-fa1-fc1)
cn = ACS$c.corrected + (fc1*b_star)/(1-fa1-fc1)
bp2 = cn - an
bbp.tild.1 = bbP.acs / bp2
}
ACS$a.corrected = an
ACS$c.corrected = cn
ACS$scat.correction.method = "mckee"
} else {
# if (instrument$HS6 ==1) {
# TO BE IMPLEMENTED
# } else {
print("No backscattering measurements to apply Mckee et al method")
print("Automatical change correction to Zaneveld method")
scat.correction = "zaneveld"
# }
}
}
# Zaneveld method
if (scat.correction == "zaneveld") {
b.nir = c.nir - a.nir
scat.factor = a.nir / b.nir
scat.factor.m = matrix(scat.factor, nrow=length(ACS$Depth), ncol=length(ACS$a.wl), byrow=F)
b_star = ACS$c.corrected - ACS$a.TScor.offset
ACS$a.corrected = ACS$a.TScor.offset - (scat.factor.m * b_star)
ACS$scat.correction.method = "zaneveld"
}
# Baseline method
if (scat.correction == "baseline") {
a.nir.m = matrix(a.nir, nrow=length(ACS$Depth), ncol=length(ACS$a.wl), byrow=F)
ACS$a.corrected = ACS$a.TScor.offset - a.nir.m
ACS$scat.correction.method = "baseline"
}
if (scat.correction == "none") {
ACS$a.corrected = ACS$a.TScor.offset
ACS$scat.correction.method = "none"
}
# ACS$a.corrected[ACS$a.corrected < 0 ] = NA
ACS$a.corrected[ACS$a.corrected > 5 ] = NA
ACS$c.corrected[ACS$c.corrected < 0 ] = NA
ACS$c.corrected[ACS$c.corrected > 50 ] = NA
}
####### Apply blank and attenuation corrections to BB9 #####
if (instrument$BB9 ==1) {
if (file.exists(parameters$blank.BB9)) {
if (BB9.DEVICE.FILE) {
print("@@@@@@@@@@@@@@@@@@@@")
print("Device file provided instead of a blank.")
print("In that case, to account for field BB9 blank (or dark measurement)")
print("Please change the offset in the device file...")
print("@@@@@@@@@@@@@@@@@@@@")
BB9$bbP.offset = BB9$bbPu
BB9$bb.offset = BB9$bbu
BB9$BetaP.offset = BB9$BetaPu
BB9$Beta.offset = BB9$Betau
BB9$blank = NA
} else {
print("Apply Blank correction to BB9")
load(file=parameters$blank.BB9)
bb.blank.m = matrix(blank.BB9, nrow=length(BB9$Timer), ncol=length(BB9$waves), byrow=T)
beta.blank.m = matrix(blank.BB9/6.9115, nrow=length(BB9$Timer), ncol=length(BB9$waves), byrow=T)
BB9$bbP.offset = BB9$bbPu - bb.blank.m
BB9$bb.offset = BB9$bbu - bb.blank.m
BB9$BetaP.offset = BB9$BetaPu - beta.blank.m
BB9$Beta.offset = BB9$Betau - beta.blank.m
BB9$blank = blank.BB9
}
} else {
print("No blank correction for BB9")
BB9$bbP.offset = BB9$bbPu
BB9$bb.offset = BB9$bbu
BB9$BetaP.offset = BB9$BetaPu
BB9$Beta.offset = BB9$Betau
BB9$blank = NA
}
############################## BB correction for attenuation
print("Apply attenuation correction to BB9")
# Compute Pure Water beta
S.fac = (1+(0.3*BB9$S/37))*1e-4 * (1+ cos(117/180*pi)*cos(117/180*pi)*((1-0.09)/(1+0.09)))
wl.fac = 1.38*(BB9$waves/500)^-4.32
S.fac.m = matrix(S.fac, nrow=length(BB9$Depth), ncol=length(BB9$waves), byrow=F)
wl.fac.m = matrix(wl.fac, nrow=length(BB9$Depth), ncol=length(BB9$waves),byrow=T)
BetaW = S.fac.m * wl.fac.m
bbW = (0.0029308*(BB9$waves/500)^-4.24)/2
bbW.m = matrix(bbW, nrow=length(BB9$Depth), ncol=length(BB9$waves), byrow=T)
Beta2bb <- 2*pi*1.1
if (instrument$ASPH == 1) {
ix.BB9.wl = rep(NA,9)
for (j in 1:9){
ix.BB9.wl[j] = which.min(abs(ASPH$wl - BB9$waves[j]))
}
dt = difftime(max(BB9$Time), min(BB9$Time), units="secs")
# use a 10 seconds time window to define the span parameter
span.BB9 = max(10 / as.numeric(dt), 0.1)
tmp = fit.with.loess(ASPH$wl[ix.BB9.wl], as.numeric(ASPH$time),
ASPH$a.corrected[,ix.BB9.wl], span.BB9, as.numeric(BB9$Time))
a.bb9 = tmp$aop.fitted
# New pathlenght estimated by Doxaran et al 2016
BB9$Beta.corrected <- BB9$Beta.offset * exp(0.01635*a.bb9) # exp(0.0391*a.bb9)
BB9$BetaP.corrected <- BB9$Beta.corrected - BetaW
BB9$bbP.corrected <- Beta2bb * BB9$BetaP.corrected
BB9$bb.corrected <- BB9$bbP.corrected + bbW.m
BB9$a <- a.bb9
BB9$sigma.correction = TRUE
}
if (instrument$ASPH == 0 & instrument$ACS == 1) {
ix.BB9.wl = rep(NA,9)
for (j in 1:9){
ix.BB9.wl[j] = which.min(abs(ACS$a.wl-BB9$waves[j]))
}
dt = difftime(max(BB9$Time), min(BB9$Time), units="secs")
# use a 10 seconds time window to define the span parameter
span.BB9 = max(max(10 / as.numeric(dt), 0.1), 0.1)
tmp = fit.with.loess(ACS$a.wl[ix.BB9.wl], ACS$Timer,
ACS$a.corrected[,ix.BB9.wl], span.BB9, BB9$Timer)
a.bb9 = tmp$aop.fitted
# clean the data in case of outliers
#a.bb9[a.bb9 > 1] = 0
a.bb9[a.bb9 < 0] = 0
# New pathlenght estimated by Doxaran et al 2016
BB9$Beta.corrected <- BB9$Beta.offset * exp(0.01635*a.bb9)
BB9$BetaP.corrected <- BB9$Beta.corrected - BetaW
BB9$bbP.corrected <- Beta2bb * BB9$BetaP.corrected
BB9$bb.corrected <- BB9$bbP.corrected + bbW.m
BB9$a <- a.bb9
BB9$sigma.correction = TRUE
}
if (instrument$ASPH == 0 & instrument$ACS == 0) {
print("WARNING: No attenuation correction to BB-9 data")
BB9$Beta.corrected <- BB9$Beta.offset
BB9$BetaP.corrected = BB9$Beta.corrected - BetaW
BB9$bbP.corrected = Beta2bb * BB9$BetaP.corrected
BB9$bb.corrected = BB9$bbP.corrected + bbW.m
BB9$a <- NA
BB9$sigma.correction = FALSE
}
}
####### Apply blank and attenuation corrections to BB3 #####
if (instrument$BB3 ==1) {
if (file.exists(parameters$blank.BB3)) {
print("Apply Blank correction to BB3")
load(file=parameters$blank.BB3)
bb.blank.m = matrix(blank.BB3, nrow=length(BB3$Timer), ncol=length(BB3$wl), byrow=T)
beta.blank.m = matrix(blank.BB3/6.9115, nrow=length(BB3$Timer), ncol=length(BB3$wl), byrow=T)
BB3$bbP.offset = BB3$bbP - bb.blank.m
BB3$bb.offset = BB3$bb - bb.blank.m
BB3$BetaP.offset = BB3$BetaP - beta.blank.m
BB3$Beta.offset = BB3$Beta - beta.blank.m
BB3$blank = blank.BB3
} else {
print("No blank correction for BB3")
BB3$bbP.offset = BB3$bbP
BB3$bb.offset = BB3$bb
BB3$BetaP.offset = BB3$BetaP
BB3$Beta.offset = BB3$Beta
BB3$blank = NA
}
############################## BB correction for attenuation
print("Apply attenuation correction to BB3")
# Compute Pure Water beta
S.fac = (1+(0.3*BB3$S/37))*1e-4 * (1+ cos(117/180*pi)*cos(117/180*pi)*((1-0.09)/(1+0.09)))
wl.fac = 1.38*(BB3$wl/500)^-4.32
S.fac.m = matrix(S.fac, nrow=length(BB3$Depth), ncol=length(BB3$wl), byrow=F)
wl.fac.m = matrix(wl.fac, nrow=length(BB3$Depth), ncol=length(BB3$wl),byrow=T)
BetaW = S.fac.m * wl.fac.m
bbW = (0.0029308*(BB3$wl/500)^-4.24)/2
bbW.m = matrix(bbW, nrow=length(BB3$Depth), ncol=length(BB3$wl), byrow=T)
Beta2bb <- 2*pi*1.1
if (instrument$ASPH == 1) {
ix.BB3.wl = rep(NA,3)
for (j in 1:3){
ix.BB3.wl[j] = which.min(abs(ASPH$wl - BB3$wl[j]))
}
dt = difftime(max(BB3$Time), min(BB3$Time), units="secs")
# use a 10 seconds time window to define the span parameter
span.BB3 = max(10 / as.numeric(dt), 0.1)
tmp = fit.with.loess(ASPH$wl[ix.BB3.wl], as.numeric(ASPH$time),
ASPH$a.corrected[,ix.BB3.wl], span.BB3, as.numeric(BB3$Time))
a.BB3 = tmp$aop.fitted
# New pathlenght estimated by Doxaran et al 2016
BB3$Beta.corrected <- BB3$Beta.offset * exp(0.01635*a.BB3) # exp(0.0391*a.BB3)
BB3$BetaP.corrected <- BB3$Beta.corrected - BetaW
BB3$bbP.corrected <- Beta2bb * BB3$BetaP.corrected
BB3$bb.corrected <- BB3$bbP.corrected + bbW.m
BB3$a <- a.BB3
BB3$sigma.correction = TRUE
}
if (instrument$ASPH == 0 & instrument$ACS == 1) {
ix.BB3.wl = rep(NA,3)
for (j in 1:3){
ix.BB3.wl[j] = which.min(abs(ACS$a.wl - BB3$wl[j]))
}
dt = difftime(max(BB3$Time), min(BB3$Time), units="secs")
# use a 10 seconds time window to define the span parameter
span.BB3 = max(10 / as.numeric(dt), 0.1)
tmp = fit.with.loess(ACS$a.wl[ix.BB3.wl], ACS$Timer,
ACS$a.corrected[,ix.BB3.wl], span.BB3, BB3$Timer)
a.BB3 = tmp$aop.fitted
# New pathlenght estimated by Doxaran et al 2016
BB3$Beta.corrected <- BB3$Beta.offset * exp(0.01635*a.BB3)
BB3$BetaP.corrected <- BB3$Beta.corrected - BetaW
BB3$bbP.corrected <- Beta2bb * BB3$BetaP.corrected
BB3$bb.corrected <- BB3$bbP.corrected + bbW.m
BB3$a <- a.BB3
BB3$sigma.correction = TRUE
}
if (instrument$ASPH == 0 & instrument$ACS == 0) {
print("WARNING: No attenuation correction to BB-9 data")
BB3$Beta.corrected <- BB3$Beta.offset
BB3$BetaP.corrected = BB3$Beta.corrected - BetaW
BB3$bbP.corrected = Beta2bb * BB3$BetaP.corrected
BB3$bb.corrected = BB3$bbP.corrected + bbW.m
BB3$a <- NA
BB3$sigma.correction = FALSE
}
}
####### Apply attenuation correction to HS6 data #####
if (instrument$HS6 ==1) {
print("Apply attenuation correction to HS6")
# FROM CALIBRATION FILE
if (HS6.CALYEAR == 2010) Sigmaexp <- c(0.139, 0.142, 0.146, 0.143, 0.146, 0.146)
if (HS6.CALYEAR == 2013) Sigmaexp <- c(0.135, 0.141, 0.143, 0.14, 0.142, 0.144)
if (HS6.CALYEAR == 2014) Sigmaexp <- c(0.137, 0.139, 0.144, 0.141, 0.141, 0.143)
if (HS6.CALYEAR == 2016) Sigmaexp <- c(0.137, 0.141, 0.145, 0.14, 0.144, 0.144)
if (HS6.CALYEAR == 2018) Sigmaexp <- c(0.137, 0.142, 0.144, 0.141, 0.143, 0.145)
Sigmaexp.m = matrix(Sigmaexp, nrow=length(HS6$depth), ncol=length(HS6$wl), byrow=T)
Beta2bb <- 6.79
# Compute Pure Water beta
S.fac = (1+(0.3*HS6$S/37))*1e-4 * (1+ cos(140/180*pi)*cos(140/180*pi)*((1-0.09)/(1+0.09)))
wl.fac = 1.38*(HS6$wl/500)^-4.32
S.fac.m = matrix(S.fac, nrow=length(HS6$depth), ncol=length(HS6$wl), byrow=F)
wl.fac.m = matrix(wl.fac, nrow=length(HS6$depth), ncol=length(HS6$wl),byrow=T)
BetaW = S.fac.m * wl.fac.m
bbW = (0.0029308*(HS6$wl/500)^-4.24)/2
bbW.m = matrix(bbW, nrow=length(HS6$depth), ncol=length(HS6$wl), byrow=T)
ix.HS6.wl = rep(NA,6)
# Match APHERE to HS
if (instrument$ASPH == 1) {
for (j in 1:6){
ix.HS6.wl[j] = which(ASPH$wl == HS6$wl[j])
}
dt = difftime(max(HS6$Time), min(HS6$Time), units="secs")
# use a 10 seconds time window to define the span parameter
span.HS6 = 10 /as.numeric(dt)
tmp = fit.with.loess(ASPH$wl[ix.HS6.wl], as.numeric(ASPH$time),
ASPH$a.corrected[,ix.HS6.wl], span.HS6, as.numeric(HS6$Time))
a.HS6 = tmp$aop.fitted
# New correction from Doxaran et al, Opt. Express 2016
Beta.measured <- HS6$betau
Beta.corrected <- Beta.measured
for (k in 1:5) {
Kbb <- a.HS6 + 4.34 * Beta2bb * (Beta.corrected - BetaW)
sigmaKbb <- exp(Sigmaexp.m*Kbb)
Beta.corrected <- sigmaKbb * Beta.measured
}
HS6$Beta.corrected <- Beta.corrected
HS6$BetaP.corrected <- Beta.corrected - BetaW
HS6$bbP.corrected <- Beta2bb * HS6$BetaP.corrected
HS6$bb.corrected <- HS6$bbP.corrected + bbW.m
HS6$a <- a.HS6
HS6$sigma.correction = TRUE
}
if (instrument$ASPH == 0 & instrument$ACS == 1) {
for (j in 1:6){
ix.HS6.wl[j] = which(ACS$a.wl == HS6$wl[j])
}
dt = difftime(max(HS6$Time), min(HS6$Time), units="secs")
# use a 10 seconds time window to define the span parameter
span.HS6 = 10 /as.numeric(dt)
tmp = fit.with.loess(ACS$a.wl[ix.HS6.wl], as.numeric(ACS$Time),
ACS$a.corrected[,ix.HS6.wl], span.HS6, as.numeric(HS6$Time))
a.HS6 = tmp$aop.fitted
# New correction from Doxaran et al, Opt. Express 2016
Beta.measured <- HS6$betau
Beta.corrected <- Beta.measured
for (k in 1:5) {
Kbb <- a.HS6 + 4.34 * Beta2bb * (Beta.corrected - BetaW)
sigmaKbb <- exp(Sigmaexp.m*Kbb)
Beta.corrected <- sigmaKbb * Beta.measured
}
HS6$Beta.corrected <- Beta.corrected
HS6$BetaP.corrected <- Beta.corrected - BetaW
HS6$bbP.corrected <- Beta2bb * HS6$BetaP.corrected
HS6$bb.corrected <- HS6$bbP.corrected + bbW.m
HS6$a <- a.HS6
HS6$sigma.correction = TRUE
}
if (instrument$ASPH == 0 & instrument$ACS == 0) {
print("WARNING: No attenuation correction to HS-6 data")
HS6$Beta.corrected <- HS6$betau
HS6$BetaP.corrected = HS6$Beta.corrected - BetaW
HS6$bbP.corrected = Beta2bb*HS6$BetaP.corrected
HS6$bb.corrected = HS6$bbP.corrected + bbW.m
HS6$a <- NA
HS6$sigma.correction = FALSE
}
}
####### Compute number of particules from LISST #####
if (instrument$LISST == 1) {
print("compute number of particules")
df = cbind(LISST$lower_size_bins*1e-6,LISST$upper_size_bins*1e-6)
D = exp(apply(log(df),1,mean))
Dm = matrix(D, nrow=length(LISST$Depth), ncol=32, byrow=T)
LISST$N_prime = (LISST$PSD*6)/pi/Dm^3
}
print(paste("Bin the data at a given depth interval:", Zint))
####### Select the portion of the profile to retain for depth interpolation#####
if (is.na(minx) & !is.null(CTD$Time)){
plot(CTD$Time, CTD$Depth)
print("Click for the begining of the cast and then ESC")
minx <- identify(CTD$Time, CTD$Depth)
print("Click for the end of the cast (at surface) and then ESC")
maxx <- identify(CTD$Time, CTD$Depth)
}
if (!is.null(CTD$Time)) {
start.profile = CTD$Time[minx]
stop.profile = CTD$Time[maxx]
CTD$ixmin = minx
CTD$ixmax = maxx
z.max = max(CTD$Depth[minx:maxx])
ix.z.max = which(CTD$Depth == z.max)[1]
CTD$ix.z.max = ix.z.max
} else {
print("No CTD data")
# This case expection was implemented for HS6 profile with no CTD
if (instrument$HS6 == 1) {
if (is.na(minx)) {
plot(HS6$Time, HS6$depth)
print("Click for the begining of the cast and then ESC")
minx <- identify(HS6$Time, HS6$depth)
print("Click for the end of the cast (at surface) and then ESC")
maxx <- identify(HS6$Time, HS6$depth)
}
start.profile = HS6$Time[minx]
stop.profile = HS6$Time[maxx]
HS6$ixmin = minx
HS6$ixmax = maxx
z.max = max(HS6$depth[minx:maxx])
ix.z.max = which(HS6$depth == z.max)[1]
HS6$ix.z.max = ix.z.max
}
}
# Set the span parameter
span.CTD.profile = depth.span/z.max
span.ASPH = (depth.span/z.max)*1.5
span.ACS = depth.span/z.max
span.BB9 = depth.span/z.max
span.BB3 = depth.span/z.max
span.LISST = depth.span/z.max
span.HS6 = depth.span/z.max
span.FLECO = depth.span/z.max
span.FLBBCD = depth.span/z.max
span.FLCHL = depth.span/z.max
time.window=c(start.profile, stop.profile)
depth.fitted=seq(Zint,z.max,Zint)
####### Interpolate CTD data ######
if (!is.null(CTD$Time)) {
print("Bin CTD....")
CTD.fitted.down=list()
df = as.data.frame(cbind(CTD$Depth[minx:ix.z.max], CTD$Temp[minx:ix.z.max]))
names(df) = c("z","T")
mod = loess(T ~ z, data=df, span=span.CTD.profile)
CTD.fitted.down$Temp = predict(mod, depth.fitted)
df = as.data.frame(cbind(CTD$Depth[minx:ix.z.max], CTD$Sal[minx:ix.z.max]))
names(df) = c("z","S")
mod = loess(S ~ z, data=df, span=span.CTD.profile)
CTD.fitted.down$Sal = predict(mod, depth.fitted)
if ((maxx - CTD$ix.z.max) < 10) {
print("No CTD upcast to fit")
CTD.fitted.up=list()
} else {
CTD.fitted.up=list()
df = as.data.frame(cbind(CTD$Depth[maxx:ix.z.max], CTD$Temp[maxx:ix.z.max]))
names(df) = c("z","T")
mod = loess(T ~ z, data=df, span=span.CTD.profile)
CTD.fitted.up$Temp = predict(mod, depth.fitted)
df = as.data.frame(cbind(CTD$Depth[maxx:ix.z.max], CTD$Sal[maxx:ix.z.max]))
names(df) = c("z","S")
mod = loess(S ~ z, data=df, span=span.CTD.profile)
CTD.fitted.up$Sal = predict(mod, depth.fitted)
}
} else {
CTD.fitted.down=list()
CTD.fitted.up=list()
}
####### Interpolate ASPH data#####
if (instrument$ASPH == 1) {
print("Bin ASPH....")
ASPH.fitted.down=list()
ASPH.fitted.up=list()
########## Get ASPH index
# find the index to start and stop profile
abs.diff.time = as.numeric(abs(difftime(ASPH$time, time.window[1])))
ASPH$ixmin = which(abs.diff.time == min(abs.diff.time))
abs.diff.time = as.numeric(abs(difftime(ASPH$time, time.window[2])))
ASPH$ixmax = which(abs.diff.time == min(abs.diff.time))
z.max = max(ASPH$depth[ASPH$ixmin:ASPH$ixmax])
ix.z.max = which(ASPH$depth == z.max)[1]
ASPH$ix.z.max = ix.z.max
ix.good.fitted.depth = which(depth.fitted < z.max)
nz=length(depth.fitted)
depth.raw = ASPH$depth[ASPH$ixmin:ix.z.max]
data.raw = ASPH$a.corrected[ASPH$ixmin:ix.z.max,]
bad.row = which(is.infinite(data.raw), arr.ind = T)[,1]
if (length(bad.row) > 0) tmp=fit.with.loess(ASPH$wl, depth.raw[-bad.row],data.raw[-bad.row,],
span=span.ASPH, depth.fitted[ix.good.fitted.depth])
else tmp=fit.with.loess(ASPH$wl, depth.raw, data.raw,
span=span.ASPH, depth.fitted[ix.good.fitted.depth])
ASPH.fitted.down$a = matrix(nrow=nz, ncol=length(ASPH$wl))
ASPH.fitted.down$a[ix.good.fitted.depth,] = tmp$aop.fitted
# remove extrapolated data for down cast
min.z = min(ASPH$depth[ASPH$ixmin:ix.z.max], na.rm=T)
max.z = max(ASPH$depth[ASPH$ixmin:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) ASPH.fitted.down$a[ix.bad,] = NA
if ((ASPH$ixmax - ix.z.max) < 4) {
print("No ASPH upcast to fit")
} else {
tmp=fit.with.loess(ASPH$wl, ASPH$depth[ASPH$ixmax:ix.z.max], ASPH$a.corrected[ASPH$ixmax:ix.z.max,],
span=span.ASPH, depth.fitted)
ASPH.fitted.up$a = tmp$aop.fitted
# remove extrapolated data for up cast
min.z = min(ASPH$depth[ASPH$ixmax:ix.z.max], na.rm=T)
max.z = max(ASPH$depth[ASPH$ixmax:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) ASPH.fitted.up$a[ix.bad,] = NA
}
} else {
ASPH.fitted.down = list()
ASPH.fitted.up = list()
}
####### Interpolate ACS data#####
if (instrument$ACS == 1) {
print("Bin ACS....")
ACS.fitted.down=list()
ACS.fitted.up=list()
########## Get ACS index
# find the index to start and stop profile
abs.diff.time = as.numeric(abs(difftime(ACS$Time, time.window[1])))
ACS$ixmin = which(abs.diff.time == min(abs.diff.time))
abs.diff.time = as.numeric(abs(difftime(ACS$Time, time.window[2])))
ACS$ixmax = which(abs.diff.time == min(abs.diff.time))
z.max = max(ACS$Depth[ACS$ixmin:ACS$ixmax],na.rm=T)
ix.z.max = which(ACS$Depth == z.max)[1]
ACS$ix.z.max = ix.z.max
# Interpolate using loess function
tmp=fit.with.loess(ACS$a.wl, ACS$Depth[ACS$ixmin:ix.z.max], ACS$a.corrected[ACS$ixmin:ix.z.max,],
span=span.ACS, depth.fitted)
ACS.fitted.down$a = tmp$aop.fitted
tmp=fit.with.loess(ACS$c.wl, ACS$Depth[ACS$ixmin:ix.z.max], ACS$c.corrected[ACS$ixmin:ix.z.max,],
span=span.ACS, depth.fitted)
ACS.fitted.down$c = tmp$aop.fitted
# remove extrapolated data for down cast
min.z = min(ACS$Depth[ACS$ixmin:ix.z.max], na.rm=T)
max.z = max(ACS$Depth[ACS$ixmin:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
ACS.fitted.down$a[ix.bad,] = NA
ACS.fitted.down$c[ix.bad,] = NA
}
if ((ACS$ixmax - ix.z.max) < 5) {
print("No ACS upcast to fit")
} else {
tmp=fit.with.loess(ACS$a.wl, ACS$Depth[ACS$ixmax:ix.z.max], ACS$a.corrected[ACS$ixmax:ix.z.max,],
span=span.ACS, depth.fitted)
ACS.fitted.up$a = tmp$aop.fitted
tmp=fit.with.loess(ACS$c.wl, ACS$Depth[ACS$ixmax:ix.z.max], ACS$c.corrected[ACS$ixmax:ix.z.max,],
span=span.ACS, depth.fitted)
ACS.fitted.up$c = tmp$aop.fitted
# remove extrapolated data for up cast
min.z = min(ACS$Depth[ACS$ixmax:ix.z.max], na.rm=T)
max.z = max(ACS$Depth[ACS$ixmax:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
ACS.fitted.up$a[ix.bad,] = NA
ACS.fitted.up$c[ix.bad,] = NA
}
}
} else {
ACS.fitted.down = list()
ACS.fitted.up = list()
}
####### Interpolate BB9 data#####
if (instrument$BB9 == 1){
print("Bin BB9....")
print("HELLO LE BUGGGG")
BB9.fitted.down=list()
BB9.fitted.up=list()
########## Get BB9 index
# find the index to start and stop profile
abs.diff.time = as.numeric(abs(difftime(BB9$Time, time.window[1])))
BB9$ixmin = which(abs.diff.time == min(abs.diff.time))[1]
abs.diff.time = as.numeric(abs(difftime(BB9$Time, time.window[2])))
BB9$ixmax = which(abs.diff.time == min(abs.diff.time))
z.max = max(BB9$Depth[BB9$ixmin:BB9$ixmax])
ix.z.max = which(BB9$Depth == z.max)[1]
BB9$ix.z.max = ix.z.max
print(BB9$ixmin)
str(BB9)
tmp=fit.with.loess(BB9$waves, BB9$Depth[BB9$ixmin:ix.z.max],
BB9$bbP.corrected[BB9$ixmin:ix.z.max,],
span=span.BB9, depth.fitted)
BB9.fitted.down$bbP = tmp$aop.fitted
tmp=fit.with.loess(BB9$waves, BB9$Depth[BB9$ixmin:ix.z.max],
BB9$bb.corrected[BB9$ixmin:ix.z.max,],
span=span.BB9, depth.fitted)
BB9.fitted.down$bb = tmp$aop.fitted
# remove extrapolated data for down cast
min.z = min(BB9$Depth[BB9$ixmin:ix.z.max], na.rm=T)
max.z = max(BB9$Depth[BB9$ixmin:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
BB9.fitted.down$bb[ix.bad,] = NA
BB9.fitted.down$bbP[ix.bad,] = NA
}
#### Compute bbp spectral slope for down cast
x = 555/BB9$waves
nz=length(depth.fitted)
bbP555=rep(0,nz)
nuP=rep(0,nz)
for (i in 1:nz) {
y = BB9.fitted.down$bbP[i,]
y[y < 0] = NA
if (any(is.na(y))) {
bbP555[i]=NA
nuP[i]=NA
} else {
z = nls(y~b*x^z, start = list(b = y[5], z = 0.5),data=data.frame(x,y), control=list(maxiter=100, warnOnly=T))
bbP555[i]=coef(z)[1]
nuP[i]=coef(z)[2]
}
}
BB9.fitted.down$bbP555 = bbP555
BB9.fitted.down$nuP = nuP
if ((BB9$ixmax - ix.z.max) < 5) {
print("No BB9 upcast to fit")
} else {
tmp=fit.with.loess(BB9$waves, BB9$Depth[BB9$ixmax:ix.z.max],
BB9$bbP.corrected[BB9$ixmax:ix.z.max,],
span=span.BB9, depth.fitted)
BB9.fitted.up$bbP = tmp$aop.fitted
tmp=fit.with.loess(BB9$waves, BB9$Depth[BB9$ixmax:ix.z.max],
BB9$bb.corrected[BB9$ixmax:ix.z.max,],
span=span.BB9, depth.fitted)
BB9.fitted.up$bb = tmp$aop.fitted
# remove extrapolated data for up cast
min.z = min(BB9$Depth[BB9$ixmax:ix.z.max], na.rm=T)
max.z = max(BB9$Depth[BB9$ixmax:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
BB9.fitted.up$bb[ix.bad,] = NA
BB9.fitted.up$bbP[ix.bad,] = NA
}
}
} else {
BB9.fitted.down = list()
BB9.fitted.up = list()
}
####### Interpolate BB3 data#####
if (instrument$BB3 == 1){
print("Bin BB3....")
BB3.fitted.down=list()
BB3.fitted.up=list()
########## Get BB3 index
# find the index to start and stop profile
abs.diff.time = as.numeric(abs(difftime(BB3$Time, time.window[1])))
BB3$ixmin = which(abs.diff.time == min(abs.diff.time))
abs.diff.time = as.numeric(abs(difftime(BB3$Time, time.window[2])))
BB3$ixmax = which(abs.diff.time == min(abs.diff.time))
z.max = max(BB3$Depth[BB3$ixmin:BB3$ixmax])
ix.z.max = which(BB3$Depth == z.max)[1]
BB3$ix.z.max = ix.z.max
tmp=fit.with.loess(BB3$wl, BB3$Depth[BB3$ixmin:ix.z.max],
BB3$bbP.corrected[BB3$ixmin:ix.z.max,],
span=span.BB3, depth.fitted)
BB3.fitted.down$bbP = tmp$aop.fitted
tmp=fit.with.loess(BB3$wl, BB3$Depth[BB3$ixmin:ix.z.max],
BB3$bb.corrected[BB3$ixmin:ix.z.max,],
span=span.BB3, depth.fitted)
BB3.fitted.down$bb = tmp$aop.fitted
# remove extrapolated data for down cast
min.z = min(BB3$Depth[BB3$ixmin:ix.z.max], na.rm=T)
max.z = max(BB3$Depth[BB3$ixmin:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
BB3.fitted.down$bb[ix.bad,] = NA
BB3.fitted.down$bbP[ix.bad,] = NA
}
#### Compute bbp spectral slope for downward cast
x = 555/BB3$wl
nz=length(depth.fitted)
bbP555=rep(0,nz)
nuP=rep(0,nz)
for (i in 1:nz) {
y = BB3.fitted.down$bbP[i,]
y[y < 0] = NA
if (any(is.na(y))) {
bbP555[i]=NA
nuP[i]=NA
} else {
z = nls(y~b*x^z, start = list(b = y[1], z = 0.5),data=data.frame(x,y), control=list(maxiter=100, warnOnly=T))
bbP555[i]=coef(z)[1]
nuP[i]=coef(z)[2]
}
}
BB3.fitted.down$bbP555 = bbP555
BB3.fitted.down$nuP = nuP
if ((BB3$ixmax - ix.z.max) < 5) {
print("No BB3 upcast to fit")
} else {
tmp=fit.with.loess(BB3$wl, BB3$Depth[BB3$ixmax:ix.z.max],
BB3$bbP.corrected[BB3$ixmax:ix.z.max,],
span=span.BB3, depth.fitted)
BB3.fitted.up$bbP = tmp$aop.fitted
tmp=fit.with.loess(BB3$wl, BB3$Depth[BB3$ixmax:ix.z.max],
BB3$bb.corrected[BB3$ixmax:ix.z.max,],
span=span.BB3, depth.fitted)
BB3.fitted.up$bb = tmp$aop.fitted
# remove extrapolated data for up cast
min.z = min(BB3$Depth[BB3$ixmax:ix.z.max], na.rm=T)
max.z = max(BB3$Depth[BB3$ixmax:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
BB3.fitted.up$bb[ix.bad,] = NA
BB3.fitted.up$bbP[ix.bad,] = NA
}
}
} else {
BB3.fitted.down = list()
BB3.fitted.up = list()
}
####### Interpolate LISST data#####
if (instrument$LISST == 1){
print("Bin LISST....")
LISST.fitted.down=list()
LISST.fitted.up=list()
########## Get LISST index
# find the index to start and stop profile
abs.diff.time = as.numeric(abs(difftime(LISST$Time, time.window[1])))
LISST$ixmin = which(abs.diff.time == min(abs.diff.time))
abs.diff.time = as.numeric(abs(difftime(LISST$Time, time.window[2])))
LISST$ixmax = which(abs.diff.time == min(abs.diff.time))
z.max = max(LISST$Depth[LISST$ixmin:LISST$ixmax])
ix.z.max = which(LISST$Depth == z.max)[1]
LISST$ix.z.max = ix.z.max
tmp=fit.with.loess(LISST$median_bins, LISST$Depth[LISST$ixmin:ix.z.max],
LISST$PSD[LISST$ixmin:ix.z.max,],
span=span.LISST, depth.fitted)
LISST.fitted.down$PSD = tmp$aop.fitted
tmp=fit.with.loess(LISST$median_bins, LISST$Depth[LISST$ixmin:ix.z.max],
LISST$N_prime[LISST$ixmin:ix.z.max,],
span=span.LISST, depth.fitted)
LISST.fitted.down$N_prime = tmp$aop.fitted
df = as.data.frame(cbind(LISST$Depth[LISST$ixmin:ix.z.max],
LISST$c670[LISST$ixmin:ix.z.max]))
names(df) = c("z","c670")
mod = loess(c670 ~ z, data=df, span=span.LISST)
LISST.fitted.down$c670 = predict(mod, depth.fitted)
# remove extrapolated data for down cast
min.z = min(LISST$Depth[LISST$ixmin:ix.z.max], na.rm=T)
max.z = max(LISST$Depth[LISST$ixmin:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
LISST.fitted.down$PSD[ix.bad,] = NA
LISST.fitted.down$N_prime[ix.bad,] = NA
LISST.fitted.down$c670[ix.bad] = NA
}
if ((LISST$ixmax - ix.z.max) < 5) {
print("No LISST upcast to fit")
} else {
tmp=fit.with.loess(LISST$median_bins, LISST$Depth[LISST$ixmax:ix.z.max],
LISST$PSD[LISST$ixmax:ix.z.max,],
span=span.LISST, depth.fitted)
LISST.fitted.up$PSD = tmp$aop.fitted
tmp=fit.with.loess(LISST$median_bins, LISST$Depth[LISST$ixmax:ix.z.max],
LISST$N_prime[LISST$ixmax:ix.z.max,],
span=span.LISST, depth.fitted)
LISST.fitted.up$N_prime = tmp$aop.fitted
df = as.data.frame(cbind(LISST$Depth[LISST$ixmax:ix.z.max],
LISST$c670[LISST$ixmax:ix.z.max]))
names(df) = c("z","c670")
mod = loess(c670 ~ z, data=df, span=span.LISST)
LISST.fitted.up$c670 = predict(mod, depth.fitted)
# remove extrapolated data for up cast
min.z = min(LISST$Depth[LISST$ixmax:ix.z.max], na.rm=T)
max.z = max(LISST$Depth[LISST$ixmax:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
LISST.fitted.up$PSD[ix.bad,] = NA
LISST.fitted.up$N_prime[ix.bad,] = NA
LISST.fitted.up$c670[ix.bad] = NA
}
}
} else {
LISST.fitted.down=list()
LISST.fitted.up=list()
}
####### Interpolate HS6 data######
if (instrument$HS6 ==1) {
print("Bin HS6....")
HS6.fitted.down=list()
HS6.fitted.up=list()
########## Get HS6 index
# find the index to start and stop profile
abs.diff.time = as.numeric(abs(difftime(HS6$Time, time.window[1])))
HS6$ixmin = which(abs.diff.time == min(abs.diff.time))[1]
abs.diff.time = as.numeric(abs(difftime(HS6$Time, time.window[2])))
HS6$ixmax = which(abs.diff.time == min(abs.diff.time))[1]
z.max = max(HS6$depth[HS6$ixmin:HS6$ixmax])
ix.z.max = which(HS6$depth == z.max)[1]
HS6$ix.z.max = ix.z.max
# Fit down cast
ix.good.fitted.depth = which(depth.fitted < z.max)
nz=length(depth.fitted)
depth.raw = HS6$depth[HS6$ixmin:ix.z.max]
data.raw = HS6$bbP.corrected[HS6$ixmin:ix.z.max,]
bad.row = which(is.infinite(data.raw), arr.ind = T)[,1]
if (length(bad.row) > 0) tmp=fit.with.loess(HS6$wl, depth.raw[-bad.row],data.raw[-bad.row,],
span=span.HS6, depth.fitted[ix.good.fitted.depth]) else
tmp=fit.with.loess(HS6$wl, depth.raw,data.raw,
span=span.HS6, depth.fitted[ix.good.fitted.depth])
HS6.fitted.down$bbP = matrix(nrow=nz, ncol=6)
HS6.fitted.down$bbP[ix.good.fitted.depth,] = tmp$aop.fitted
data.raw = HS6$bb.corrected[HS6$ixmin:ix.z.max,]
bad.row = which(is.infinite(data.raw), arr.ind = T)[,1]
if (length(bad.row) > 0) tmp=fit.with.loess(HS6$wl, depth.raw[-bad.row],data.raw[-bad.row,],
span=span.HS6, depth.fitted[ix.good.fitted.depth]) else
tmp=fit.with.loess(HS6$wl, depth.raw,data.raw,
span=span.HS6, depth.fitted[ix.good.fitted.depth])
HS6.fitted.down$bb = matrix(nrow=nz, ncol=6)
HS6.fitted.down$bb[ix.good.fitted.depth,] = tmp$aop.fitted
# Fluorescence channels
data.raw = HS6$fluo[HS6$ixmin:ix.z.max,]
bad.row = which(is.infinite(data.raw), arr.ind = T)[,1]
if (length(bad.row) > 0) tmp=fit.with.loess(c(470,700), depth.raw[-bad.row], data.raw[-bad.row,],
span=span.HS6, depth.fitted[ix.good.fitted.depth]) else
tmp=fit.with.loess(c(470,700), depth.raw,data.raw,
span=span.HS6, depth.fitted[ix.good.fitted.depth])
HS6.fitted.down$FDOM = rep(NA,nz)
HS6.fitted.down$FDOM[ix.good.fitted.depth] = tmp$aop.fitted[,1]
HS6.fitted.down$FCHL = rep(NA,nz)
HS6.fitted.down$FCHL[ix.good.fitted.depth] = tmp$aop.fitted[,2]
# remove extrapolated data for down cast
min.z = min(HS6$depth[HS6$ixmin:ix.z.max], na.rm=T)
max.z = max(HS6$depth[HS6$ixmin:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
HS6.fitted.down$FDOM[ix.bad] = NA
HS6.fitted.down$FCHL[ix.bad] = NA
HS6.fitted.down$bbP[ix.bad,] = NA
HS6.fitted.down$bb[ix.bad,] = NA
}
#### Compute bbp and bb spectral slope
x = 555/HS6$wl
nz=length(depth.fitted)
bbP555.down =rep(0,nz)
nuP.down =rep(0,nz)
for (i in 1:nz) {
y = HS6.fitted.down$bbP[i,]
y[y < 0] = NA
y[y > 0.1] = NA
if (any(is.na(y))) {
bbP555.down[i]=NA
nuP.down[i]=NA
} else {
model = nls(y~b*x^z, start = list(b = y[3], z = 1),data=data.frame(x,y), control=list(maxiter=100, warnOnly=T))
bbP555.down[i]=coef(model)[1]
nuP.down[i]=coef(model)[2]
}
}
HS6.fitted.down$bbP555 = bbP555.down
HS6.fitted.down$nuP = nuP.down
if ((HS6$ixmax - ix.z.max) < 5) {
print("No HS6 upcast to fit")
} else {
#### Fit the upcast
depth.raw = HS6$depth[HS6$ixmax:ix.z.max]
data.raw = HS6$bbP.corrected[HS6$ixmax:ix.z.max,]
bad.row = which(is.infinite(data.raw), arr.ind = T)[,1]
if (length(bad.row) > 0){
tmp=fit.with.loess(HS6$wl, depth.raw[-bad.row] ,
data.raw[-bad.row,],
span=span.HS6, depth.fitted)
} else {
tmp=fit.with.loess(HS6$wl, depth.raw,
data.raw,
span=span.HS6, depth.fitted)
}
HS6.fitted.up$bbP = tmp$aop.fitted
data.raw = HS6$bb.corrected[HS6$ixmax:ix.z.max,]
bad.row = which(is.infinite(data.raw), arr.ind = T)[,1]
if (length(bad.row) > 0){
tmp=fit.with.loess(HS6$wl, depth.raw[-bad.row] ,
data.raw[-bad.row,],
span=span.HS6, depth.fitted)
} else {
tmp=fit.with.loess(HS6$wl, depth.raw,
data.raw,
span=span.HS6, depth.fitted)
}
HS6.fitted.up$bb = tmp$aop.fitted
data.raw = HS6$fluo[HS6$ixmax:ix.z.max,]
bad.row = which(is.infinite(data.raw), arr.ind = T)[,1]
if (length(bad.row) > 0){
tmp=fit.with.loess(c(470,700), depth.raw[-bad.row] ,
data.raw[-bad.row,],
span=span.HS6, depth.fitted)
} else {
tmp=fit.with.loess(c(470,700), depth.raw,
data.raw,
span=span.HS6, depth.fitted)
}
HS6.fitted.up$FDOM = tmp$aop.fitted[,1]
HS6.fitted.up$FCHL = tmp$aop.fitted[,2]
# remove extrapolated data for up cast
min.z = min(HS6$depth[HS6$ixmax:ix.z.max], na.rm=T)
max.z = max(HS6$depth[HS6$ixmax:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
HS6.fitted.up$FDOM[ix.bad] = NA
HS6.fitted.up$FCHL[ix.bad] = NA
HS6.fitted.up$bbP[ix.bad,] = NA
HS6.fitted.up$bb[ix.bad,] = NA
}
#### Compute bbp and bb spectral slope
x = 555/HS6$wl
nz=length(depth.fitted)
bbP555.up =rep(0,nz)
nuP.up =rep(0,nz)
for (i in 1:nz) {
y = HS6.fitted.up$bbP[i,]
y[y < 0] = NA
y[y > 0.1] = NA
if (any(is.na(y))) {
bbP555.up[i]=NA
nuP.up[i]=NA
} else {
model = nls(y~b*x^z, start = list(b = y[3], z = 1),data=data.frame(x,y), control=list(maxiter=100, warnOnly=T))
bbP555.up[i]=coef(model)[1]
nuP.up[i]=coef(model)[2]
}
}
HS6.fitted.up$bbP555 = bbP555.up
HS6.fitted.up$nuP = nuP.up
}
} else {
HS6.fitted.down = list()
HS6.fitted.up = list()
}
####### Interpolate FLECO data#####
if (instrument$FLECO ==1) {
print("Bin FLECO....")
FLECO.fitted.down=list()
FLECO.fitted.up=list()
########## Get FLECO index
# find the index to start and stop profile
abs.diff.time = as.numeric(abs(difftime(FLECO$Time, time.window[1])))
FLECO$ixmin = which(abs.diff.time == min(abs.diff.time))
abs.diff.time = as.numeric(abs(difftime(FLECO$Time, time.window[2])))
FLECO$ixmax = which(abs.diff.time == min(abs.diff.time))
z.max = max(FLECO$Depth[FLECO$ixmin:FLECO$ixmax])
ix.z.max = which(FLECO$Depth == z.max)[1]
FLECO$ix.z.max = ix.z.max
# Interpolate using loess function
tmp=fit.with.loess(FLECO$em, FLECO$Depth[FLECO$ixmin:ix.z.max], FLECO$FL[FLECO$ixmin:ix.z.max,],
span=span.FLECO, depth.fitted)
FLECO.fitted.down$FL = tmp$aop.fitted
# remove extrapolated data for down cast
min.z = min(FLECO$depth[FLECO$ixmin:ix.z.max], na.rm=T)
max.z = max(FLECO$depth[FLECO$ixmin:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
FLECO.fitted.down$FL[ix.bad,] = NA
}
if ((FLECO$ixmax - ix.z.max) < 5) {
print("No FLECO upcast to fit")
} else {
tmp=fit.with.loess(FLECO$em, FLECO$Depth[FLECO$ixmax:ix.z.max], FLECO$FL[FLECO$ixmax:ix.z.max,],
span=span.FLECO, depth.fitted)
FLECO.fitted.up$FL = tmp$aop.fitted
# remove extrapolated data for up cast
min.z = min(FLECO$depth[FLECO$ixmax:ix.z.max], na.rm=T)
max.z = max(FLECO$depth[FLECO$ixmax:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
FLECO.fitted.up$FL[ix.bad,] = NA
}
}
} else {
FLECO.fitted.down = list()
FLECO.fitted.up = list()
}
####### Interpolate FLBBCD data #####
if (instrument$FLBBCD ==1) {
print("Bin FLBBCD....")
FLBBCD.fitted.down=list()
FLBBCD.fitted.up=list()
########## Get FLBBCD index
# find the index to start and stop profile
abs.diff.time = as.numeric(abs(difftime(FLBBCD$Time, time.window[1])))
FLBBCD$ixmin = which(abs.diff.time == min(abs.diff.time))
abs.diff.time = as.numeric(abs(difftime(FLBBCD$Time, time.window[2])))
FLBBCD$ixmax = which(abs.diff.time == min(abs.diff.time))
z.max = max(FLBBCD$Depth[FLBBCD$ixmin:FLBBCD$ixmax])
ix.z.max = which(FLBBCD$Depth == z.max)[1]
FLBBCD$ix.z.max = ix.z.max
# Interpolate using loess function down cast
depth.raw = FLBBCD$Depth[FLBBCD$ixmin:ix.z.max]
# Fluorescence channels
data.raw = cbind(FLBBCD$FCHL[FLBBCD$ixmin:ix.z.max],FLBBCD$FDOM[FLBBCD$ixmin:ix.z.max])
bad.row = which(is.infinite(data.raw), arr.ind = T)[,1]
if (length(bad.row) > 0){
tmp=fit.with.loess(c(470,700), depth.raw[-bad.row] ,
data.raw[-bad.row,],
span=span.FLBBCD, depth.fitted)
} else {
tmp=fit.with.loess(c(470,700), depth.raw,
data.raw,
span=span.FLBBCD, depth.fitted)
}
FLBBCD.fitted.down$FCHL = tmp$aop.fitted[,1]
FLBBCD.fitted.down$FDOM = tmp$aop.fitted[,2]
# BB channel
data.raw = cbind(FLBBCD$bb700[FLBBCD$ixmin:ix.z.max],FLBBCD$bbP700[FLBBCD$ixmin:ix.z.max])
bad.row = which(is.infinite(data.raw), arr.ind = T)[,1]
if (length(bad.row) > 0){
tmp=fit.with.loess(c(700,700), depth.raw[-bad.row] ,
data.raw[-bad.row,],
span=span.FLBBCD, depth.fitted)
} else {
tmp=fit.with.loess(c(700,700), depth.raw,
data.raw,
span=span.FLBBCD, depth.fitted)
}
FLBBCD.fitted.down$bb700 = tmp$aop.fitted[,1]
FLBBCD.fitted.down$bbP700 = tmp$aop.fitted[,2]
# remove extrapolated data for down cast
min.z = min(FLBBCD$depth[FLBBCD$ixmin:ix.z.max], na.rm=T)
max.z = max(FLBBCD$depth[FLBBCD$ixmin:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
FLBBCD.fitted.down$FCHL[ix.bad] = NA
FLBBCD.fitted.down$FDOM[ix.bad] = NA
FLBBCD.fitted.down$bb700[ix.bad] = NA
FLBBCD.fitted.down$bbp700[ix.bad] = NA
}
if ((FLBBCD$ixmax - ix.z.max) < 5) {
print("No FLBBCD upcast to fit")
} else {
##### up cast
depth.raw = FLBBCD$Depth[FLBBCD$ixmax:ix.z.max]
# Fluorescence channels
data.raw = cbind(FLBBCD$FCHL[FLBBCD$ixmax:ix.z.max],FLBBCD$FDOM[FLBBCD$ixmax:ix.z.max])
bad.row = which(is.infinite(data.raw), arr.ind = T)[,1]
if (length(bad.row) > 0){
tmp=fit.with.loess(c(470,700), depth.raw[-bad.row] ,
data.raw[-bad.row,],
span=span.FLBBCD, depth.fitted)
} else {
tmp=fit.with.loess(c(470,700), depth.raw,
data.raw,
span=span.FLBBCD, depth.fitted)
}
FLBBCD.fitted.up$FCHL = tmp$aop.fitted[,1]
FLBBCD.fitted.up$FDOM = tmp$aop.fitted[,2]
# BB channel
data.raw = cbind(FLBBCD$bb700[FLBBCD$ixmax:ix.z.max],FLBBCD$bbP700[FLBBCD$ixmax:ix.z.max])
bad.row = which(is.infinite(data.raw), arr.ind = T)[,1]
if (length(bad.row) > 0){
tmp=fit.with.loess(c(700,700), depth.raw[-bad.row] ,
data.raw[-bad.row,],
span=span.FLBBCD, depth.fitted)
} else {
tmp=fit.with.loess(c(700,700), depth.raw,
data.raw,
span=span.FLBBCD, depth.fitted)
}
FLBBCD.fitted.up$bb700 = tmp$aop.fitted[,1]
FLBBCD.fitted.up$bbP700 = tmp$aop.fitted[,2]
# remove extrapolated data for down cast
min.z = min(FLBBCD$depth[FLBBCD$ixmax:ix.z.max], na.rm=T)
max.z = max(FLBBCD$depth[FLBBCD$ixmax:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
FLBBCD.fitted.up$FCHL[ix.bad] = NA
FLBBCD.fitted.up$FDOM[ix.bad] = NA
FLBBCD.fitted.up$bb700[ix.bad] = NA
FLBBCD.fitted.up$bbp700[ix.bad] = NA
}
}
} else {
FLBBCD.fitted.down = list()
FLBBCD.fitted.up = list()
}
####### Interpolate FLCHL data #####
if (instrument$FLCHL ==1) {
print("Bin FLCHL...")
FLCHL.fitted.down=list()
FLCHL.fitted.up=list()
########## Get FLCHL index
# find the index to start and stop profile
abs.diff.time = as.numeric(abs(difftime(FLCHL$Time, time.window[1])))
FLCHL$ixmin = which(abs.diff.time == min(abs.diff.time))
abs.diff.time = as.numeric(abs(difftime(FLCHL$Time, time.window[2])))
FLCHL$ixmax = which(abs.diff.time == min(abs.diff.time))
z.max = max(FLCHL$Depth[FLCHL$ixmin:FLCHL$ixmax])
ix.z.max = which(FLCHL$Depth == z.max)[1]
FLCHL$ix.z.max = ix.z.max
# Interpolate using loess function
df = as.data.frame(cbind(FLCHL$Depth[FLCHL$ixmin:ix.z.max], FLCHL$FCHL[FLCHL$ixmin:ix.z.max]))
names(df) = c("z","FCHL")
mod = loess(FCHL ~ z, data=df, span=span.FLCHL)
FLCHL.fitted.down$FCHL = predict(mod, depth.fitted)
# remove extrapolated data for down cast
min.z = min(FLCHL$Depth[FLCHL$ixmin:ix.z.max], na.rm=T)
max.z = max(FLCHL$Depth[FLCHL$ixmin:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
FLCHL.fitted.down$FCHL[ix.bad] = NA
}
if ((FLCHL$ixmax - ix.z.max) < 5) {
print("No FLCHL upcast to fit")
} else {
# Interpolate using loess function
df = as.data.frame(cbind(FLCHL$Depth[FLCHL$ixmax:ix.z.max], FLCHL$FCHL[FLCHL$ixmax:ix.z.max]))
names(df) = c("z","FCHL")
mod = loess(FCHL ~ z, data=df, span=span.FLCHL)
FLCHL.fitted.up$FCHL = predict(mod, depth.fitted)
# remove extrapolated data for up cast
min.z = min(FLCHL$Depth[FLCHL$ixmax:ix.z.max], na.rm=T)
max.z = max(FLCHL$Depth[FLCHL$ixmax:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
FLCHL.fitted.up$FCHL[ix.bad] = NA
}
}
} else {
FLCHL.fitted.down = list()
FLCHL.fitted.up = list()
}
####### PUT TOGETHER ALL PARAMETERS######
IOP = list(time.window=time.window, Time0.CTD = Time0.CTD,
ASPH=ASPH, ACS=ACS, BB9=BB9, CTD=CTD, LISST=LISST,
HS6=HS6, FLECO=FLECO, BB3=BB3, FLBBCD=FLBBCD, FLCHL=FLCHL)
IOP.fitted.down = list(Depth = depth.fitted, ASPH=ASPH.fitted.down,
ACS=ACS.fitted.down,
BB9=BB9.fitted.down,
CTD=CTD.fitted.down,
LISST=LISST.fitted.down,
HS6=HS6.fitted.down,
FLECO=FLECO.fitted.down,
BB3=BB3.fitted.down,
FLBBCD=FLBBCD.fitted.down,
FLCHL=FLCHL.fitted.down)
IOP.fitted.up = list(Depth = depth.fitted, ASPH=ASPH.fitted.up,
ACS=ACS.fitted.up,
BB9=BB9.fitted.up,
CTD=CTD.fitted.up,
LISST=LISST.fitted.up,
HS6=HS6.fitted.up,
FLECO=FLECO.fitted.up,
BB3=BB3.fitted.up,
FLBBCD=FLBBCD.fitted.up,
FLCHL=FLCHL.fitted.up)
print("Save data into RData format: IOP, IOP.fitted.down,IOP.fitted.up")
save(IOP,file="IOP.RData")
save(IOP.fitted.down,file="IOP.fitted.down.RData")
save(IOP.fitted.up,file="IOP.fitted.up.RData")
# update cast.info.file
print("UPDATING FILE: cast.info.dat")
df = data.frame(
parameters$lon,
parameters$lat,
parameters$cast,
Time0.CTD,
Time0.LISST,
minx,
maxx,
parameters$Zint,
parameters$asph.skip,
parameters$maxbb,
parameters$Ndepth.to.plot,
parameters$depth.interval.for.smoothing)
names(df) <- c( "lon", "lat", "cast", "Time0.CTD", "Time0.LISST",
"minx", "maxx", "Zint", "asph.skip",
"maxbb", "Ndepth.to.plot","depth.interval.for.smoothing")
write.table(df, file="cast.info.dat", sep=",", quote=F, row.names = F)
print("Return IOP data object")
return(IOP)
}
belasi01/Riops documentation built on Sept. 5, 2022, 6:38 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions correct.merge.IOP.profile
Documented in correct.merge.IOP.profile
#' Apply corrections and merge IOP data from an optical package
#'
#' Apply corrections and merge data coming from various instruments for a
#' given vertical profile in the water column.
#'
#' @param instrument is a list or a data frame of instruments available. For example,
#' instrument = list(ASPH=1, HS6=1, CTD.UQAR=1)
#' @param parameters is a list or a data frame with the processing
#' parameters provided by the user. When called from
#' \code{\link{process.IOPs}}, the parameters are red in cast.info.dat
#' and cal.info.dat. A detailed description of each parameter is provided
#' below (Details section)
#'
#' @return
#' It retruns a large list with all the raw and corrected data of each instrument
#' The list is also save as IOP.RData format. One or tow additional RData are save
#' to store the interpolated data for the down and up cast respectively
#' (IOP.fitted.down.RData and IOP.fitted.up.RData). The upcast
#' is optional and depends on the user input (parameter maxx).
#'
#' @details
#' This is the heart of the \pkg{Riops} package. All the data are red,
#' matched using a unique time frame, corrected and interpolated on a common
#' grid of equally spaced depth. The detailed description of the correction
#' is not described here (To be done in a future version). Here is a summary
#' of the minimum require to run the code succesfully.
#'
#' The processing parameters are stored in two ascii files named cast.info.dat
#' and cal.info.dat. The former file includes the following parameters
#' \itemize{
#' \item{\strong{lon}} { is the longitude of the station. It is not used in the processing
#' it will be in the PDF reporting (see \code{\link{create.report}})}
#' \item{\strong{lat}} { in the latitude of the station. It is not used in the processing
#' it will be in the PDF reporting (see \code{\link{create.report}})}
#' \item{\strong{cast}} { is a character string corresponding to either MINIDAS cast number
#' or file extension generated by the WAP from DH4 archive file (e.g. "001")}
#' \item{\strong{Time0.CTD}}{ is a character string for CTD start time. It will be transformed
#' into a time stamp using \code{\link{as.POSIXct}} and should respect this format "2015-05-06 12:18:01".
#' However, if NA the program it check in other instruments (e.g. ASPH or BB3) to find the starting
#' time. If no absolute time is found, the user will be asked to add a time stamp manually in the
#' cast.info.dat file)}
#' \item{\strong{Time0.LISST}}{ is a character string for LISST starting time. It will be transformed
#' into a time stamp using \code{\link{as.POSIXct}} and should respect this format "2015-05-06 12:18:01".
#' However, if NA the program it check in other instruments (e.g. ASPH or BB3) to find the starting
#' time. If no absolute time is found, the user will be asked to add a time stamp manually in the
#' cast.info.dat file)}
#' \item{\strong{minx}}{ is the index of the CTD timer corresponding to the start of the profile.
#' If NA, the user will be prompted to select click on the start of the profile on a plot
#' showing the depth versus time. This a call to the function \code{\link{identify}}. Once
#' this is done, the user can adjust the begining of the cast by changing the values in
#' the cast.info.dat file afterward (won't be prompted if minx!=NA).}
#' \item{\strong{maxx}}{ is the index of the CTD corresponding to the End of the profile.
#' Similar to \strong{minx}.}
#' \item{\strong{Zint}}{ is the depth interval for profile smooting
#' (e.g. 1 will result in a fitted profile every 1 meter)}
#' \item{\strong{depth.interval.for.smoothing}}{ is the depth interval, in meter, that will
#' be used to smooth data with the loess function. It will be converted into a "span" value.
#' (see \code{\link{loess}})}
#' \item{\strong{asph.skip}}{ is the number or records to skip at the begining of the a-sphere file.
#' This is necessary only when the a-sphere file contains 2 or more profiles. This
#' happens some time when the MiniDAS reset the cast number. Otherwise 0 will fit
#' most situations.}
#' \item{\strong{maxbb}}{ is a parameter for the backscattering plots generated in the PDF report.}
#' \item{\strong{Ndepth.to.plot}}{ is the number of spectra along the depth profile to
#' plot in the report. 8 to 12 usually is enough}
#'}
#'
#' The calibration information are stored in the cal.info.dat.
#' The file includes the following information
#' (NOTE: any field can be ommited with out problem):
#'
#' \itemize{
#' \item{\strong{Tref.ASPH}} { is the temperature of pure water used by
#' Hobilabs for the ASPH calibration.
#' In 2010 and 2014 it was 13.2,
#' in 2013 it was 19,
#' in 2016 it was 14.4.}
#' \item{\strong{HS6.CALYEAR}}{ is the year of the HS6 calibration.}
#' \item{\strong{Tref.ACS}}{ is the temperature of pure water used by
#' Wetlabs for the ACS calibration. The IML ACs calibration of 2013 was 20.3}
#' \item{\strong{scat.correction}}{ is the method for the scattering correction of AC-S abssorption
#' (eg : "mckee","zaneveld", "baseline","none")}
#' \item{\strong{blank.ASPH}}{ is a string for the path of the blank for ASPH as created by
#' \code{\link{analyse.ASPH.blank}}.}
#' \item{\strong{blank.ACS}}{ is a string for the path of the blank for ACS as created by
#' \code{\link{analyse.ACs.blank}}.}
#' \item{\strong{blank.BB9}}{ is a string for the path of the blank for BB9.}
#' \item{\strong{blank.BB3}}{ is a string for the path of the blank for BB3.}
#' }
#'
#'
#' @author Simon Belanger
#'
#' @export
#'
#'
correct.merge.IOP.profile <- function(instrument, parameters){
####### Get prepared to process data #####
print("Get prepared to process data")
path = parameters$path
cast = parameters$cast
Time0.CTD = parameters$Time0.CTD
Time0.LISST = parameters$Time0.LISST
minx= parameters$minx
maxx= parameters$maxx
Zint= parameters$Zint
asph.skip = parameters$asph.skip
depth.span= parameters$depth.interval.for.smoothing
# Check which instrument is missing in the list
if (is.null(instrument$CTD.UQAR)) instrument$CTD.UQAR <- 0
if (is.null(instrument$HS6)) instrument$HS6 <- 0
if (is.null(instrument$FLECO)) instrument$FLECO <- 0
if (is.null(instrument$ASPH)) instrument$ASPH <- 0
if (is.null(instrument$CTD.DH4)) instrument$CTD.DH4 <- 0
if (is.null(instrument$ACS)) instrument$ACS <- 0
if (is.null(instrument$BB9)) instrument$BB9 <- 0
if (is.null(instrument$BB3)) instrument$BB3 <- 0
if (is.null(instrument$LISST)) instrument$LISST <- 0
if (is.null(instrument$FLBBCD)) instrument$FLBBCD <- 0
if (is.null(instrument$FLCHL)) instrument$FLCHL <- 0
# Check information from cal.info
if (!is.null(parameters$Tref.ASPH)) {
Tref.ASPH = parameters$Tref.ASPH
} else {
if (instrument$ASPH == 1) {
print("WARNING: A-sphere calibration temperature of PW must be provided!!!")
print("Add Tref.ASPH in the cal.info.dat file")
print("Stop processing")
retrun(0)
}
Tref.ASPH = NA
}
if (!is.null(parameters$Tref.ACS) & !is.null(parameters$scat.correction)) {
Tref.ACS = parameters$Tref.ACS
scat.correction = parameters$scat.correction
} else {
if (instrument$ACS == 1) {
print("WARNING: AC-s calibration temperature of PW and scat.correction method must be provided!!!")
print("Add Tref.ACS and scat.correction in the cal.info.dat file")
print("Stop processing")
retrun(0)
}
Tref.ACS = NA
scat.correction = NA
}
if (!is.null(parameters$HS6.CALYEAR) & !is.na(parameters$HS6.CALYEAR)) {
HS6.CALYEAR = parameters$HS6.CALYEAR
print(paste("Hydroscat-6 calibration year is", HS6.CALYEAR))
} else {
print(paste("No Hydroscat-6 calibration year provided"))
if (instrument$HS6 == 1) {
print("WARNING: HS-6 calibration year must be provided!!!")
print("Add HS6.CALYEAR in the cal.info.dat file")
print("Stop processing")
return(0)
}
HS6.CALYEAR = NA
}
if (!is.null(parameters$blank.BB9) & !is.na(parameters$blank.BB9)) {
index_ext = length(unlist(strsplit(parameters$blank.BB9, "[.]")))
ext = unlist(strsplit(parameters$blank.BB9, "[.]"))[index_ext]
if (ext == "dev") BB9.DEVICE.FILE = TRUE else BB9.DEVICE.FILE = FALSE
}
if (!is.null(parameters$blank.BB3) & !is.na(parameters$blank.BB3)) {
index_ext = length(unlist(strsplit(parameters$blank.BB3, "[.]")))
ext = unlist(strsplit(parameters$blank.BB3, "[.]"))[index_ext]
if (ext == "dev") BB3.DEVICE.FILE = TRUE else BB3.DEVICE.FILE = FALSE
}
setwd(path)
####### Read all data#########
print("Read all data")
# This was added in order to process previous instrument file
if (!is.null(instrument$CTD.IML)) instrument$CTD.DH4 = instrument$CTD.IML
if (instrument$CTD.UQAR == 1) {
filen = paste("CTD",cast, ".TXT", sep="")
# If the standard CTD file do not exist,
# check for a cnv file instead
if (file.exists(filen)) {
CTD = read.CTD(filen)
} else {
filen = list.files( pattern = "*\\.cnv")
if (file.exists(filen)) {
print("Reading CNV file using oce read function")
CTD = list()
CTD.OCE = read.ctd.sbe(filen)
# Extract data from the CTD object and store as list
if (is.na(Time0.CTD)) Time0.CTD = as.POSIXct(format(CTD.OCE[["startTime"]], tz="GMT"), tz="GMT")
if (is.na(Time0.CTD)) {
print("No CTD time available")
print("Please add a time string for Time0.CTD")
print("in cast.info.dat")
print("STOP PROCESSING")
return(0)
}
if (is.numeric(CTD.OCE[["time"]])) { # This condition was added because
#the "time" in the cnv files are now output
#in POSIXct by the oce package since 2018.
CTD$Time = Time0.CTD + CTD.OCE[["time"]]
} else CTD$Time = CTD.OCE[["time"]]
CTD$Temp = CTD.OCE[["temperature"]]
CTD$Sal = CTD.OCE[["salinity"]]
if (!is.null(CTD.OCE[["par"]])) CTD$PAR<- CTD.OCE[["par"]]
if (!is.null(CTD.OCE[["turbidity"]])) CTD$turbidity <- CTD.OCE[["turbidity"]]
if (!is.null(CTD.OCE[["fluorescence"]])) CTD$fluorescence <- CTD.OCE[["fluorescence"]]
if (!is.null(CTD.OCE[["depth"]])) CTD$Depth= swDepth(CTD.OCE[["depth"]]) else CTD$Depth= swDepth(CTD.OCE[["pressure"]])
if (!is.null(CTD.OCE[["sigmaTheta"]])) CTD$sigmaTheta=CTD.OCE[["sigmaTheta"]]
}
}
} else {
print("No CTD.UQAR to process")
if (instrument$CTD.DH4 == 1) {
print("Read Port Assossiation for DH4")
default.DH4.ports.file <- paste( Sys.getenv("R_IOPs_DATA_DIR"), "DH4.ports.dat", sep = "/")
DH4.ports.file <- paste(parameters$path, "DH4.ports.dat", sep = "/")
if(!file.exists(DH4.ports.file)) {
file.copy(from = default.DH4.ports.file, to = DH4.ports.file)
cat("EDIT file", DH4.ports.file, "and CUSTOMIZE IT\n")
stop()
}
DH4.ports <- read.table(DH4.ports.file, sep=";", header = T,
colClasses = c("character", "character"))
ix.inst = which(DH4.ports$Sensor == "CTD")
pattern = paste("_", DH4.ports$Port[ix.inst],"_", sep="")
filen = list.files( pattern = pattern)
if(!file.exists(filen)) {
print(paste("CTD file not found in", getwd()))
print("Check whether or not the port number is good")
print("Use capital letters for sensor name")
print(DH4.ports)
stop()
} else CTD = read.CTD.DH4(filen)
} else {
print("No CTD to process")
# print("Abort processing")
CTD = list()
# return(1)
}
}
if (instrument$ASPH == 1) {
filen = list.files( pattern = "CST.*\\.txt")
if (length(filen) == 0) filen = list.files( pattern = "ASPH.*\\.txt")
if (length(filen) > 1) {
ix.filen = which.min(str_length(filen))
filen = filen[ix.filen]
}
ASPH = read.ASPH(filen, skip=asph.skip) # modifiy by Simon Belanger
} else {
print("No ASPH to process")
ASPH = list()
# This condition is added because the CTD from IML has no absolute time information
# So the absolute time is given by the ASPH instrument.
# if (instrument$CTD.DH4 == 1) {
# print("Abort processing")
# return(0)
# }
}
if (instrument$HS6 == 1) {
filen = paste("./HS6",cast,".dat", sep="")
if (!file.exists(filen)) filen = list.files( pattern = ".*\\HS.*\\.dat")
HS6 = read.HS6(filen)
} else {
print("No Hydroscat-6 to process")
HS6 = list()
}
if (instrument$FLECO == 1) {
FLECO = read.FLECO(paste("./FL",cast,".TXT", sep=""))
} else {
print("No fluorescence ECOtriplet to process")
FLECO = list()
}
if (instrument$BB9 == 1) {
# BB9 attached to a DH4...
if (instrument$CTD.DH4 == 1) {
ix.inst = which(DH4.ports$Sensor == "BB9")
pattern = paste("_", DH4.ports$Port[ix.inst],"_", sep="")
filen = list.files( pattern = pattern)
if(!file.exists(filen)) {
print(paste("BB9 file not found in", getwd()))
print("Check whether or not the port number is good")
print("Use capital letters for sensor name")
print(DH4.ports)
stop()
} else BB9 = read.BB9(filen)
# extract the port number
BB9.port = as.numeric(DH4.ports$Port[ix.inst])
} else {
# BB9 in stand alone mode.... data could be provided in calibrated format,
# or in raw format which the prefered method.
pattern = "BB9"
filen = list.files( pattern = pattern)
if(!file.exists(filen)) {
print(paste("BB9 file not found in", getwd()))
print("BB9 should appear in the file name...")
stop()
} else {
# Check the extension. If raw, then apply calibration from the device file.
index_ext = length(unlist(strsplit(filen, "[.]"))) # for station names with periods, ex. G604.5
ext = unlist(strsplit(filen, "[.]"))[index_ext]
if (ext == "raw") {
print("BB9 deployed in Stand Alone mode and raw format")
BB9.raw <- read.BB9(filen, raw = TRUE)
if (BB9.DEVICE.FILE) {
BB9 <- apply.ECO.cal(BB9.raw, dev.file=parameters$blank.BB9, dark.file=NA,ECO.type="BB9")
} else {
print("Raw BB9 detected but no device file provided!!! ")
print("Please put the path and file name of the device file")
print("in the cal.info.dat for the >blank.BB9< parameter")
BB9.DEVICE.FILE == FALSE
}
} else {
print("BB9 deployed in Stand Alone mode and eng format")
BB9 <- read.BB9(filen)
}
}
}
} else {
print("No BB9 to process")
BB9 = list()
}
if (instrument$ACS == 1) {
# ACS attached to a DH4...
if (instrument$CTD.DH4 == 1) {
ix.inst = which(DH4.ports$Sensor == "ACS")
pattern = paste("_", DH4.ports$Port[ix.inst],"_", sep="")
filen = list.files( pattern = pattern)
if(!file.exists(filen)) {
print(paste("ACS file not found in", getwd()))
print("Check whether or not the port number is good")
print("Use capital letters for sensor name")
print(DH4.ports)
stop()
} else ACS = read.ACs(filen)
# extract the port number from name
ACS.port = as.numeric(DH4.ports$Port[ix.inst])
} else {
# ACS in stand alone mode.... data could be provided in calibrated format,
# or in raw format which the prefered method.
pattern = "ACS"
filen = list.files( pattern = pattern)
if(!file.exists(filen)) {
print(paste("ACS file not found in", getwd()))
print("ACS should appear in the file name...")
stop()
} else {
ACS <- read.ACs(filen)
}
}
} else {
print("No ACS to process")
ACS = list()
}
if (instrument$LISST == 1){
ix.inst = which(DH4.ports$Sensor == "LISST")
pattern = paste("_", DH4.ports$Port[ix.inst],"_", sep="")
filen = list.files( pattern = pattern)
if(!file.exists(filen)) {
print("Check for LISST in stand alone mode")
filen = list.files( pattern = "L.*\\.asc")
if (!file.exists(filen)) {
print(paste("LISST file not found in", getwd()))
print("Check whether or not the port number is good")
print("Use capital letters for sensor name")
print(DH4.ports)
stop()
} else LISST = read.LISST(filen)
} else LISST = read.LISST(filen)
} else {
print("No LISST to process")
LISST = list()
}
if (instrument$BB3 == 1){
ix.inst = which(DH4.ports$Sensor == "BB3")
pattern = paste("_", DH4.ports$Port[ix.inst],"_", sep="")
filen = list.files( pattern = pattern)
if(!file.exists(filen)) {
print(paste("BB3 file not found in", getwd()))
print("Check whether or not the port number is good")
print("Use capital letters for sensor name")
print(DH4.ports)
stop()
} else BB3 = read.BB3(filen)
BB3.port = as.numeric(DH4.ports$Port[ix.inst])
} else {
print("No BB3 to process")
BB3 = list()
}
if (instrument$FLBBCD == 1){
ix.inst = which(DH4.ports$Sensor == "FLBBCD")
pattern = paste("_", DH4.ports$Port[ix.inst],"_", sep="")
filen = list.files( pattern = pattern)
if(!file.exists(filen)) {
print(paste("FLBBCD file not found in", getwd()))
print("Check whether or not the port number is good")
print("Use capital letters for sensor name")
print(DH4.ports)
stop()
} else FLBBCD = read.FLBBCD(filen)
FLBBCD.port = as.numeric(DH4.ports$Port[ix.inst])
} else {
print("No FLBBCD to process")
FLBBCD = list()
}
if (instrument$FLCHL == 1){
ix.inst = which(DH4.ports$Sensor == "FLCHL")
pattern = paste("_", DH4.ports$Port[ix.inst],"_", sep="")
filen = list.files( pattern = pattern)
if(!file.exists(filen)) {
print(paste("FLCHL file not found in", getwd()))
print("Check whether or not the port number is good")
print("Use capital letters for sensor name")
print(DH4.ports)
stop()
} else FLCHL = read.FLCHL(filen)
FLCHL.port = as.numeric(DH4.ports$Port[ix.inst])
} else {
print("No FLCHL to process")
FLCHL = list()
}
####### Adjusting time stamp to each sensor; two cases are considered####
######## for 1) IML optical package and 2) UQAR optical package, or SBE CTD
if (instrument$CTD.DH4 == 1) {
# Reads Time OFFSET for each sensor
filen = list.files( pattern = "_TO.")
TO = read.table(filen, header=T)
if (is.na(Time0.CTD) & instrument$ASPH == 1) {
#Time0.CTD = find.Time0.CTD(path,cast)
# Find time0 from ASPH to apply time stamp to CTD, BB9 and ACS
print("Create Time Stamp from ASPHERE")
ix.maxZ.CTD = which.max(CTD$Depth)
ix.maxZ.ASPH = which.max(ASPH$depth)
Time0.CTD = ASPH$time[ix.maxZ.ASPH]- CTD$Timer[ix.maxZ.CTD]/1000
# Plot
plot(CTD$Time, CTD$Depth, pch=19,cex=0.5, main=" depth matching (adjust cast.info if needed)")
points(ASPH$time, ASPH$depth, col=2, pch=19, cex=0.2)
legend("topright", c("CTD", "ASPH"), col=c(1,2), pch=c(19,19))
}
if (is.na(Time0.CTD) & instrument$BB3 == 1) {
#Time0.CTD = find.Time0.CTD(path,cast)
# Find time0 from ASPH to apply time stamp to CTD, BB9 and ACS
print("Create Time Stamp from BB3")
Time0.CTD = BB3$Time[1]
}
if (is.na(Time0.CTD)) {
print("No instrument have absolute time string")
print("You need to provide a POSXCT time string")
print("for Time0.CTD and Time0.LISST in the cast.info.dat file")
print("(e.g.2015-05-06 12:18:01 ")
stop()
}
CTD$Time = Time0.CTD + CTD$Timer/1000
# NOTE the Timer, when using the DH4, are in milliseconds. This is why we divide by 1000
if (instrument$ACS == 1) ACS$Time = Time0.CTD + ACS$Timer/1000 - TO$offset[ACS.port]/1000
if (instrument$BB9 == 1) BB9$Time = Time0.CTD + BB9$Timer/1000 - TO$offset[BB9.port]/1000
if (instrument$BB3 == 1) BB3$Time = Time0.CTD + BB3$Timer/1000 - TO$offset[BB3.port]/1000
if (instrument$FLBBCD == 1) FLBBCD$Time = Time0.CTD + FLBBCD$Timer/1000 - TO$offset[FLBBCD.port]/1000
if (instrument$FLCHL == 1) FLCHL$Time = Time0.CTD + FLCHL$Timer/1000 - TO$offset[FLCHL.port]/1000
dt = difftime(max(CTD$Time), min(CTD$Time), units="secs")
# use a 10 seconds time window to define the span parameter for smoothing data
span.CTD = max(10 / as.numeric(dt), 0.1)
# Find time0 from ASPH to apply time stamp to LISST
if (instrument$LISST == 1) {
if (is.na(Time0.LISST) & instrument$ASPH == 1) {
# Check if LISST is in stand alone mode or attached to DH4
ix.inst = which(DH4.ports$Sensor == "LISST")
if(length(ix.inst) == 0) {
ix.maxZ.LISST = which.max(LISST$Depth)
ix.maxZ.ASPH = which.max(ASPH$depth)
# the 1.02932 is the asquisition rate in seconds of the LISST in stand alone mode.
Time0.LISST = ASPH$time[ix.maxZ.ASPH] - ix.maxZ.LISST*1.029432
LISST$Time = Time0.LISST + (1:length(LISST$time))*1.029432
} else {
ix.maxZ.LISST = which.max(LISST$Depth)
ix.maxZ.ASPH = which.max(ASPH$depth)
Time0.LISST = ASPH$time[ix.maxZ.ASPH] - ix.maxZ.LISST
LISST$Time = Time0.LISST + (1:length(LISST$time))
}
}
ix.inst = which(DH4.ports$Sensor == "LISST")
if(length(ix.inst) == 0) {
# the LISST time step in stand alone mode is 1.0294, while it is equal to 1 when attached to DH4
LISST$Time = Time0.LISST + (1:length(LISST$time))*1.029432
} else LISST$Time = Time0.LISST + (1:length(LISST$time))
plot(LISST$Time, LISST$Depth, pch=19,cex=0.5, main=" depth matching (adjust cast.info if needed)")
points(ASPH$time, ASPH$depth, col=2, pch=19, cex=0.2)
legend("topright", c("LISST", "ASPH"), col=c(1,2), pch=c(19,19))
}
}
if (instrument$CTD.UQAR == 1) {
# Assume CTD has the correct time
if (!is.na(Time0.CTD)) {
diffCTD = difftime(CTD$Time[1], Time0.CTD, units="secs")
CTD$Time = CTD$Time + diffCTD
} else {
Time0.CTD = CTD$Time[1]
diffCTD = difftime(CTD$Time[1], Time0.CTD, units="secs") # so difftime ==0
}
print(diffCTD)
CTD$Time <- CTD$Time + diffCTD
#### modification for processing time series
if (parameters$timeserie == 1) {
print("Correct time difference based on the time when each instruments hit the water")
plot(CTD$Time, CTD$Depth)
print("Click for the CTD time stamp")
ixCTD = identify(CTD$Time, CTD$Depth)
} else {
ix.maxZ.CTD = which.max(CTD$Depth)
}
if (instrument$ASPH == 1) {
ix.maxZ.ASPH = which.max(ASPH$depth)
if (parameters$timeserie == 1) {
plot(ASPH$time, ASPH$depth)
print("Click for the ASPH time stamp")
ixASPH = identify(ASPH$time, ASPH$depth)
diffASPH = difftime(CTD$Time[ixCTD], ASPH$time[ixASPH], units="secs")
} else {
diffASPH = difftime(CTD$Time[ix.maxZ.CTD], ASPH$time[ix.maxZ.ASPH], units="secs")
}
ASPH$time = ASPH$time + diffASPH
plot(ASPH$time, ASPH$depth, pch=19,cex=0.5, main="Check depth matching",
sub = "Adjust Time0.CTD in cast.info if needed")
points(CTD$Time, CTD$Depth, col=2, pch=19, cex=0.2)
legend("topright", c("ASPH", "CTD"), col=c(1,2), pch=c(19,19))
text(ASPH$time[ix.maxZ.ASPH], 0.8*ASPH$depth[ix.maxZ.ASPH],
paste('Time diff:', as.character(signif(diffASPH,3))))
png("CTDvsASPH_depth.png", height = 5, width = 7, units = "in", res=300)
plot(ASPH$time, ASPH$depth, pch=19,cex=0.5, main="Check depth matching",
sub = "Adjust Time0.CTD in cast.info if needed")
points(CTD$Time, CTD$Depth, col=2, pch=19, cex=0.2)
legend("topright", c("ASPH", "CTD"), col=c(1,2), pch=c(19,19))
text(ASPH$time[ix.maxZ.ASPH], 0.8*ASPH$depth[ix.maxZ.ASPH],
paste('Time diff:', as.character(signif(diffASPH,3))))
dev.off()
}
if (instrument$HS6 == 1) {
ix.maxZ.HS6 = which.max(HS6$depth)
if (parameters$timeserie == 1) {
print("Correct time difference based on the time when each instruments hit the water")
plot(HS6$Time, HS6$depth)
print("Click for the HS6 time stamp")
ixHS6 = identify(HS6$Time, HS6$depth)
plot(CTD$Time, CTD$Depth)
print("Click for the CTD time stamp")
ixCTD = identify(CTD$Time, CTD$Depth)
diffHS6 = difftime(CTD$Time[ixCTD], HS6$Time[ixHS6], units="secs")
HS6$Time = HS6$Time + diffHS6
} else {
diff.z = abs(HS6$depth[ix.maxZ.HS6] - CTD$Depth[ix.maxZ.CTD])
if (diff.z < 1) {
diffHS6 = difftime(CTD$Time[ix.maxZ.CTD], HS6$Time[ix.maxZ.HS6], units="secs")
HS6$Time = HS6$Time + diffHS6
} else {
print("Correct time difference based on the time when each instruments hit the water")
plot(HS6$Time, HS6$depth)
print("Click for the HS6 time stamp")
ixHS6 = identify(HS6$Time, HS6$depth)
plot(CTD$Time, CTD$Depth)
print("Click for the CTD time stamp")
ixCTD = identify(CTD$Time, CTD$Depth)
diffHS6 = difftime(CTD$Time[ixCTD], HS6$Time[ixHS6], units="secs")
HS6$Time = HS6$Time + diffHS6
}
}
plot(HS6$Time, HS6$depth, pch=19,cex=0.5, main="Check depth matching",
sub = "Adjust Time0.CTD in cast.info if needed")
points(CTD$Time, CTD$Depth, col=2, pch=19, cex=0.2)
legend("topright", c("HS6", "CTD"), col=c(1,2), pch=c(19,19))
text(HS6$Time[ix.maxZ.HS6], 0.8*HS6$depth[ix.maxZ.HS6],
paste('Time diff:', as.character(signif(diffHS6,3))))
png("CTDvsHS6_depth.png", height = 5, width = 7, units = "in", res=300)
plot(HS6$Time, HS6$depth, pch=19,cex=0.5, main="Check depth matching",
sub = "Adjust Time0.CTD in cast.info if needed")
points(CTD$Time, CTD$Depth, col=2, pch=19, cex=0.2)
legend("topright", c("HS6", "CTD"), col=c(1,2), pch=c(19,19))
text(HS6$Time[ix.maxZ.HS6], 0.8*HS6$depth[ix.maxZ.HS6],
paste('Time diff:', as.character(signif(diffHS6,3))))
dev.off()
}
if (instrument$FLECO ==1) {
print("FOR ECOFL, hit when the instrument leave or enter the water")
plot(FLECO$Time, FLECO$FL[,3])
print("Click for the FLECO time stamp")
ixFLECO = identify(FLECO$Time, FLECO$FL[,3])
plot(CTD$Time, CTD$Sal)
print("Click for the CTD time stamp")
ixCTD = identify(CTD$Time, CTD$Sal)
diffFLECO = difftime(CTD$Time[ixCTD], FLECO$Time[ixFLECO], units="secs")
FLECO$Time = FLECO$Time + diffFLECO
}
# New cases added when CTD and BB9 or ACS are deployed in Stand alone mode
if (instrument$BB9 == 1 & instrument$CTD.DH4 == 0) {
BB9$Time <- Time0.CTD + BB9$Timer # Here the timer is in second
plot(BB9$Time, BB9$Betau[,2])
print("Click for the BB9 time stamp when exiting the water")
ixBB9 = identify(BB9$Time, BB9$Betau[,2])
plot(CTD$Time, CTD$Sal)
print("Click for the CTD time stamp when exiting the water")
ixCTD = identify(CTD$Time, CTD$Sal)
diffBB9 = difftime(CTD$Time[ixCTD], BB9$Time[ixBB9], units="secs")
BB9$Time = BB9$Time + diffBB9
}
if (instrument$ACS == 1 & instrument$CTD.DH4 == 0) {
ACS$Time <- Time0.CTD + ACS$Timer/1000 # Here the timer is in millsecond
plot(ACS$Time, ACS$c[,10])
print("Click for the ACS time stamp when exiting the water")
ixACS = identify(ACS$Time, ACS$c[,10])
plot(CTD$Time, CTD$Sal)
print("Click for the CTD time stamp when exiting the water")
ixCTD = identify(CTD$Time, CTD$Sal)
diffACS = difftime(CTD$Time[ixCTD], ACS$Time[ixACS], units="secs")
ACS$Time = ACS$Time + diffACS
}
dt = difftime(max(CTD$Time), min(CTD$Time), units="secs")
# use a 10 seconds time window to define the span parameter
span.CTD = max(10 / as.numeric(dt), 0.1)
}
if (is.null(CTD$Time)) {
print("No CTD data; no time to adjust")
### Check if HS6 and ASPH are present and adjust the depth...
if (instrument$ASPH == 1 && instrument$HS6 == 1) {
ix.maxZ.ASPH = which.max(ASPH$depth)
ix.maxZ.HS6 = which.max(HS6$depth)
diffASPH = difftime(HS6$Time[ix.maxZ.HS6], ASPH$time[ix.maxZ.ASPH], units="secs")
ASPH$time = ASPH$time + diffASPH
plot(ASPH$time, ASPH$depth, pch=19,cex=0.5, main="Check depth matching",
sub = "Adjust Time0.CTD in cast.info if needed")
points(HS6$Time, HS6$depth, col=2, pch=19, cex=0.2)
legend("topright", c("ASPH", "HS6"), col=c(1,2), pch=c(19,19))
text(ASPH$time[ix.maxZ.ASPH], 0.8*ASPH$depth[ix.maxZ.ASPH],
paste('Time diff:', as.character(signif(diffASPH,3))))
}
}
####### Add Depth to instruments without pressure sensors#####
# Match BB9, BB3, FLBBCD, FLECO and ACs to obtain depth from the CTD using time stamp.
if (instrument$ACS == 1 & !is.null(CTD$Time)) {
print("Add depth to ACS")
df = as.data.frame(cbind(CTD$Time, CTD$Depth))
names(df) = c("t","z")
mod = loess(z ~ t, data=df, span=span.CTD)
ACS$Depth = predict(mod, ACS$Time)
}
if (instrument$BB9 == 1 & !is.null(CTD$Time)){
print("Add depth to BB9")
df = as.data.frame(cbind(CTD$Time, CTD$Depth))
names(df) = c("t","z")
mod = loess(z ~ t, data=df, span=span.CTD)
BB9$Depth = predict(mod, BB9$Time)
}
if (instrument$BB3 == 1 & !is.null(CTD$Time)){
print("Add depth to BB3")
df = as.data.frame(cbind(CTD$Time, CTD$Depth))
names(df) = c("t","z")
mod = loess(z ~ t, data=df, span=span.CTD)
BB3$Depth = predict(mod, BB3$Time)
}
if (instrument$FLECO == 1 & !is.null(CTD$Time)) {
print("Add depth to FLECO")
df = as.data.frame(cbind(CTD$Time, CTD$Depth))
names(df) = c("t","z")
mod = loess(z ~ t, data=df, span=span.CTD)
FLECO$Depth = predict(mod, FLECO$Time)
}
if (instrument$FLECO == 1 & is.null(CTD$Time)) {
print("NO CTD available to set the depth to FLECO")
print("Check for depth from another instrument")
if (instrument$HS6 == 1) {
df = as.data.frame(cbind(HS6$Time, HS6$depth))
names(df) = c("t","z")
dt = difftime(max(HS6$Time), min(HS6$Time), units="secs")
# use a 10 seconds time window to define the span parameter
span.HS6 = max(10 / as.numeric(dt), 0.1)
mod = loess(z ~ t, data=df, span=span.HS6)
FLECO$Depth = predict(mod, FLECO$Time)
}
if (instrument$HS6 == 0 && instrument$ASPH == 1) {
df = as.data.frame(cbind(ASPH$time, ASPH$depth))
names(df) = c("t","z")
dt = difftime(max(ASPH$time), min(ASPH$time), units="secs")
# use a 10 seconds time window to define the span parameter
span.ASPH = max(10 / as.numeric(dt), 0.1)
mod = loess(z ~ t, data=df, span=span.ASPH)
FLECO$Depth = predict(mod, FLECO$Time)
}
if (instrument$HS6 == 0 && instrument$ASPH == 0) {
print("No depth sensor available for FLECO")
print("Set instrument to 0")
instrument$FLECO <- 0
}
}
if (instrument$FLBBCD == 1 & !is.null(CTD$Time)) {
print("Add depth to FLBBCD")
df = as.data.frame(cbind(CTD$Time, CTD$Depth))
names(df) = c("t","z")
mod = loess(z ~ t, data=df, span=span.CTD)
FLBBCD$Depth = predict(mod, FLBBCD$Time)
}
if (instrument$FLCHL == 1 & !is.null(CTD$Time)) {
print("Add depth to FLCHL")
df = as.data.frame(cbind(CTD$Time, CTD$Depth))
names(df) = c("t","z")
mod = loess(z ~ t, data=df, span=span.CTD)
FLCHL$Depth = predict(mod, FLCHL$Time)
}
####### Add Temperature and/or Salinity to instruments for various correction #####
# Add Temperature and Salinity in the ACS file for TS correction
if (instrument$ACS == 1 & !is.null(CTD$Time)) {
df = as.data.frame(cbind(CTD$Time, CTD$Temp))
names(df) = c("t","T")
mod = loess(T ~ t, data=df, span=span.CTD)
ACS$T = predict(mod, ACS$Time)
df = as.data.frame(cbind(CTD$Time, CTD$Sal))
names(df) = c("t","S")
mod = loess(S ~ t, data=df, span=span.CTD)
ACS$S = predict(mod, ACS$Time)
}
# Add Temperature and Salinity in the ASPH file for TS correction
if (instrument$ASPH == 1 & !is.null(CTD$Time)) {
df = as.data.frame(cbind(CTD$Time, CTD$Temp))
names(df) = c("t","T")
mod = loess(T ~ t, data=df, span=span.CTD)
ASPH$T = predict(mod, ASPH$time)
df = as.data.frame(cbind(CTD$Time, CTD$Sal))
names(df) = c("t","S")
mod = loess(S ~ t, data=df, span=span.CTD)
ASPH$S = predict(mod, ASPH$time)
}
# Add Salinity in the BB9 file for betaP calculation
if (instrument$BB9 == 1 & !is.null(CTD$Time)) {
df = as.data.frame(cbind(CTD$Time, CTD$Sal))
names(df) = c("t","S")
mod = loess(S ~ t, data=df, span=span.CTD)
BB9$S = predict(mod, BB9$Time)
}
# Add Salinity in the BB3 file for betaP calculation
if (instrument$BB3 == 1 & !is.null(CTD$Time)) {
df = as.data.frame(cbind(CTD$Time, CTD$Sal))
names(df) = c("t","S")
mod = loess(S ~ t, data=df, span=span.CTD)
BB3$S = predict(mod, BB3$Time)
}
# Add Salinity in the FLBBCD file for betaP calculation
if (instrument$FLBBCD == 1 & !is.null(CTD$Time)) {
df = as.data.frame(cbind(CTD$Time, CTD$Sal))
names(df) = c("t","S")
mod = loess(S ~ t, data=df, span=span.CTD)
FLBBCD$S = predict(mod, FLBBCD$Time)
}
# Add Salinity in the HS6 file for betaP calculation
if (instrument$HS6 == 1 & !is.null(CTD$Time)) {
df = as.data.frame(cbind(CTD$Time, CTD$Sal))
names(df) = c("t","S")
mod = loess(S ~ t, data=df, span=span.CTD)
HS6$S = predict(mod, HS6$Time)
}
# This case exceptions have been writen to manage data with no CTD
if (instrument$ASPH == 1 & is.null(CTD$Time)) {
Salinity = as.numeric(readline("No CTD data; Enter a salinity value for beta correction of ASPH: "))
ASPH$S = rep(Salinity, length(ASPH$time))
Temperature = as.numeric(readline("No CTD data; Enter a temperature value for beta correction of ASPH: "))
ASPH$T = rep(Temperature, length(ASPH$time))
}
if (instrument$ACS == 1 & is.null(CTD$Time)) {
Salinity = as.numeric(readline("No CTD data; Enter a salinity value for beta correction of ACS: "))
ACS$S = rep(Salinity, length(ACS$Time))
Temperature = as.numeric(readline("No CTD data; Enter a Temperature value for beta correction of ACS: "))
ACS$T = rep(Temperature, length(ACS$Time))
}
if (instrument$HS6 == 1 & is.null(CTD$Time)) {
Salinity = as.numeric(readline("No CTD data; Enter a salinity value for beta correction of HS6: "))
HS6$S = rep(Salinity, length(HS6$Time))
}
if (instrument$BB9 == 1 & is.null(CTD$Time)) {
Salinity = as.numeric(readline("No CTD data; Enter a salinity value for beta correction of BB9: "))
BB9$S = rep(Salinity, length(BB9$Time))
}
if (instrument$BB3 == 1 & is.null(CTD$Time)) {
Salinity = as.numeric(readline("No CTD data; Enter a salinity value for beta correction of BB3: "))
BB3$S = rep(Salinity, length(BB3$Time))
}
if (instrument$FLBBCD == 1 & is.null(CTD$Time)) {
Salinity = as.numeric(readline("No CTD data; Enter a salinity value for beta correction of FLBBCD: "))
FLBBCD$S = rep(Salinity, length(FLBBCD$Time))
}
####### Apply TS correction to ASPH and/or ACS ##############################
if (instrument$ASPH == 1) {
# create a matrix of coefficients for ASPH
print("Apply TS correction to ASPH")
PsiT = spline(Tdf$wavelength, Tdf$PsiT, xout=ASPH$wl,method="natural")$y
PsiS = spline(Sdf$wavelength, Sdf$PsiS_ASW, xout=ASPH$wl,method="natural")$y
PsiTm = matrix(PsiT, nrow=length(ASPH$depth), ncol=length(ASPH$wl), byrow=T)
delta_T = (ASPH$T-Tref.ASPH)
delta_Tm = matrix(delta_T, nrow=length(ASPH$depth), ncol=length(ASPH$wl), byrow=F)
Tcor_factor = PsiTm * delta_Tm
PsiSm = matrix(PsiS, nrow=length(ASPH$depth), ncol=length(ASPH$wl), byrow=T)
Sm = matrix(ASPH$S, nrow=length(ASPH$depth), ncol=length(ASPH$wl), byrow=F)
Scor_factor = PsiSm * Sm
# Apply correction to ASPH
ASPH$a.TScor = ASPH$a - Tcor_factor - Scor_factor
}
if (instrument$ACS ==1) {
print("Apply TS correction to ACS a and c beams")
PsiT = spline(TS4.cor.df$wavelength, TS4.cor.df$PsiT, xout=ACS$a.wl,method="natural")$y
PsiS = spline(TS4.cor.df$wavelength, TS4.cor.df$PsiSa, xout=ACS$a.wl,method="natural")$y
PsiSc = spline(TS4.cor.df$wavelength, TS4.cor.df$PsiSc, xout=ACS$c.wl,method="natural")$y
PsiTm = matrix(PsiT, nrow=length(ACS$Depth), ncol=length(ACS$a.wl), byrow=T)
delta_T = (ACS$T-Tref.ACS)
delta_Tm = matrix(delta_T, nrow=length(ACS$Depth), ncol=length(ACS$a.wl), byrow=F)
Tcor_factor = PsiTm * delta_Tm
# create a matrix of coefficients for ACs
PsiSm = matrix(PsiS, nrow=length(ACS$Depth), ncol=length(ACS$a.wl), byrow=T)
PsiScm = matrix(PsiSc, nrow=length(ACS$Depth), ncol=length(ACS$c.wl), byrow=T)
Sm = matrix(ACS$S, nrow=length(ACS$Depth), ncol=length(ACS$a.wl), byrow=T)
Scor_factor.a = PsiSm * Sm
Scor_factor.c = PsiScm * Sm
# Apply correction to ACS
ACS$a.TScor = ACS$a - Tcor_factor - Scor_factor.a
ACS$c.TScor = ACS$c - Tcor_factor - Scor_factor.c
}
####### Apply blank correction to ASC and ASPH####
if (instrument$ASPH == 1) {
if (file.exists(parameters$blank.ASPH)){
print("Apply blank correction to ASPH")
load(file=parameters$blank.ASPH)
ASPH.blank.m = matrix(ASPH.blank$smooth, nrow=length(ASPH$depth), ncol=length(ASPH$wl), byrow=T)
ASPH$a.corrected = ASPH$a.TScor -ASPH.blank.m
} else {
print("No blank correction to ASPH")
ASPH$a.corrected = ASPH$a.TScor
}
}
if (instrument$ACS ==1) {
if (file.exists(parameters$blank.ACS)){
print("Apply blank correction to ACS")
load(file=parameters$blank.ACS)
ACS.blank.a.m = matrix(ACS.blank$a.smooth, nrow=length(ACS$Depth), ncol=length(ACS$a.wl), byrow=T)
ACS.blank.c.m = matrix(ACS.blank$c.smooth, nrow=length(ACS$Depth), ncol=length(ACS$c.wl), byrow=T)
ACS$a.TScor.offset = ACS$a.TScor - ACS.blank.a.m
ACS$c.corrected = ACS$c.TScor - ACS.blank.c.m
} else {
print("No blank correction to ACS")
ACS$a.TScor.offset = ACS$a.TScor
ACS$c.corrected = ACS$c.TScor
}
####### Apply scattering correction to ACS ####
print("Apply scattering correction to ACS absorption")
ix.NIR = which(ACS$c.wl > 720)
a.nir = apply(ACS$a.TScor.offset[,ix.NIR],1,mean, na.rm=T)
c.nir = apply(ACS$c.corrected[,ix.NIR],1,mean, na.rm=T)
#a.nir[a.nir < 0] = NA
#c.nir[c.nir < 0] = NA
#a.nir[a.nir > 5] = NA
#c.nir[c.nir > 50] = NA
# Mckee et al Opt. Express 2008 method
if (scat.correction == "mckee") {
if (instrument$BB9 ==1) {
# First correct BB9 data using ACS absorption corrected using zaneveld
b.nir = c.nir - a.nir
scat.factor = a.nir / b.nir
scat.factor[scat.factor < 0] = 0
scat.factor.m = matrix(scat.factor, nrow=length(ACS$Depth), ncol=length(ACS$a.wl), byrow=F)
b_star = ACS$c.corrected - ACS$a.TScor.offset
ACS$a.corrected = ACS$a.TScor.offset - (scat.factor.m * b_star)
# correct bb9 data
# Compute Pure Water beta
S.fac = (1+(0.3*BB9$S/37))*1e-4 * (1+ cos(117/180*pi)*cos(117/180*pi)*((1-0.09)/(1+0.09)))
wl.fac = 1.38*(BB9$waves/500)^-4.32
S.fac.m = matrix(S.fac, nrow=length(BB9$Depth), ncol=length(BB9$waves), byrow=F)
wl.fac.m = matrix(wl.fac, nrow=length(BB9$Depth), ncol=length(BB9$waves),byrow=T)
BetaW = S.fac.m * wl.fac.m
bbW = (0.0029308*(BB9$waves/500)^-4.24)/2
bbW.m = matrix(bbW, nrow=length(BB9$Depth), ncol=length(BB9$waves), byrow=T)
Beta2bb <- 2*pi*1.1
ix.BB9.wl = rep(NA,9)
for (j in 1:9){
ix.BB9.wl[j] = which.min(abs(ACS$a.wl-BB9$waves[j]))
}
dt = difftime(max(BB9$Time), min(BB9$Time), units="secs")
# use a 10 seconds time window to define the span parameter
span.BB9 = max(10 / as.numeric(dt), 0.1)
tmp = fit.with.loess(ACS$a.wl[ix.BB9.wl], as.numeric(ACS$Time),
ACS$a.corrected[,ix.BB9.wl], span.BB9, as.numeric(BB9$Time))
a.bb9 = tmp$aop.fitted
# New pathlenght estimated by Doxaran et al 2016
BB9$Beta.corrected <- BB9$Betau * exp(0.01635*a.bb9)
BB9$BetaP.corrected <- BB9$Beta.corrected - BetaW
BB9$bbP.corrected <- Beta2bb * BB9$BetaP.corrected
# Extrapolate bbp to match ACS data matrix.
# Only use the extremity to avoid extrapolation artefact
# A nls fit could be implemented eventually
good.ix <- which(!is.na(BB9$bbP.corrected[,1]))
bbP.acs = t(apply(BB9$bbP.corrected[good.ix,c(2,9)],1,
function(y){
res = spline(BB9$waves[c(2,9)],y,xout=ACS$c.wl, method="natural")$y
}))
# now interpolate for depth
dt = difftime(max(ACS$Time), min(ACS$Time), units="secs")
span.ACS = max(10 / as.numeric(dt), 0.1)
tmp = fit.with.loess(ACS$c.wl, as.numeric(BB9$Time[good.ix]),
bbP.acs, span.ACS, as.numeric(ACS$Time))
bbP.acs = tmp$aop.fitted # overwrite
bbp.tild.0 = bbP.acs / (1.5*b_star) # from Mckee et al. (Fig. 4)
fa0 = apply(bbp.tild.0, 2, function(x){
y = 2.699e-3 + 4.636*x - 3.746e1 *x^2 + 3.177e2*x^3 -1.166e3*x^4})
fc0 = apply(bbp.tild.0, 2, function(x){
y = 6.809e-3/(8.502e-3+x)-(1.918e-2)})
bp1 = b_star / (1 - fa0 - fc0) # from Mckee et al. (eq. 7)
bbp.tild.1 = bbP.acs / bp1 # get a first approximation of bbp tild
for (j in 1:10) { # Replace the condition by a fixed number of iterations.
# This is to speed up the processing
prev.bbp.tild = bbp.tild.1
fa1 = apply(bbp.tild.1, 2, function(x){
y = 2.699e-3 + 4.636*x - 3.746e1 *x^2 + 3.177e2*x^3 -1.166e3*x^4})
fc1 = apply(bbp.tild.1, 2, function(x){
y = 6.809e-3/(8.502e-3+x)-(1.918e-2)})
an = ACS$a.TScor.offset - (fa1*b_star)/(1-fa1-fc1)
cn = ACS$c.corrected + (fc1*b_star)/(1-fa1-fc1)
bp2 = cn - an
bbp.tild.1 = bbP.acs / bp2
}
ACS$a.corrected = an
ACS$c.corrected = cn
ACS$scat.correction.method = "mckee"
} else {
# if (instrument$HS6 ==1) {
# TO BE IMPLEMENTED
# } else {
print("No backscattering measurements to apply Mckee et al method")
print("Automatical change correction to Zaneveld method")
scat.correction = "zaneveld"
# }
}
}
# Zaneveld method
if (scat.correction == "zaneveld") {
b.nir = c.nir - a.nir
scat.factor = a.nir / b.nir
scat.factor.m = matrix(scat.factor, nrow=length(ACS$Depth), ncol=length(ACS$a.wl), byrow=F)
b_star = ACS$c.corrected - ACS$a.TScor.offset
ACS$a.corrected = ACS$a.TScor.offset - (scat.factor.m * b_star)
ACS$scat.correction.method = "zaneveld"
}
# Baseline method
if (scat.correction == "baseline") {
a.nir.m = matrix(a.nir, nrow=length(ACS$Depth), ncol=length(ACS$a.wl), byrow=F)
ACS$a.corrected = ACS$a.TScor.offset - a.nir.m
ACS$scat.correction.method = "baseline"
}
if (scat.correction == "none") {
ACS$a.corrected = ACS$a.TScor.offset
ACS$scat.correction.method = "none"
}
# ACS$a.corrected[ACS$a.corrected < 0 ] = NA
ACS$a.corrected[ACS$a.corrected > 5 ] = NA
ACS$c.corrected[ACS$c.corrected < 0 ] = NA
ACS$c.corrected[ACS$c.corrected > 50 ] = NA
}
####### Apply blank and attenuation corrections to BB9 #####
if (instrument$BB9 ==1) {
if (file.exists(parameters$blank.BB9)) {
if (BB9.DEVICE.FILE) {
print("@@@@@@@@@@@@@@@@@@@@")
print("Device file provided instead of a blank.")
print("In that case, to account for field BB9 blank (or dark measurement)")
print("Please change the offset in the device file...")
print("@@@@@@@@@@@@@@@@@@@@")
BB9$bbP.offset = BB9$bbPu
BB9$bb.offset = BB9$bbu
BB9$BetaP.offset = BB9$BetaPu
BB9$Beta.offset = BB9$Betau
BB9$blank = NA
} else {
print("Apply Blank correction to BB9")
load(file=parameters$blank.BB9)
bb.blank.m = matrix(blank.BB9, nrow=length(BB9$Timer), ncol=length(BB9$waves), byrow=T)
beta.blank.m = matrix(blank.BB9/6.9115, nrow=length(BB9$Timer), ncol=length(BB9$waves), byrow=T)
BB9$bbP.offset = BB9$bbPu - bb.blank.m
BB9$bb.offset = BB9$bbu - bb.blank.m
BB9$BetaP.offset = BB9$BetaPu - beta.blank.m
BB9$Beta.offset = BB9$Betau - beta.blank.m
BB9$blank = blank.BB9
}
} else {
print("No blank correction for BB9")
BB9$bbP.offset = BB9$bbPu
BB9$bb.offset = BB9$bbu
BB9$BetaP.offset = BB9$BetaPu
BB9$Beta.offset = BB9$Betau
BB9$blank = NA
}
############################## BB correction for attenuation
print("Apply attenuation correction to BB9")
# Compute Pure Water beta
S.fac = (1+(0.3*BB9$S/37))*1e-4 * (1+ cos(117/180*pi)*cos(117/180*pi)*((1-0.09)/(1+0.09)))
wl.fac = 1.38*(BB9$waves/500)^-4.32
S.fac.m = matrix(S.fac, nrow=length(BB9$Depth), ncol=length(BB9$waves), byrow=F)
wl.fac.m = matrix(wl.fac, nrow=length(BB9$Depth), ncol=length(BB9$waves),byrow=T)
BetaW = S.fac.m * wl.fac.m
bbW = (0.0029308*(BB9$waves/500)^-4.24)/2
bbW.m = matrix(bbW, nrow=length(BB9$Depth), ncol=length(BB9$waves), byrow=T)
Beta2bb <- 2*pi*1.1
if (instrument$ASPH == 1) {
ix.BB9.wl = rep(NA,9)
for (j in 1:9){
ix.BB9.wl[j] = which.min(abs(ASPH$wl - BB9$waves[j]))
}
dt = difftime(max(BB9$Time), min(BB9$Time), units="secs")
# use a 10 seconds time window to define the span parameter
span.BB9 = max(10 / as.numeric(dt), 0.1)
tmp = fit.with.loess(ASPH$wl[ix.BB9.wl], as.numeric(ASPH$time),
ASPH$a.corrected[,ix.BB9.wl], span.BB9, as.numeric(BB9$Time))
a.bb9 = tmp$aop.fitted
# New pathlenght estimated by Doxaran et al 2016
BB9$Beta.corrected <- BB9$Beta.offset * exp(0.01635*a.bb9) # exp(0.0391*a.bb9)
BB9$BetaP.corrected <- BB9$Beta.corrected - BetaW
BB9$bbP.corrected <- Beta2bb * BB9$BetaP.corrected
BB9$bb.corrected <- BB9$bbP.corrected + bbW.m
BB9$a <- a.bb9
BB9$sigma.correction = TRUE
}
if (instrument$ASPH == 0 & instrument$ACS == 1) {
ix.BB9.wl = rep(NA,9)
for (j in 1:9){
ix.BB9.wl[j] = which.min(abs(ACS$a.wl-BB9$waves[j]))
}
dt = difftime(max(BB9$Time), min(BB9$Time), units="secs")
# use a 10 seconds time window to define the span parameter
span.BB9 = max(max(10 / as.numeric(dt), 0.1), 0.1)
tmp = fit.with.loess(ACS$a.wl[ix.BB9.wl], ACS$Timer,
ACS$a.corrected[,ix.BB9.wl], span.BB9, BB9$Timer)
a.bb9 = tmp$aop.fitted
# clean the data in case of outliers
#a.bb9[a.bb9 > 1] = 0
a.bb9[a.bb9 < 0] = 0
# New pathlenght estimated by Doxaran et al 2016
BB9$Beta.corrected <- BB9$Beta.offset * exp(0.01635*a.bb9)
BB9$BetaP.corrected <- BB9$Beta.corrected - BetaW
BB9$bbP.corrected <- Beta2bb * BB9$BetaP.corrected
BB9$bb.corrected <- BB9$bbP.corrected + bbW.m
BB9$a <- a.bb9
BB9$sigma.correction = TRUE
}
if (instrument$ASPH == 0 & instrument$ACS == 0) {
print("WARNING: No attenuation correction to BB-9 data")
BB9$Beta.corrected <- BB9$Beta.offset
BB9$BetaP.corrected = BB9$Beta.corrected - BetaW
BB9$bbP.corrected = Beta2bb * BB9$BetaP.corrected
BB9$bb.corrected = BB9$bbP.corrected + bbW.m
BB9$a <- NA
BB9$sigma.correction = FALSE
}
}
####### Apply blank and attenuation corrections to BB3 #####
if (instrument$BB3 ==1) {
if (file.exists(parameters$blank.BB3)) {
print("Apply Blank correction to BB3")
load(file=parameters$blank.BB3)
bb.blank.m = matrix(blank.BB3, nrow=length(BB3$Timer), ncol=length(BB3$wl), byrow=T)
beta.blank.m = matrix(blank.BB3/6.9115, nrow=length(BB3$Timer), ncol=length(BB3$wl), byrow=T)
BB3$bbP.offset = BB3$bbP - bb.blank.m
BB3$bb.offset = BB3$bb - bb.blank.m
BB3$BetaP.offset = BB3$BetaP - beta.blank.m
BB3$Beta.offset = BB3$Beta - beta.blank.m
BB3$blank = blank.BB3
} else {
print("No blank correction for BB3")
BB3$bbP.offset = BB3$bbP
BB3$bb.offset = BB3$bb
BB3$BetaP.offset = BB3$BetaP
BB3$Beta.offset = BB3$Beta
BB3$blank = NA
}
############################## BB correction for attenuation
print("Apply attenuation correction to BB3")
# Compute Pure Water beta
S.fac = (1+(0.3*BB3$S/37))*1e-4 * (1+ cos(117/180*pi)*cos(117/180*pi)*((1-0.09)/(1+0.09)))
wl.fac = 1.38*(BB3$wl/500)^-4.32
S.fac.m = matrix(S.fac, nrow=length(BB3$Depth), ncol=length(BB3$wl), byrow=F)
wl.fac.m = matrix(wl.fac, nrow=length(BB3$Depth), ncol=length(BB3$wl),byrow=T)
BetaW = S.fac.m * wl.fac.m
bbW = (0.0029308*(BB3$wl/500)^-4.24)/2
bbW.m = matrix(bbW, nrow=length(BB3$Depth), ncol=length(BB3$wl), byrow=T)
Beta2bb <- 2*pi*1.1
if (instrument$ASPH == 1) {
ix.BB3.wl = rep(NA,3)
for (j in 1:3){
ix.BB3.wl[j] = which.min(abs(ASPH$wl - BB3$wl[j]))
}
dt = difftime(max(BB3$Time), min(BB3$Time), units="secs")
# use a 10 seconds time window to define the span parameter
span.BB3 = max(10 / as.numeric(dt), 0.1)
tmp = fit.with.loess(ASPH$wl[ix.BB3.wl], as.numeric(ASPH$time),
ASPH$a.corrected[,ix.BB3.wl], span.BB3, as.numeric(BB3$Time))
a.BB3 = tmp$aop.fitted
# New pathlenght estimated by Doxaran et al 2016
BB3$Beta.corrected <- BB3$Beta.offset * exp(0.01635*a.BB3) # exp(0.0391*a.BB3)
BB3$BetaP.corrected <- BB3$Beta.corrected - BetaW
BB3$bbP.corrected <- Beta2bb * BB3$BetaP.corrected
BB3$bb.corrected <- BB3$bbP.corrected + bbW.m
BB3$a <- a.BB3
BB3$sigma.correction = TRUE
}
if (instrument$ASPH == 0 & instrument$ACS == 1) {
ix.BB3.wl = rep(NA,3)
for (j in 1:3){
ix.BB3.wl[j] = which.min(abs(ACS$a.wl - BB3$wl[j]))
}
dt = difftime(max(BB3$Time), min(BB3$Time), units="secs")
# use a 10 seconds time window to define the span parameter
span.BB3 = max(10 / as.numeric(dt), 0.1)
tmp = fit.with.loess(ACS$a.wl[ix.BB3.wl], ACS$Timer,
ACS$a.corrected[,ix.BB3.wl], span.BB3, BB3$Timer)
a.BB3 = tmp$aop.fitted
# New pathlenght estimated by Doxaran et al 2016
BB3$Beta.corrected <- BB3$Beta.offset * exp(0.01635*a.BB3)
BB3$BetaP.corrected <- BB3$Beta.corrected - BetaW
BB3$bbP.corrected <- Beta2bb * BB3$BetaP.corrected
BB3$bb.corrected <- BB3$bbP.corrected + bbW.m
BB3$a <- a.BB3
BB3$sigma.correction = TRUE
}
if (instrument$ASPH == 0 & instrument$ACS == 0) {
print("WARNING: No attenuation correction to BB-9 data")
BB3$Beta.corrected <- BB3$Beta.offset
BB3$BetaP.corrected = BB3$Beta.corrected - BetaW
BB3$bbP.corrected = Beta2bb * BB3$BetaP.corrected
BB3$bb.corrected = BB3$bbP.corrected + bbW.m
BB3$a <- NA
BB3$sigma.correction = FALSE
}
}
####### Apply attenuation correction to HS6 data #####
if (instrument$HS6 ==1) {
print("Apply attenuation correction to HS6")
# FROM CALIBRATION FILE
if (HS6.CALYEAR == 2010) Sigmaexp <- c(0.139, 0.142, 0.146, 0.143, 0.146, 0.146)
if (HS6.CALYEAR == 2013) Sigmaexp <- c(0.135, 0.141, 0.143, 0.14, 0.142, 0.144)
if (HS6.CALYEAR == 2014) Sigmaexp <- c(0.137, 0.139, 0.144, 0.141, 0.141, 0.143)
if (HS6.CALYEAR == 2016) Sigmaexp <- c(0.137, 0.141, 0.145, 0.14, 0.144, 0.144)
if (HS6.CALYEAR == 2018) Sigmaexp <- c(0.137, 0.142, 0.144, 0.141, 0.143, 0.145)
Sigmaexp.m = matrix(Sigmaexp, nrow=length(HS6$depth), ncol=length(HS6$wl), byrow=T)
Beta2bb <- 6.79
# Compute Pure Water beta
S.fac = (1+(0.3*HS6$S/37))*1e-4 * (1+ cos(140/180*pi)*cos(140/180*pi)*((1-0.09)/(1+0.09)))
wl.fac = 1.38*(HS6$wl/500)^-4.32
S.fac.m = matrix(S.fac, nrow=length(HS6$depth), ncol=length(HS6$wl), byrow=F)
wl.fac.m = matrix(wl.fac, nrow=length(HS6$depth), ncol=length(HS6$wl),byrow=T)
BetaW = S.fac.m * wl.fac.m
bbW = (0.0029308*(HS6$wl/500)^-4.24)/2
bbW.m = matrix(bbW, nrow=length(HS6$depth), ncol=length(HS6$wl), byrow=T)
ix.HS6.wl = rep(NA,6)
# Match APHERE to HS
if (instrument$ASPH == 1) {
for (j in 1:6){
ix.HS6.wl[j] = which(ASPH$wl == HS6$wl[j])
}
dt = difftime(max(HS6$Time), min(HS6$Time), units="secs")
# use a 10 seconds time window to define the span parameter
span.HS6 = 10 /as.numeric(dt)
tmp = fit.with.loess(ASPH$wl[ix.HS6.wl], as.numeric(ASPH$time),
ASPH$a.corrected[,ix.HS6.wl], span.HS6, as.numeric(HS6$Time))
a.HS6 = tmp$aop.fitted
# New correction from Doxaran et al, Opt. Express 2016
Beta.measured <- HS6$betau
Beta.corrected <- Beta.measured
for (k in 1:5) {
Kbb <- a.HS6 + 4.34 * Beta2bb * (Beta.corrected - BetaW)
sigmaKbb <- exp(Sigmaexp.m*Kbb)
Beta.corrected <- sigmaKbb * Beta.measured
}
HS6$Beta.corrected <- Beta.corrected
HS6$BetaP.corrected <- Beta.corrected - BetaW
HS6$bbP.corrected <- Beta2bb * HS6$BetaP.corrected
HS6$bb.corrected <- HS6$bbP.corrected + bbW.m
HS6$a <- a.HS6
HS6$sigma.correction = TRUE
}
if (instrument$ASPH == 0 & instrument$ACS == 1) {
for (j in 1:6){
ix.HS6.wl[j] = which(ACS$a.wl == HS6$wl[j])
}
dt = difftime(max(HS6$Time), min(HS6$Time), units="secs")
# use a 10 seconds time window to define the span parameter
span.HS6 = 10 /as.numeric(dt)
tmp = fit.with.loess(ACS$a.wl[ix.HS6.wl], as.numeric(ACS$Time),
ACS$a.corrected[,ix.HS6.wl], span.HS6, as.numeric(HS6$Time))
a.HS6 = tmp$aop.fitted
# New correction from Doxaran et al, Opt. Express 2016
Beta.measured <- HS6$betau
Beta.corrected <- Beta.measured
for (k in 1:5) {
Kbb <- a.HS6 + 4.34 * Beta2bb * (Beta.corrected - BetaW)
sigmaKbb <- exp(Sigmaexp.m*Kbb)
Beta.corrected <- sigmaKbb * Beta.measured
}
HS6$Beta.corrected <- Beta.corrected
HS6$BetaP.corrected <- Beta.corrected - BetaW
HS6$bbP.corrected <- Beta2bb * HS6$BetaP.corrected
HS6$bb.corrected <- HS6$bbP.corrected + bbW.m
HS6$a <- a.HS6
HS6$sigma.correction = TRUE
}
if (instrument$ASPH == 0 & instrument$ACS == 0) {
print("WARNING: No attenuation correction to HS-6 data")
HS6$Beta.corrected <- HS6$betau
HS6$BetaP.corrected = HS6$Beta.corrected - BetaW
HS6$bbP.corrected = Beta2bb*HS6$BetaP.corrected
HS6$bb.corrected = HS6$bbP.corrected + bbW.m
HS6$a <- NA
HS6$sigma.correction = FALSE
}
}
####### Compute number of particules from LISST #####
if (instrument$LISST == 1) {
print("compute number of particules")
df = cbind(LISST$lower_size_bins*1e-6,LISST$upper_size_bins*1e-6)
D = exp(apply(log(df),1,mean))
Dm = matrix(D, nrow=length(LISST$Depth), ncol=32, byrow=T)
LISST$N_prime = (LISST$PSD*6)/pi/Dm^3
}
print(paste("Bin the data at a given depth interval:", Zint))
####### Select the portion of the profile to retain for depth interpolation#####
if (is.na(minx) & !is.null(CTD$Time)){
plot(CTD$Time, CTD$Depth)
print("Click for the begining of the cast and then ESC")
minx <- identify(CTD$Time, CTD$Depth)
print("Click for the end of the cast (at surface) and then ESC")
maxx <- identify(CTD$Time, CTD$Depth)
}
if (!is.null(CTD$Time)) {
start.profile = CTD$Time[minx]
stop.profile = CTD$Time[maxx]
CTD$ixmin = minx
CTD$ixmax = maxx
z.max = max(CTD$Depth[minx:maxx])
ix.z.max = which(CTD$Depth == z.max)[1]
CTD$ix.z.max = ix.z.max
} else {
print("No CTD data")
# This case expection was implemented for HS6 profile with no CTD
if (instrument$HS6 == 1) {
if (is.na(minx)) {
plot(HS6$Time, HS6$depth)
print("Click for the begining of the cast and then ESC")
minx <- identify(HS6$Time, HS6$depth)
print("Click for the end of the cast (at surface) and then ESC")
maxx <- identify(HS6$Time, HS6$depth)
}
start.profile = HS6$Time[minx]
stop.profile = HS6$Time[maxx]
HS6$ixmin = minx
HS6$ixmax = maxx
z.max = max(HS6$depth[minx:maxx])
ix.z.max = which(HS6$depth == z.max)[1]
HS6$ix.z.max = ix.z.max
}
}
# Set the span parameter
span.CTD.profile = depth.span/z.max
span.ASPH = (depth.span/z.max)*1.5
span.ACS = depth.span/z.max
span.BB9 = depth.span/z.max
span.BB3 = depth.span/z.max
span.LISST = depth.span/z.max
span.HS6 = depth.span/z.max
span.FLECO = depth.span/z.max
span.FLBBCD = depth.span/z.max
span.FLCHL = depth.span/z.max
time.window=c(start.profile, stop.profile)
depth.fitted=seq(Zint,z.max,Zint)
####### Interpolate CTD data ######
if (!is.null(CTD$Time)) {
print("Bin CTD....")
CTD.fitted.down=list()
df = as.data.frame(cbind(CTD$Depth[minx:ix.z.max], CTD$Temp[minx:ix.z.max]))
names(df) = c("z","T")
mod = loess(T ~ z, data=df, span=span.CTD.profile)
CTD.fitted.down$Temp = predict(mod, depth.fitted)
df = as.data.frame(cbind(CTD$Depth[minx:ix.z.max], CTD$Sal[minx:ix.z.max]))
names(df) = c("z","S")
mod = loess(S ~ z, data=df, span=span.CTD.profile)
CTD.fitted.down$Sal = predict(mod, depth.fitted)
if ((maxx - CTD$ix.z.max) < 10) {
print("No CTD upcast to fit")
CTD.fitted.up=list()
} else {
CTD.fitted.up=list()
df = as.data.frame(cbind(CTD$Depth[maxx:ix.z.max], CTD$Temp[maxx:ix.z.max]))
names(df) = c("z","T")
mod = loess(T ~ z, data=df, span=span.CTD.profile)
CTD.fitted.up$Temp = predict(mod, depth.fitted)
df = as.data.frame(cbind(CTD$Depth[maxx:ix.z.max], CTD$Sal[maxx:ix.z.max]))
names(df) = c("z","S")
mod = loess(S ~ z, data=df, span=span.CTD.profile)
CTD.fitted.up$Sal = predict(mod, depth.fitted)
}
} else {
CTD.fitted.down=list()
CTD.fitted.up=list()
}
####### Interpolate ASPH data#####
if (instrument$ASPH == 1) {
print("Bin ASPH....")
ASPH.fitted.down=list()
ASPH.fitted.up=list()
########## Get ASPH index
# find the index to start and stop profile
abs.diff.time = as.numeric(abs(difftime(ASPH$time, time.window[1])))
ASPH$ixmin = which(abs.diff.time == min(abs.diff.time))
abs.diff.time = as.numeric(abs(difftime(ASPH$time, time.window[2])))
ASPH$ixmax = which(abs.diff.time == min(abs.diff.time))
z.max = max(ASPH$depth[ASPH$ixmin:ASPH$ixmax])
ix.z.max = which(ASPH$depth == z.max)[1]
ASPH$ix.z.max = ix.z.max
ix.good.fitted.depth = which(depth.fitted < z.max)
nz=length(depth.fitted)
depth.raw = ASPH$depth[ASPH$ixmin:ix.z.max]
data.raw = ASPH$a.corrected[ASPH$ixmin:ix.z.max,]
bad.row = which(is.infinite(data.raw), arr.ind = T)[,1]
if (length(bad.row) > 0) tmp=fit.with.loess(ASPH$wl, depth.raw[-bad.row],data.raw[-bad.row,],
span=span.ASPH, depth.fitted[ix.good.fitted.depth])
else tmp=fit.with.loess(ASPH$wl, depth.raw, data.raw,
span=span.ASPH, depth.fitted[ix.good.fitted.depth])
ASPH.fitted.down$a = matrix(nrow=nz, ncol=length(ASPH$wl))
ASPH.fitted.down$a[ix.good.fitted.depth,] = tmp$aop.fitted
# remove extrapolated data for down cast
min.z = min(ASPH$depth[ASPH$ixmin:ix.z.max], na.rm=T)
max.z = max(ASPH$depth[ASPH$ixmin:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) ASPH.fitted.down$a[ix.bad,] = NA
if ((ASPH$ixmax - ix.z.max) < 4) {
print("No ASPH upcast to fit")
} else {
tmp=fit.with.loess(ASPH$wl, ASPH$depth[ASPH$ixmax:ix.z.max], ASPH$a.corrected[ASPH$ixmax:ix.z.max,],
span=span.ASPH, depth.fitted)
ASPH.fitted.up$a = tmp$aop.fitted
# remove extrapolated data for up cast
min.z = min(ASPH$depth[ASPH$ixmax:ix.z.max], na.rm=T)
max.z = max(ASPH$depth[ASPH$ixmax:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) ASPH.fitted.up$a[ix.bad,] = NA
}
} else {
ASPH.fitted.down = list()
ASPH.fitted.up = list()
}
####### Interpolate ACS data#####
if (instrument$ACS == 1) {
print("Bin ACS....")
ACS.fitted.down=list()
ACS.fitted.up=list()
########## Get ACS index
# find the index to start and stop profile
abs.diff.time = as.numeric(abs(difftime(ACS$Time, time.window[1])))
ACS$ixmin = which(abs.diff.time == min(abs.diff.time))
abs.diff.time = as.numeric(abs(difftime(ACS$Time, time.window[2])))
ACS$ixmax = which(abs.diff.time == min(abs.diff.time))
z.max = max(ACS$Depth[ACS$ixmin:ACS$ixmax],na.rm=T)
ix.z.max = which(ACS$Depth == z.max)[1]
ACS$ix.z.max = ix.z.max
# Interpolate using loess function
tmp=fit.with.loess(ACS$a.wl, ACS$Depth[ACS$ixmin:ix.z.max], ACS$a.corrected[ACS$ixmin:ix.z.max,],
span=span.ACS, depth.fitted)
ACS.fitted.down$a = tmp$aop.fitted
tmp=fit.with.loess(ACS$c.wl, ACS$Depth[ACS$ixmin:ix.z.max], ACS$c.corrected[ACS$ixmin:ix.z.max,],
span=span.ACS, depth.fitted)
ACS.fitted.down$c = tmp$aop.fitted
# remove extrapolated data for down cast
min.z = min(ACS$Depth[ACS$ixmin:ix.z.max], na.rm=T)
max.z = max(ACS$Depth[ACS$ixmin:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
ACS.fitted.down$a[ix.bad,] = NA
ACS.fitted.down$c[ix.bad,] = NA
}
if ((ACS$ixmax - ix.z.max) < 5) {
print("No ACS upcast to fit")
} else {
tmp=fit.with.loess(ACS$a.wl, ACS$Depth[ACS$ixmax:ix.z.max], ACS$a.corrected[ACS$ixmax:ix.z.max,],
span=span.ACS, depth.fitted)
ACS.fitted.up$a = tmp$aop.fitted
tmp=fit.with.loess(ACS$c.wl, ACS$Depth[ACS$ixmax:ix.z.max], ACS$c.corrected[ACS$ixmax:ix.z.max,],
span=span.ACS, depth.fitted)
ACS.fitted.up$c = tmp$aop.fitted
# remove extrapolated data for up cast
min.z = min(ACS$Depth[ACS$ixmax:ix.z.max], na.rm=T)
max.z = max(ACS$Depth[ACS$ixmax:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
ACS.fitted.up$a[ix.bad,] = NA
ACS.fitted.up$c[ix.bad,] = NA
}
}
} else {
ACS.fitted.down = list()
ACS.fitted.up = list()
}
####### Interpolate BB9 data#####
if (instrument$BB9 == 1){
print("Bin BB9....")
print("HELLO LE BUGGGG")
BB9.fitted.down=list()
BB9.fitted.up=list()
########## Get BB9 index
# find the index to start and stop profile
abs.diff.time = as.numeric(abs(difftime(BB9$Time, time.window[1])))
BB9$ixmin = which(abs.diff.time == min(abs.diff.time))[1]
abs.diff.time = as.numeric(abs(difftime(BB9$Time, time.window[2])))
BB9$ixmax = which(abs.diff.time == min(abs.diff.time))
z.max = max(BB9$Depth[BB9$ixmin:BB9$ixmax])
ix.z.max = which(BB9$Depth == z.max)[1]
BB9$ix.z.max = ix.z.max
print(BB9$ixmin)
str(BB9)
tmp=fit.with.loess(BB9$waves, BB9$Depth[BB9$ixmin:ix.z.max],
BB9$bbP.corrected[BB9$ixmin:ix.z.max,],
span=span.BB9, depth.fitted)
BB9.fitted.down$bbP = tmp$aop.fitted
tmp=fit.with.loess(BB9$waves, BB9$Depth[BB9$ixmin:ix.z.max],
BB9$bb.corrected[BB9$ixmin:ix.z.max,],
span=span.BB9, depth.fitted)
BB9.fitted.down$bb = tmp$aop.fitted
# remove extrapolated data for down cast
min.z = min(BB9$Depth[BB9$ixmin:ix.z.max], na.rm=T)
max.z = max(BB9$Depth[BB9$ixmin:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
BB9.fitted.down$bb[ix.bad,] = NA
BB9.fitted.down$bbP[ix.bad,] = NA
}
#### Compute bbp spectral slope for down cast
x = 555/BB9$waves
nz=length(depth.fitted)
bbP555=rep(0,nz)
nuP=rep(0,nz)
for (i in 1:nz) {
y = BB9.fitted.down$bbP[i,]
y[y < 0] = NA
if (any(is.na(y))) {
bbP555[i]=NA
nuP[i]=NA
} else {
z = nls(y~b*x^z, start = list(b = y[5], z = 0.5),data=data.frame(x,y), control=list(maxiter=100, warnOnly=T))
bbP555[i]=coef(z)[1]
nuP[i]=coef(z)[2]
}
}
BB9.fitted.down$bbP555 = bbP555
BB9.fitted.down$nuP = nuP
if ((BB9$ixmax - ix.z.max) < 5) {
print("No BB9 upcast to fit")
} else {
tmp=fit.with.loess(BB9$waves, BB9$Depth[BB9$ixmax:ix.z.max],
BB9$bbP.corrected[BB9$ixmax:ix.z.max,],
span=span.BB9, depth.fitted)
BB9.fitted.up$bbP = tmp$aop.fitted
tmp=fit.with.loess(BB9$waves, BB9$Depth[BB9$ixmax:ix.z.max],
BB9$bb.corrected[BB9$ixmax:ix.z.max,],
span=span.BB9, depth.fitted)
BB9.fitted.up$bb = tmp$aop.fitted
# remove extrapolated data for up cast
min.z = min(BB9$Depth[BB9$ixmax:ix.z.max], na.rm=T)
max.z = max(BB9$Depth[BB9$ixmax:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
BB9.fitted.up$bb[ix.bad,] = NA
BB9.fitted.up$bbP[ix.bad,] = NA
}
}
} else {
BB9.fitted.down = list()
BB9.fitted.up = list()
}
####### Interpolate BB3 data#####
if (instrument$BB3 == 1){
print("Bin BB3....")
BB3.fitted.down=list()
BB3.fitted.up=list()
########## Get BB3 index
# find the index to start and stop profile
abs.diff.time = as.numeric(abs(difftime(BB3$Time, time.window[1])))
BB3$ixmin = which(abs.diff.time == min(abs.diff.time))
abs.diff.time = as.numeric(abs(difftime(BB3$Time, time.window[2])))
BB3$ixmax = which(abs.diff.time == min(abs.diff.time))
z.max = max(BB3$Depth[BB3$ixmin:BB3$ixmax])
ix.z.max = which(BB3$Depth == z.max)[1]
BB3$ix.z.max = ix.z.max
tmp=fit.with.loess(BB3$wl, BB3$Depth[BB3$ixmin:ix.z.max],
BB3$bbP.corrected[BB3$ixmin:ix.z.max,],
span=span.BB3, depth.fitted)
BB3.fitted.down$bbP = tmp$aop.fitted
tmp=fit.with.loess(BB3$wl, BB3$Depth[BB3$ixmin:ix.z.max],
BB3$bb.corrected[BB3$ixmin:ix.z.max,],
span=span.BB3, depth.fitted)
BB3.fitted.down$bb = tmp$aop.fitted
# remove extrapolated data for down cast
min.z = min(BB3$Depth[BB3$ixmin:ix.z.max], na.rm=T)
max.z = max(BB3$Depth[BB3$ixmin:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
BB3.fitted.down$bb[ix.bad,] = NA
BB3.fitted.down$bbP[ix.bad,] = NA
}
#### Compute bbp spectral slope for downward cast
x = 555/BB3$wl
nz=length(depth.fitted)
bbP555=rep(0,nz)
nuP=rep(0,nz)
for (i in 1:nz) {
y = BB3.fitted.down$bbP[i,]
y[y < 0] = NA
if (any(is.na(y))) {
bbP555[i]=NA
nuP[i]=NA
} else {
z = nls(y~b*x^z, start = list(b = y[1], z = 0.5),data=data.frame(x,y), control=list(maxiter=100, warnOnly=T))
bbP555[i]=coef(z)[1]
nuP[i]=coef(z)[2]
}
}
BB3.fitted.down$bbP555 = bbP555
BB3.fitted.down$nuP = nuP
if ((BB3$ixmax - ix.z.max) < 5) {
print("No BB3 upcast to fit")
} else {
tmp=fit.with.loess(BB3$wl, BB3$Depth[BB3$ixmax:ix.z.max],
BB3$bbP.corrected[BB3$ixmax:ix.z.max,],
span=span.BB3, depth.fitted)
BB3.fitted.up$bbP = tmp$aop.fitted
tmp=fit.with.loess(BB3$wl, BB3$Depth[BB3$ixmax:ix.z.max],
BB3$bb.corrected[BB3$ixmax:ix.z.max,],
span=span.BB3, depth.fitted)
BB3.fitted.up$bb = tmp$aop.fitted
# remove extrapolated data for up cast
min.z = min(BB3$Depth[BB3$ixmax:ix.z.max], na.rm=T)
max.z = max(BB3$Depth[BB3$ixmax:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
BB3.fitted.up$bb[ix.bad,] = NA
BB3.fitted.up$bbP[ix.bad,] = NA
}
}
} else {
BB3.fitted.down = list()
BB3.fitted.up = list()
}
####### Interpolate LISST data#####
if (instrument$LISST == 1){
print("Bin LISST....")
LISST.fitted.down=list()
LISST.fitted.up=list()
########## Get LISST index
# find the index to start and stop profile
abs.diff.time = as.numeric(abs(difftime(LISST$Time, time.window[1])))
LISST$ixmin = which(abs.diff.time == min(abs.diff.time))
abs.diff.time = as.numeric(abs(difftime(LISST$Time, time.window[2])))
LISST$ixmax = which(abs.diff.time == min(abs.diff.time))
z.max = max(LISST$Depth[LISST$ixmin:LISST$ixmax])
ix.z.max = which(LISST$Depth == z.max)[1]
LISST$ix.z.max = ix.z.max
tmp=fit.with.loess(LISST$median_bins, LISST$Depth[LISST$ixmin:ix.z.max],
LISST$PSD[LISST$ixmin:ix.z.max,],
span=span.LISST, depth.fitted)
LISST.fitted.down$PSD = tmp$aop.fitted
tmp=fit.with.loess(LISST$median_bins, LISST$Depth[LISST$ixmin:ix.z.max],
LISST$N_prime[LISST$ixmin:ix.z.max,],
span=span.LISST, depth.fitted)
LISST.fitted.down$N_prime = tmp$aop.fitted
df = as.data.frame(cbind(LISST$Depth[LISST$ixmin:ix.z.max],
LISST$c670[LISST$ixmin:ix.z.max]))
names(df) = c("z","c670")
mod = loess(c670 ~ z, data=df, span=span.LISST)
LISST.fitted.down$c670 = predict(mod, depth.fitted)
# remove extrapolated data for down cast
min.z = min(LISST$Depth[LISST$ixmin:ix.z.max], na.rm=T)
max.z = max(LISST$Depth[LISST$ixmin:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
LISST.fitted.down$PSD[ix.bad,] = NA
LISST.fitted.down$N_prime[ix.bad,] = NA
LISST.fitted.down$c670[ix.bad] = NA
}
if ((LISST$ixmax - ix.z.max) < 5) {
print("No LISST upcast to fit")
} else {
tmp=fit.with.loess(LISST$median_bins, LISST$Depth[LISST$ixmax:ix.z.max],
LISST$PSD[LISST$ixmax:ix.z.max,],
span=span.LISST, depth.fitted)
LISST.fitted.up$PSD = tmp$aop.fitted
tmp=fit.with.loess(LISST$median_bins, LISST$Depth[LISST$ixmax:ix.z.max],
LISST$N_prime[LISST$ixmax:ix.z.max,],
span=span.LISST, depth.fitted)
LISST.fitted.up$N_prime = tmp$aop.fitted
df = as.data.frame(cbind(LISST$Depth[LISST$ixmax:ix.z.max],
LISST$c670[LISST$ixmax:ix.z.max]))
names(df) = c("z","c670")
mod = loess(c670 ~ z, data=df, span=span.LISST)
LISST.fitted.up$c670 = predict(mod, depth.fitted)
# remove extrapolated data for up cast
min.z = min(LISST$Depth[LISST$ixmax:ix.z.max], na.rm=T)
max.z = max(LISST$Depth[LISST$ixmax:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
LISST.fitted.up$PSD[ix.bad,] = NA
LISST.fitted.up$N_prime[ix.bad,] = NA
LISST.fitted.up$c670[ix.bad] = NA
}
}
} else {
LISST.fitted.down=list()
LISST.fitted.up=list()
}
####### Interpolate HS6 data######
if (instrument$HS6 ==1) {
print("Bin HS6....")
HS6.fitted.down=list()
HS6.fitted.up=list()
########## Get HS6 index
# find the index to start and stop profile
abs.diff.time = as.numeric(abs(difftime(HS6$Time, time.window[1])))
HS6$ixmin = which(abs.diff.time == min(abs.diff.time))[1]
abs.diff.time = as.numeric(abs(difftime(HS6$Time, time.window[2])))
HS6$ixmax = which(abs.diff.time == min(abs.diff.time))[1]
z.max = max(HS6$depth[HS6$ixmin:HS6$ixmax])
ix.z.max = which(HS6$depth == z.max)[1]
HS6$ix.z.max = ix.z.max
# Fit down cast
ix.good.fitted.depth = which(depth.fitted < z.max)
nz=length(depth.fitted)
depth.raw = HS6$depth[HS6$ixmin:ix.z.max]
data.raw = HS6$bbP.corrected[HS6$ixmin:ix.z.max,]
bad.row = which(is.infinite(data.raw), arr.ind = T)[,1]
if (length(bad.row) > 0) tmp=fit.with.loess(HS6$wl, depth.raw[-bad.row],data.raw[-bad.row,],
span=span.HS6, depth.fitted[ix.good.fitted.depth]) else
tmp=fit.with.loess(HS6$wl, depth.raw,data.raw,
span=span.HS6, depth.fitted[ix.good.fitted.depth])
HS6.fitted.down$bbP = matrix(nrow=nz, ncol=6)
HS6.fitted.down$bbP[ix.good.fitted.depth,] = tmp$aop.fitted
data.raw = HS6$bb.corrected[HS6$ixmin:ix.z.max,]
bad.row = which(is.infinite(data.raw), arr.ind = T)[,1]
if (length(bad.row) > 0) tmp=fit.with.loess(HS6$wl, depth.raw[-bad.row],data.raw[-bad.row,],
span=span.HS6, depth.fitted[ix.good.fitted.depth]) else
tmp=fit.with.loess(HS6$wl, depth.raw,data.raw,
span=span.HS6, depth.fitted[ix.good.fitted.depth])
HS6.fitted.down$bb = matrix(nrow=nz, ncol=6)
HS6.fitted.down$bb[ix.good.fitted.depth,] = tmp$aop.fitted
# Fluorescence channels
data.raw = HS6$fluo[HS6$ixmin:ix.z.max,]
bad.row = which(is.infinite(data.raw), arr.ind = T)[,1]
if (length(bad.row) > 0) tmp=fit.with.loess(c(470,700), depth.raw[-bad.row], data.raw[-bad.row,],
span=span.HS6, depth.fitted[ix.good.fitted.depth]) else
tmp=fit.with.loess(c(470,700), depth.raw,data.raw,
span=span.HS6, depth.fitted[ix.good.fitted.depth])
HS6.fitted.down$FDOM = rep(NA,nz)
HS6.fitted.down$FDOM[ix.good.fitted.depth] = tmp$aop.fitted[,1]
HS6.fitted.down$FCHL = rep(NA,nz)
HS6.fitted.down$FCHL[ix.good.fitted.depth] = tmp$aop.fitted[,2]
# remove extrapolated data for down cast
min.z = min(HS6$depth[HS6$ixmin:ix.z.max], na.rm=T)
max.z = max(HS6$depth[HS6$ixmin:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
HS6.fitted.down$FDOM[ix.bad] = NA
HS6.fitted.down$FCHL[ix.bad] = NA
HS6.fitted.down$bbP[ix.bad,] = NA
HS6.fitted.down$bb[ix.bad,] = NA
}
#### Compute bbp and bb spectral slope
x = 555/HS6$wl
nz=length(depth.fitted)
bbP555.down =rep(0,nz)
nuP.down =rep(0,nz)
for (i in 1:nz) {
y = HS6.fitted.down$bbP[i,]
y[y < 0] = NA
y[y > 0.1] = NA
if (any(is.na(y))) {
bbP555.down[i]=NA
nuP.down[i]=NA
} else {
model = nls(y~b*x^z, start = list(b = y[3], z = 1),data=data.frame(x,y), control=list(maxiter=100, warnOnly=T))
bbP555.down[i]=coef(model)[1]
nuP.down[i]=coef(model)[2]
}
}
HS6.fitted.down$bbP555 = bbP555.down
HS6.fitted.down$nuP = nuP.down
if ((HS6$ixmax - ix.z.max) < 5) {
print("No HS6 upcast to fit")
} else {
#### Fit the upcast
depth.raw = HS6$depth[HS6$ixmax:ix.z.max]
data.raw = HS6$bbP.corrected[HS6$ixmax:ix.z.max,]
bad.row = which(is.infinite(data.raw), arr.ind = T)[,1]
if (length(bad.row) > 0){
tmp=fit.with.loess(HS6$wl, depth.raw[-bad.row] ,
data.raw[-bad.row,],
span=span.HS6, depth.fitted)
} else {
tmp=fit.with.loess(HS6$wl, depth.raw,
data.raw,
span=span.HS6, depth.fitted)
}
HS6.fitted.up$bbP = tmp$aop.fitted
data.raw = HS6$bb.corrected[HS6$ixmax:ix.z.max,]
bad.row = which(is.infinite(data.raw), arr.ind = T)[,1]
if (length(bad.row) > 0){
tmp=fit.with.loess(HS6$wl, depth.raw[-bad.row] ,
data.raw[-bad.row,],
span=span.HS6, depth.fitted)
} else {
tmp=fit.with.loess(HS6$wl, depth.raw,
data.raw,
span=span.HS6, depth.fitted)
}
HS6.fitted.up$bb = tmp$aop.fitted
data.raw = HS6$fluo[HS6$ixmax:ix.z.max,]
bad.row = which(is.infinite(data.raw), arr.ind = T)[,1]
if (length(bad.row) > 0){
tmp=fit.with.loess(c(470,700), depth.raw[-bad.row] ,
data.raw[-bad.row,],
span=span.HS6, depth.fitted)
} else {
tmp=fit.with.loess(c(470,700), depth.raw,
data.raw,
span=span.HS6, depth.fitted)
}
HS6.fitted.up$FDOM = tmp$aop.fitted[,1]
HS6.fitted.up$FCHL = tmp$aop.fitted[,2]
# remove extrapolated data for up cast
min.z = min(HS6$depth[HS6$ixmax:ix.z.max], na.rm=T)
max.z = max(HS6$depth[HS6$ixmax:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
HS6.fitted.up$FDOM[ix.bad] = NA
HS6.fitted.up$FCHL[ix.bad] = NA
HS6.fitted.up$bbP[ix.bad,] = NA
HS6.fitted.up$bb[ix.bad,] = NA
}
#### Compute bbp and bb spectral slope
x = 555/HS6$wl
nz=length(depth.fitted)
bbP555.up =rep(0,nz)
nuP.up =rep(0,nz)
for (i in 1:nz) {
y = HS6.fitted.up$bbP[i,]
y[y < 0] = NA
y[y > 0.1] = NA
if (any(is.na(y))) {
bbP555.up[i]=NA
nuP.up[i]=NA
} else {
model = nls(y~b*x^z, start = list(b = y[3], z = 1),data=data.frame(x,y), control=list(maxiter=100, warnOnly=T))
bbP555.up[i]=coef(model)[1]
nuP.up[i]=coef(model)[2]
}
}
HS6.fitted.up$bbP555 = bbP555.up
HS6.fitted.up$nuP = nuP.up
}
} else {
HS6.fitted.down = list()
HS6.fitted.up = list()
}
####### Interpolate FLECO data#####
if (instrument$FLECO ==1) {
print("Bin FLECO....")
FLECO.fitted.down=list()
FLECO.fitted.up=list()
########## Get FLECO index
# find the index to start and stop profile
abs.diff.time = as.numeric(abs(difftime(FLECO$Time, time.window[1])))
FLECO$ixmin = which(abs.diff.time == min(abs.diff.time))
abs.diff.time = as.numeric(abs(difftime(FLECO$Time, time.window[2])))
FLECO$ixmax = which(abs.diff.time == min(abs.diff.time))
z.max = max(FLECO$Depth[FLECO$ixmin:FLECO$ixmax])
ix.z.max = which(FLECO$Depth == z.max)[1]
FLECO$ix.z.max = ix.z.max
# Interpolate using loess function
tmp=fit.with.loess(FLECO$em, FLECO$Depth[FLECO$ixmin:ix.z.max], FLECO$FL[FLECO$ixmin:ix.z.max,],
span=span.FLECO, depth.fitted)
FLECO.fitted.down$FL = tmp$aop.fitted
# remove extrapolated data for down cast
min.z = min(FLECO$depth[FLECO$ixmin:ix.z.max], na.rm=T)
max.z = max(FLECO$depth[FLECO$ixmin:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
FLECO.fitted.down$FL[ix.bad,] = NA
}
if ((FLECO$ixmax - ix.z.max) < 5) {
print("No FLECO upcast to fit")
} else {
tmp=fit.with.loess(FLECO$em, FLECO$Depth[FLECO$ixmax:ix.z.max], FLECO$FL[FLECO$ixmax:ix.z.max,],
span=span.FLECO, depth.fitted)
FLECO.fitted.up$FL = tmp$aop.fitted
# remove extrapolated data for up cast
min.z = min(FLECO$depth[FLECO$ixmax:ix.z.max], na.rm=T)
max.z = max(FLECO$depth[FLECO$ixmax:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
FLECO.fitted.up$FL[ix.bad,] = NA
}
}
} else {
FLECO.fitted.down = list()
FLECO.fitted.up = list()
}
####### Interpolate FLBBCD data #####
if (instrument$FLBBCD ==1) {
print("Bin FLBBCD....")
FLBBCD.fitted.down=list()
FLBBCD.fitted.up=list()
########## Get FLBBCD index
# find the index to start and stop profile
abs.diff.time = as.numeric(abs(difftime(FLBBCD$Time, time.window[1])))
FLBBCD$ixmin = which(abs.diff.time == min(abs.diff.time))
abs.diff.time = as.numeric(abs(difftime(FLBBCD$Time, time.window[2])))
FLBBCD$ixmax = which(abs.diff.time == min(abs.diff.time))
z.max = max(FLBBCD$Depth[FLBBCD$ixmin:FLBBCD$ixmax])
ix.z.max = which(FLBBCD$Depth == z.max)[1]
FLBBCD$ix.z.max = ix.z.max
# Interpolate using loess function down cast
depth.raw = FLBBCD$Depth[FLBBCD$ixmin:ix.z.max]
# Fluorescence channels
data.raw = cbind(FLBBCD$FCHL[FLBBCD$ixmin:ix.z.max],FLBBCD$FDOM[FLBBCD$ixmin:ix.z.max])
bad.row = which(is.infinite(data.raw), arr.ind = T)[,1]
if (length(bad.row) > 0){
tmp=fit.with.loess(c(470,700), depth.raw[-bad.row] ,
data.raw[-bad.row,],
span=span.FLBBCD, depth.fitted)
} else {
tmp=fit.with.loess(c(470,700), depth.raw,
data.raw,
span=span.FLBBCD, depth.fitted)
}
FLBBCD.fitted.down$FCHL = tmp$aop.fitted[,1]
FLBBCD.fitted.down$FDOM = tmp$aop.fitted[,2]
# BB channel
data.raw = cbind(FLBBCD$bb700[FLBBCD$ixmin:ix.z.max],FLBBCD$bbP700[FLBBCD$ixmin:ix.z.max])
bad.row = which(is.infinite(data.raw), arr.ind = T)[,1]
if (length(bad.row) > 0){
tmp=fit.with.loess(c(700,700), depth.raw[-bad.row] ,
data.raw[-bad.row,],
span=span.FLBBCD, depth.fitted)
} else {
tmp=fit.with.loess(c(700,700), depth.raw,
data.raw,
span=span.FLBBCD, depth.fitted)
}
FLBBCD.fitted.down$bb700 = tmp$aop.fitted[,1]
FLBBCD.fitted.down$bbP700 = tmp$aop.fitted[,2]
# remove extrapolated data for down cast
min.z = min(FLBBCD$depth[FLBBCD$ixmin:ix.z.max], na.rm=T)
max.z = max(FLBBCD$depth[FLBBCD$ixmin:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
FLBBCD.fitted.down$FCHL[ix.bad] = NA
FLBBCD.fitted.down$FDOM[ix.bad] = NA
FLBBCD.fitted.down$bb700[ix.bad] = NA
FLBBCD.fitted.down$bbp700[ix.bad] = NA
}
if ((FLBBCD$ixmax - ix.z.max) < 5) {
print("No FLBBCD upcast to fit")
} else {
##### up cast
depth.raw = FLBBCD$Depth[FLBBCD$ixmax:ix.z.max]
# Fluorescence channels
data.raw = cbind(FLBBCD$FCHL[FLBBCD$ixmax:ix.z.max],FLBBCD$FDOM[FLBBCD$ixmax:ix.z.max])
bad.row = which(is.infinite(data.raw), arr.ind = T)[,1]
if (length(bad.row) > 0){
tmp=fit.with.loess(c(470,700), depth.raw[-bad.row] ,
data.raw[-bad.row,],
span=span.FLBBCD, depth.fitted)
} else {
tmp=fit.with.loess(c(470,700), depth.raw,
data.raw,
span=span.FLBBCD, depth.fitted)
}
FLBBCD.fitted.up$FCHL = tmp$aop.fitted[,1]
FLBBCD.fitted.up$FDOM = tmp$aop.fitted[,2]
# BB channel
data.raw = cbind(FLBBCD$bb700[FLBBCD$ixmax:ix.z.max],FLBBCD$bbP700[FLBBCD$ixmax:ix.z.max])
bad.row = which(is.infinite(data.raw), arr.ind = T)[,1]
if (length(bad.row) > 0){
tmp=fit.with.loess(c(700,700), depth.raw[-bad.row] ,
data.raw[-bad.row,],
span=span.FLBBCD, depth.fitted)
} else {
tmp=fit.with.loess(c(700,700), depth.raw,
data.raw,
span=span.FLBBCD, depth.fitted)
}
FLBBCD.fitted.up$bb700 = tmp$aop.fitted[,1]
FLBBCD.fitted.up$bbP700 = tmp$aop.fitted[,2]
# remove extrapolated data for down cast
min.z = min(FLBBCD$depth[FLBBCD$ixmax:ix.z.max], na.rm=T)
max.z = max(FLBBCD$depth[FLBBCD$ixmax:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
FLBBCD.fitted.up$FCHL[ix.bad] = NA
FLBBCD.fitted.up$FDOM[ix.bad] = NA
FLBBCD.fitted.up$bb700[ix.bad] = NA
FLBBCD.fitted.up$bbp700[ix.bad] = NA
}
}
} else {
FLBBCD.fitted.down = list()
FLBBCD.fitted.up = list()
}
####### Interpolate FLCHL data #####
if (instrument$FLCHL ==1) {
print("Bin FLCHL...")
FLCHL.fitted.down=list()
FLCHL.fitted.up=list()
########## Get FLCHL index
# find the index to start and stop profile
abs.diff.time = as.numeric(abs(difftime(FLCHL$Time, time.window[1])))
FLCHL$ixmin = which(abs.diff.time == min(abs.diff.time))
abs.diff.time = as.numeric(abs(difftime(FLCHL$Time, time.window[2])))
FLCHL$ixmax = which(abs.diff.time == min(abs.diff.time))
z.max = max(FLCHL$Depth[FLCHL$ixmin:FLCHL$ixmax])
ix.z.max = which(FLCHL$Depth == z.max)[1]
FLCHL$ix.z.max = ix.z.max
# Interpolate using loess function
df = as.data.frame(cbind(FLCHL$Depth[FLCHL$ixmin:ix.z.max], FLCHL$FCHL[FLCHL$ixmin:ix.z.max]))
names(df) = c("z","FCHL")
mod = loess(FCHL ~ z, data=df, span=span.FLCHL)
FLCHL.fitted.down$FCHL = predict(mod, depth.fitted)
# remove extrapolated data for down cast
min.z = min(FLCHL$Depth[FLCHL$ixmin:ix.z.max], na.rm=T)
max.z = max(FLCHL$Depth[FLCHL$ixmin:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
FLCHL.fitted.down$FCHL[ix.bad] = NA
}
if ((FLCHL$ixmax - ix.z.max) < 5) {
print("No FLCHL upcast to fit")
} else {
# Interpolate using loess function
df = as.data.frame(cbind(FLCHL$Depth[FLCHL$ixmax:ix.z.max], FLCHL$FCHL[FLCHL$ixmax:ix.z.max]))
names(df) = c("z","FCHL")
mod = loess(FCHL ~ z, data=df, span=span.FLCHL)
FLCHL.fitted.up$FCHL = predict(mod, depth.fitted)
# remove extrapolated data for up cast
min.z = min(FLCHL$Depth[FLCHL$ixmax:ix.z.max], na.rm=T)
max.z = max(FLCHL$Depth[FLCHL$ixmax:ix.z.max], na.rm=T)
ix.bad = which(depth.fitted < min.z | depth.fitted > max.z)
if (length(ix.bad) > 0) {
FLCHL.fitted.up$FCHL[ix.bad] = NA
}
}
} else {
FLCHL.fitted.down = list()
FLCHL.fitted.up = list()
}
####### PUT TOGETHER ALL PARAMETERS######
IOP = list(time.window=time.window, Time0.CTD = Time0.CTD,
ASPH=ASPH, ACS=ACS, BB9=BB9, CTD=CTD, LISST=LISST,
HS6=HS6, FLECO=FLECO, BB3=BB3, FLBBCD=FLBBCD, FLCHL=FLCHL)
IOP.fitted.down = list(Depth = depth.fitted, ASPH=ASPH.fitted.down,
ACS=ACS.fitted.down,
BB9=BB9.fitted.down,
CTD=CTD.fitted.down,
LISST=LISST.fitted.down,
HS6=HS6.fitted.down,
FLECO=FLECO.fitted.down,
BB3=BB3.fitted.down,
FLBBCD=FLBBCD.fitted.down,
FLCHL=FLCHL.fitted.down)
IOP.fitted.up = list(Depth = depth.fitted, ASPH=ASPH.fitted.up,
ACS=ACS.fitted.up,
BB9=BB9.fitted.up,
CTD=CTD.fitted.up,
LISST=LISST.fitted.up,
HS6=HS6.fitted.up,
FLECO=FLECO.fitted.up,
BB3=BB3.fitted.up,
FLBBCD=FLBBCD.fitted.up,
FLCHL=FLCHL.fitted.up)
print("Save data into RData format: IOP, IOP.fitted.down,IOP.fitted.up")
save(IOP,file="IOP.RData")
save(IOP.fitted.down,file="IOP.fitted.down.RData")
save(IOP.fitted.up,file="IOP.fitted.up.RData")
# update cast.info.file
print("UPDATING FILE: cast.info.dat")
df = data.frame(
parameters$lon,
parameters$lat,
parameters$cast,
Time0.CTD,
Time0.LISST,
minx,
maxx,
parameters$Zint,
parameters$asph.skip,
parameters$maxbb,
parameters$Ndepth.to.plot,
parameters$depth.interval.for.smoothing)
names(df) <- c( "lon", "lat", "cast", "Time0.CTD", "Time0.LISST",
"minx", "maxx", "Zint", "asph.skip",
"maxbb", "Ndepth.to.plot","depth.interval.for.smoothing")
write.table(df, file="cast.info.dat", sep=",", quote=F, row.names = F)
print("Return IOP data object")
return(IOP)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.