#' Compute the total absorption for COPS processing of LuZ profile
#'
#' Compute the total absorption for COPS instrument self-shadow
#' correction of LuZ. It is run only after a preliminary processing
#' of the COPS data using the \pkg{Cops} package.
#'
#' @param pathIOP is the path where IOPs data are located
#' (i.e. IOPs.fitted.down.RData or IOPs.fitted.up.RData). These files
#' are generated by \code{\link{correct.merge.IOP.profile}}.
#' IMPORTANT NOTE: a COPS folder with files preliminarly processed using
#' the \pkg{Cops} package must be present in the parent folder "YYYYMMDD_StationXXX".
#' The \pkg{Cops} processing creates a file named absorption.cops.dat
#' that will be automatically edited.
#' @param cast is either "down" or "up" indicating whether the down cast
#' or the up cast will be used to average the non-water absorption.
#' The default is "down".
#' @param depth.interval is the depth interval for which the non-water
#' absorption will be avaraged.
#' The default is c(0.75,2.1)
#' @param a.instrument indicates whether the non-water absorption will be
#' taken from the a-sphere (ASPH) or the ac-s (ACS).
#' The default is "ASPH".
#'
#' @details
#' The program will open a RData file located in ../COPS/BIN folder to get
#' the wavelengths available on the COPS. It will edit the file
#' ../COPS/absorption.cops.dat and produce a PNG figure showing the spectral
#' absorption that will be used in the COPS processing for instrument self-shadow
#' correction.
#'
#' The pure water absorption is taken from a composite table
#' (Kou etal 1993, Pope and Fry 1997, etc).
#'
#' @return It returns a vector of total absorption.
#'
#' @seealso \code{\link{correct.merge.IOP.profile}} and
#' \code{\link{IOPs.go}}
#'
#' @author Simon Belanger
#' @export
compute.aTOT.for.COPS <- function(pathIOP=".", cast="down", depth.interval=c(0.75,2.1), instrument = "ASPH") {
setwd(pathIOP)
# Reading IOP files
if (cast == "down" ) {
load(file="IOP.fitted.down.RData")
IOP.fitted = IOP.fitted.down
} else {
load(file="IOP.fitted.up.RData")
IOP.fitted = IOP.fitted.up
}
load(file="IOP.RData") # this is to get the wavelenghts
# Finding the COPS files avaiblable
file.cops = list.files(path="../COPS/BIN/", pattern = "*\\.RData")
# Reading the first one to get the cops wavelength
if (length(file.cops) > 1) {
load(file=paste("../COPS/BIN/",file.cops[1],sep=""))
lambda = cops$Ed0.waves
} else {
print("NO COPS FILES PROCESSED FOR THIS STATION")
return(0)
}
# get Pure water absorption at COPS wavelengths
a.w.cops = spectral.aw(lambda)
# get non-water absorption
ix = which(IOP.fitted$Depth > depth.interval[1] & IOP.fitted$Depth < depth.interval[2])
if (instrument == "ASPH") {
if (!is.null(IOP.fitted$ASPH$a)) {
a.tot.w = apply(IOP.fitted$ASPH$a[ix,],2,mean, na.rm=T)
a.tot.w.cops = spline(IOP$ASPH$wl, a.tot.w, xout=lambda, method="natural")$y
# for lambda < 360:
ix.cops.UV = which(lambda < 360)
ix.wl365.asph = which(IOP$ASPH$wl == 365)
ix.wl400.asph = which(IOP$ASPH$wl == 400)
df = as.data.frame(cbind(IOP$ASPH$wl[ix.wl365.asph:ix.wl400.asph], a.tot.w[ix.wl365.asph:ix.wl400.asph]))
names(df) <- c("wl", "a")
model <- nls(a~a400*exp(-S*(wl-400)), start=list(a400=a.tot.w[ix.wl400.asph], S=0.02),data=df)
a.tot.w.cops[ix.cops.UV] = predict(model, newdata=list(wl=lambda[ix.cops.UV]))
# for NIR nm
# for lambda > 750:
ix.cops.NIR = which(lambda > 750)
ix.wl760.asph = which(IOP$ASPH$wl == 760)
df = as.data.frame(cbind(IOP$ASPH$wl[ix.wl365.asph:ix.wl760.asph], a.tot.w[ix.wl365.asph:ix.wl760.asph]))
names(df) <- c("wl", "a")
model <- nls(a~a400*exp(-S*(wl-400)), start=list(a400=a.tot.w[ix.wl400.asph], S=0.02),data=df)
a.tot.w.cops[ix.cops.NIR] = predict(model, newdata=list(wl=lambda[ix.cops.NIR]))
} else {
print("NO ASPH AVAILABLE IN IOPs")
print("IF ACS IS AVAILABLE, SET instrument PARAMETER TO ACS")
print("ABORT PROCESSING")
return(0)
}
}
if (instrument == "ACS") {
if (!is.null(IOP.fitted$ACS$a)) {
a.tot.w = apply(IOP.fitted$ACS$a[ix,],2,mean, na.rm=T)
a.tot.w.cops = spline(IOP$ACS$a.wl, a.tot.w, xout=lambda, method="natural")$y
# for lambda < min lambda available:
ix.cops.UV = which(lambda < IOP$ACS$a.wl[1])
ix.wl420.acs = which.min(abs(IOP$ACS$a.wl-420))
ix.wl400.acs = which.min(abs(IOP$ACS$a.wl-400))
df = as.data.frame(cbind(IOP$ACS$a.wl[ix.wl400.acs:ix.wl420.acs], a.tot.w[ix.wl400.acs:ix.wl420.acs]))
names(df) <- c("wl", "a")
model <- nls(a~a400*exp(-S*(wl-400)), start=list(a400=a.tot.w[ix.wl400.acs], S=0.02),data=df)
a.tot.w.cops[ix.cops.UV] = predict(model, newdata=list(wl=lambda[ix.cops.UV]))
# for NIR nm assume 0
ix.cops.NIR = which(lambda > max(IOP$ACS$a.wl))
a.tot.w.cops[ix.cops.NIR] = 0
} else {
print("NO ACS AVAILABLE IN IOPs")
print("IF ASPH IS AVAILABLE, SET instrument PARAMETER TO ASPH")
print("ABORT PROCESSING")
return(0)
}
}
a.tot.w.cops[a.tot.w.cops < 0] = 0
a.tot.cops = a.tot.w.cops + a.w.cops
# Plot and save results
png(filename = "absorption.for.cops.png")
plot(lambda, a.tot.cops,
xlab = "Wavelenght", ylab="Absorption", pch=19,
ylim=c(0, max(a.tot.cops)))
if (instrument == "ASPH") lines(IOP$ASPH$wl, a.tot.w, col=2, lwd=3)
if (instrument == "ACS") lines(IOP$ACS$a.wl, a.tot.w, col=2, lwd=3)
lines(300:900, spectral.aw(300:900), col=4, lwd=4)
dev.off()
# write results in file absorption.cops.dat
df = read.table(file="../COPS/absorption.cops.dat", sep=";")
nfile = length(df$V1) - 1
a.mat = matrix(a.tot.cops, nrow=nfile, ncol=length(lambda), byrow=T)
df.a = rbind(lambda,a.mat)
df.final = as.data.frame(cbind(df$V1, as.data.frame(df.a)))
write.table(df.final, quote=F, file="../COPS/absorption.cops.dat", sep=";", col.names = F, row.names=F)
return(a.tot.cops)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.