R/compute.aTOT.for.COPS.R

In belasi01/Riops: IOPs data processing

#' 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)
}