shadow.correction <- function(instr, cops, SB=NA) {
mymessage("shadow correction")
chl <- cops$chl
SHADOW.CORRECTION <- cops$SHADOW.CORRECTION
SHADOW.CORRECTION.FROM.KD <- NA
absorption.waves <- cops$absorption.waves
absorption.values <- cops$absorption.values
date.mean <- cops$date.mean
sunzen <- cops$sunzen
longitude <- cops$longitude
latitude <- cops$latitude
waves <- cops[[paste(instr, "waves", sep = ".")]]
julian.day <- as.numeric(format(date.mean, format = "%j"))
if(!is.na(chl) && chl != 999) { #### Use Morel and Maritorena model to
#### get absorption from CHL, but if 999 then use the Kd derived absorption
# 350-700
chlTMP <- chl
chlTMP <- min(chlTMP, 9.99999)
chlTMP <- max(chlTMP, 0.03001)
wavesTMP <- waves
wavesTMP[wavesTMP > 349.999 & wavesTMP < 350.001] <- 350.001
wavesTMP[wavesTMP > 699.999 & wavesTMP < 700.001] <- 699.999
i.wavesSUP <- wavesTMP < 350 | wavesTMP > 700
wavesTMP[i.wavesSUP] <- 550
aR <- as.vector(popt.fv.a(wavesTMP, chlTMP)) * radius.instrument.optics[instr]
aR[i.wavesSUP] <- NA
# > 700
wavesTMP <- waves
wavesTMP[wavesTMP > 699.999 & wavesTMP < 700.001] <- 699.999
wavesTMP[wavesTMP > 864.999] <- 864.999
i.wavesIR <- wavesTMP > 700
aR[i.wavesIR] <- popt.f.aw(wavesTMP[i.wavesIR]) * radius.instrument.optics[instr]
# < 350
wavesTMP <- waves
wavesTMP[wavesTMP > 349.999 & wavesTMP < 350.001] <- 350.001
i.wavesUV <- wavesTMP < 350
aR[i.wavesUV] <- popt.f.a(350.001, chlTMP) * radius.instrument.optics[instr]
} else {
if (!is.na(chl) && chl == 999) { #### CHL == 999 mean used Kd-derived absorption
##### This correction needs an estimation of Kd and R.
##### It uses linear fit if valid, and loess estimation
##### if any linear fit failed (for either Lu and Ed!)
##### Begin with LuZ
if (!is.null(cops$LuZ.0m.linear)) {
Q = 4 # this is an approximation of the Q-Factor for the
# calulation of the sub-surface irradiance reflectance
#### Check for a valid Kd with linear interpolation
ix.LuZ.0m.valid = which(!is.na(cops$LuZ.0m))
if (anyNA(cops$K.EdZ.surf[ix.LuZ.0m.valid]) |
anyNA(cops$K.LuZ.surf[ix.LuZ.0m.valid])) {
# check for the depth used for the LuZ extrapolation
#max.depth <- max(cops$LuZ.Z.interval, na.rm = T)
#ix.max.depth <- which.min(abs(cops$depth.fitted - max.depth))
ix.2.5 <- which.min(abs(cops$depth.fitted - 2.5))
Kd = cops$K0.EdZ.fitted[ix.2.5,]
Ed.0m = cops$EdZ.0m
R = (cops$LuZ.0m*Q)/Ed.0m
SHADOW.CORRECTION.FROM.KD <- "LOESS"
} else {
Kd = cops$K.EdZ.surf
Ed.0m = cops$EdZ.0m.linear
R = (cops$LuZ.0m.linear*Q)/Ed.0m
SHADOW.CORRECTION.FROM.KD <- "LINEAR"
}
# I think this condition is not necessary any more...
#if (length(ix.LuZ.0m.valid) == 0) {
### This is when SHADOW.CORRECTION == 999
### but the linear fit failed at all bands....
### The Kd for the to 2.5 meter is assumed.
# ix.2.5 <- which.min(abs(cops$depth.fitted - 2.5))
# Kd = cops$K0.EdZ.fitted[ix.2.5,]
# Ed.0m = cops$EdZ.0m
# R = (cops$LuZ.0m*Q)/Ed.0m
# SHADOW.CORRECTION.FROM.KD <- "LOESS"
#}
}
##### The with EuZ
if (!is.null(cops$EuZ.0m.linear)) {
#### Check for a valid Kd with linear interpolation
ix.EuZ.0m.valid = which(!is.na(cops$EuZ.0m))
if (anyNA(cops$K.EdZ.surf[ix.EuZ.0m.valid]) |
anyNA(cops$K.EuZ.surf[ix.LuZ.0m.valid])) {
# check for the depth used for the LuZ extrapolation
ix.2.5 <- which.min(abs(cops$depth.fitted - 2.5))
Kd = cops$K0.EdZ.fitted[ix.2.5,]
Ed.0m = cops$EdZ.0m
R = (cops$EuZ.0m)/Ed.0m
SHADOW.CORRECTION.FROM.KD <- "LOESS"
} else {
Kd = cops$K.EdZ.surf
Ed.0m = cops$EdZ.0m.linear
R = cops$EuZ.0m.linear/Ed.0m
SHADOW.CORRECTION.FROM.KD <- "LINEAR"
}
if (length(ix.EuZ.0m.valid) == 0) {
### This is when SHADOW.CORRECTION == 999
### but the linear fit failed at all bands....
### The Kd for the to 2.5 meter is assumed.
ix.2.5 <- which.min(abs(cops$depth.fitted - 2.5))
Kd = cops$K0.EdZ.fitted[ix.2.5,]
Ed.0m = cops$EdZ.0m
R = (cops$EuZ.0m)/Ed.0m
SHADOW.CORRECTION.FROM.KD <- "LOESS"
}
}
# From Morel and Maritorena JGR 2001 eq 8'
absorption.values <- Kd*0.9*(1-R)/(1+2.25*R)
cops$absorption.values <- absorption.values
aR <- absorption.values * radius.instrument.optics[instr]
} else {
agree.waves <- isTRUE(all.equal(absorption.waves, waves, scale = 1, tolerance = 1))
if(!agree.waves) {
cat("wavelengths in file absorption.cops.dat do not agree with wavelengths of", instr, "instrument\n")
cat(instr, ":", waves, "\n")
cat("file absorption.cops.dat", ":", absorption.waves, "\n")
stop()
}
aR <- absorption.values * radius.instrument.optics[instr]
}
}
names(aR) <- waves
### Add the case when the Shadow Band data is available (Simon Belanger 2018/12/15)
if (is.na(SB)[1]) {
# This code was added to improve the Edif and Edir estimation using Gregg and Carder model
# The visibility is reduced to match the measured Ed0.0p obtained from the COPS reference sensor
print ("No Bioshade measurements")
print ("Compute Ed0.dif/Ed0.tot using Gregg and Carder model")
julian.day <- as.numeric(format(cops$date.mean, format = "%j"))
visibility <- 25
egc <- GreggCarder.f(julian.day, longitude, latitude, sunzen, lam.sel = waves, Vi=visibility)
ix.490 <- which.min(abs(waves - 490))
ratio = egc$Ed[ix.490]*100/cops$Ed0.0p[ix.490]
while (ratio > 1.05 & visibility > 0.5) {
egc <- GreggCarder.f(julian.day, longitude, latitude, sunzen,lam.sel = waves, Vi=visibility)
ratio = egc$Ed[ix.490]*100/cops$Ed0.0p[ix.490]
visibility = visibility - 0.5
}
Edif <- egc$Edif*100 # factor 100 to convert into same COPS units
Edir <- egc$Edir*100
} else {
Edif <- SB$Ed0.dif
Edir <- SB$Ed0.tot - Edif
}
###
ratio.edsky.edsun <- Edif / Edir
epss <- shadow.epsilon(instr, aR, sunzen, ratio.edsky.edsun)
eps.sun <- epss$eps.sun
eps.sky <- epss$eps.sky
eps <- epss$eps
correction <- 1 - eps
corr.names <- c("shad.aR", "shad.Edif", "shad.Edir", "shad.ratio.edsky.edsun",
"shad.eps.sun", "shad.eps.sky", "shad.eps", "shad.correction")
corr.names <- paste(instr, corr.names, sep = ".")
corr.names <- c(corr.names,"absorption.values", "absorption.waves", "SHADOW.CORRECTION.FROM.KD")
corr.list <- list(aR, Edif, Edir, ratio.edsky.edsun, eps.sun,
eps.sky, eps, correction, absorption.values, waves, SHADOW.CORRECTION.FROM.KD)
names(corr.list) <- corr.names
corr.list
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.