process.EdZ <- function(cops.raw,
cops.dd,
cops.black,
cops.Ed0,
EXTRAPOLATION.0m = TRUE) {
mymessage(paste("Processing EdZ", " ..."), head = "+")
aop <- cops.raw$EdZ
correction <- cops.Ed0$Ed0.correction
aop <- aop * correction
waves <- as.numeric(cops.raw$EdZ.waves)
Depth <- cops.dd$Depth + delta.capteur.optics["EdZ"]
depth.fitted <- cops.dd$depth.fitted
if(!is.null(cops.dd[["EdZ.tilt"]])) {
tilt <- cops.dd[["EdZ.tilt"]]
} else {
tilt <- cops.dd[["EuZ.tilt"]]
}
if(!is.null(cops.black)) {
black <- cops.black$EdZ
aop <- t(t(aop) - black)
}
# tilt limitation
valid.tilt <- tilt < tiltmax.optics["EdZ"]
aop.all <- aop[cops.dd$Depth.good & valid.tilt,] # make a backup
Depth.all <- Depth[cops.dd$Depth.good & valid.tilt]
Depth.kept <- cops.dd$Depth.good & valid.tilt & Depth > sub.surface.removed.layer.optics["EdZ"]
Depth <- Depth[Depth.kept]
tilt <- tilt[Depth.kept]
aop <- aop[Depth.kept, ]
##### Remove radiometric measurement below the detection limit
aop.all[aop.all < 0] <-0
for (w in 1:length(waves)) {
aop[aop[,w] <= EdZ.detect.lim[w],w] <- NA
}
# METHOD loess to fit this aop
if (EXTRAPOLATION.0m) {
idx.depth.0 <- which.min(abs(depth.fitted - 0))
# idx.depth.0 <- 1 ### Il me semble que ca donne la même chose...
} else {
idx.depth.0 <- which.min(abs(depth.fitted - Depth[1])) # ice version
}
depth.0 <- depth.fitted[idx.depth.0]
span <- depth.interval.for.smoothing.optics["EdZ"]
fitted <- fit.with.loess(waves, Depth, log(aop), span, depth.fitted,
idx.depth.0,
span.wave.correction = TRUE, DEPTH.SPAN = TRUE ,
minimum.obs = 10)
aop.fitted <- exp(fitted$aop.fitted)
aop.0 <- exp(fitted$aop.0)
K <- compute.K(depth.fitted, idx.depth.0, aop.0, aop.fitted)
KZ.fitted <- K$KZ.fitted
K0.fitted <- K$K0.fitted
#### Clean calculated AOP
n.fitted <- length(depth.fitted)
for (w in 1:length(waves)) {
print(paste("Clean AOP for EdZ ", waves[w]))
if (!all(is.na(aop.fitted[,w]))) { # if all NA, then the AOPs allready equal to NA
# Apply a spline on raw data for further flaging on the AOP
tmp = smooth.spline(Depth.all, aop.all[,w], spar=0.2) # spar was added to make the data smoother
aop.spline = spline(tmp, xout = depth.fitted, method = 'natural')$y
# remove bad data
KZ.fitted[(aop.fitted[(2:n.fitted),w] <= EdZ.detect.lim[w] |
aop.spline[2:n.fitted] <= EdZ.detect.lim[w] &
depth.fitted[2:n.fitted] > Depth[1]),w] <- NA
K0.fitted[(aop.fitted[(2:n.fitted),w] <= EdZ.detect.lim[w] |
aop.spline[2:n.fitted] <= EdZ.detect.lim[w]&
depth.fitted[2:n.fitted] > Depth[1]),w] <- NA
aop.fitted[(aop.fitted[,w] <= EdZ.detect.lim[w] |
aop.spline <= EdZ.detect.lim[w]&
depth.fitted > Depth[1]),w] <- NA
}
if (!is.na(aop.0[w])) {
if (aop.0[w] <= EdZ.detect.lim[w]) aop.0[w] <- NA
}
}
# Extrapolate EdZ to 0- using linear method
K <- compute.Ksurf.linear(Depth, aop,
instrument = "EdZ",
delta.depth= linear.fit.max.delta.depth.optics["EdZ"],
r2.threshold=linear.fit.Rsquared.threshold.optics["EdZ"])
r2 <- K$r2
K.surf <- K$Kx
Z.interval <- K$Z.interval
ix.Z.interval <- K$ix.Z.interval
if (EXTRAPOLATION.0m) {
EdZ.0m.linear <- K$X.0m
} else {
EdZ.0m.linear <- K$X.0m * exp(-depth.0*K.surf)
}
KolmolSmirnov.p.value<-K$KolmolSmirnov.p.value
EdZ.fitted <- matrix(NA, ncol=dim(aop.fitted)[2], nrow=dim(aop.fitted)[1])
for (i in 1:length(waves)) {
if (!is.na(EdZ.0m.linear[i])) {
EdZ.fitted[idx.depth.0:n.fitted,i] <- K$X.0m[i] * exp(-depth.fitted[idx.depth.0:n.fitted]*K.surf[i])
}
}
# Compute PAR at depth
PAR <- compute.PARz(depth.fitted, waves, aop.fitted, aop.0,
f.PAR=c(0.001, 0.01, 0.05, 0.1,0.5))
PLOT.LINEAR <- !all(is.na(EdZ.0m.linear))
if (!PLOT.LINEAR) Z.interval=5
# PLOT
aop.cols <- rainbow.modified(length(waves))
if(INTERACTIVE) x11(width = win.width, height = win.height)
matplot(aop.all, Depth.all, type = "p", log = "x",
ylim = c(max(Z.interval,na.rm = T)+1,0),
#ylim = rev(range(c(0, Depth))),
xlim=c(min(EdZ.detect.lim),max(aop, na.rm=T)), pch = ".", cex = 1,
ylab="Depth (m)",
xlab = expression(E[d]*z* ~~ "("*mu*W.*cm^{-2}*.nm^{-1}*")"),
col = aop.cols)
grid(col = 1)
matplot(aop.fitted, depth.fitted, type = "l",
lty = 1, lwd = 2, add = TRUE, col = aop.cols)
points(aop.0, rep(depth.0, length(aop.0)), col = aop.cols)
if (PLOT.LINEAR) matplot(EdZ.fitted, depth.fitted, type = "l",
lty = 2, lwd = 2, add = TRUE, col = aop.cols)
if (PLOT.LINEAR) points(diag(aop[ix.Z.interval,]), Z.interval, cex=1.5,
pch = 19,col = aop.cols)
par(xpd = TRUE)
legend(10^par("usr")[1], par("usr")[4], legend = waves, xjust = 0, yjust = 0, lty = 1, lwd = 2, col = aop.cols, ncol = ceiling(length(waves) / 2), cex = 0.75)
par(xpd = FALSE)
if (!PLOT.LINEAR) text(aop.0[8], 0, "LINEAR INTERPOLATION FAILED, SHOULD YOU RELAX THE TILT THRESOLD?", pos=4)
if(INTERACTIVE) x11(width = win.width, height = win.height)
mai.old1 <- par("mai")
mgp.old1 <- par("mgp")
par(mfrow = c(4, 5))
par(mai = par("mai") / 2, mgp = par("mgp") / 2)
mai.old2 <- par("mai")
mgp.old2 <- par("mgp")
for(i in 1:length(waves)) {
if (length(which(!is.na(aop[,i]))) > 0) {
plot(aop.all[, i], Depth.all, type = "p", log = "x",
xlim = range(aop[aop[, i] > 0, i], aop.0[i], na.rm = TRUE),
ylim = rev(range(Depth, depth.fitted)),
pch = ".",
xlab = "", ylab = "", axes = FALSE, frame.plot = TRUE,
main = substitute(E[d]*z~x~"("*mu*W.*cm^{-2}*.nm^{-1}*")",list(x = waves[i])))
grid(col = 1)
lines(aop.fitted[, i], depth.fitted, col = 2)
points(aop.0[i], depth.0, pch = 20, col = 4)
abline(v=EdZ.detect.lim[i], col="orange", lwd=3)
abline(h=sub.surface.removed.layer.optics["EdZ"], col="green", lwd=3)
axis(1)
axis(2)
}
}
par(mai = mai.old2)
par(mgp = mgp.old2)
par(mfrow = c(1, 1))
par(mai = mai.old1)
par(mgp = mgp.old1)
##### check surface extrapolation
if(INTERACTIVE) x11(width = win.width, height = win.height)
mai.old1 <- par("mai")
mgp.old1 <- par("mgp")
par(mfrow = c(4, 5))
par(mai = par("mai") / 2, mgp = par("mgp") / 2)
mai.old2 <- par("mai")
mgp.old2 <- par("mgp")
for(i in 1:length(waves)) {
if (length(which(!is.na(aop[,i]))) > 0) {
if (PLOT.LINEAR) {
my.xlim = range(aop[, i], aop.0[i], EdZ.0m.linear[i], na.rm = TRUE)
} else {
my.xlim = range(aop[, i], aop.0[i], na.rm = TRUE)
}
plot(aop.all[, i], Depth.all, type = "p", log = "x",
xlim = my.xlim,
ylim = c(max(Z.interval,na.rm = T)+0.5,0),
pch = ".", xlab = "", ylab = "",
axes = FALSE, frame.plot = TRUE,
main = substitute(E[d]*z~x~r^2==r.2,list(x = waves[i], r.2=signif(r2[i],3))))
grid(col = 1)
lines(aop.fitted[, i], depth.fitted, col = "red")
points(aop.0[i], depth.0, pch = 1, cex=1.5, col = "red")
if (PLOT.LINEAR) lines(EdZ.fitted[, i], depth.fitted, col = "blue")
if (PLOT.LINEAR) points(EdZ.0m.linear[i],depth.0, pch = 1, cex=1.5, col = "blue")
if (PLOT.LINEAR) points(aop[ix.Z.interval[i],i], Z.interval[i], cex=1.5,
pch = 19,col ="blue")
abline(v=EdZ.detect.lim[i], col="orange", lwd=2)
abline(h=sub.surface.removed.layer.optics["EdZ"], col="green", lwd=1)
axis(1)
axis(2)
}
}
par(mai = mai.old2)
par(mgp = mgp.old2)
par(mfrow = c(1, 1))
par(mai = mai.old1)
par(mgp = mgp.old1)
##### Same as above but 4 selected wavelenghts
if(INTERACTIVE) x11(width = win.width, height = win.height)
mai.old1 <- par("mai")
mgp.old1 <- par("mgp")
par(mfrow = c(2, 2))
par(mai = par("mai") / 2, mgp = par("mgp") / 2)
mai.old2 <- par("mai")
mgp.old2 <- par("mgp")
#### Select 4 wavelengths
if (PLOT.LINEAR) {
ix.w = which(!is.na(r2))
} else ix.w <- 1:19
for(i in floor(seq(ix.w[1], ix.w[length(ix.w)], length.out = 4))) {
if (length(which(!is.na(aop[,i]))) > 0) {
if (PLOT.LINEAR) {
my.xlim = range(aop[, i], aop.0[i], EdZ.0m.linear[i], na.rm = TRUE)
} else {
my.xlim = range(aop[, i], aop.0[i], na.rm = TRUE)
}
plot(aop.all[, i], Depth.all, type = "p", log = "x",
xlim = my.xlim,
ylim = c(max(Z.interval,na.rm = T)+0.5,0),
pch = 19, xlab = "", ylab = "",
axes = FALSE, frame.plot = TRUE,
main = substitute(E[d]*z~x~"("*mu*W.*cm^{-2}*.nm^{-1}*")"~r^2==r.2,list(x = waves[i], r.2=signif(r2[i],3))))
grid(col = 1)
lines(aop.fitted[, i], depth.fitted, col = "red")
points(aop.0[i], depth.0, pch = 1, cex=1.5, col = "red")
if (PLOT.LINEAR) lines(EdZ.fitted[, i], depth.fitted, col = "blue")
if (PLOT.LINEAR) points(EdZ.0m.linear[i],depth.0, pch = 1, cex=1.8, col = "blue")
if (PLOT.LINEAR) points(aop[ix.Z.interval[i],i], Z.interval[i], cex=1.8,
pch = 19,col ="blue")
abline(v=EdZ.detect.lim[i], col="orange", lwd=3)
abline(h=sub.surface.removed.layer.optics["EdZ"], col="green", lwd=2)
if (i == floor(seq(ix.w[1], ix.w[length(ix.w)], length.out = 4))[4])
legend("bottomright", c("LOESS", "Linear"),
lwd=c(3,3), col=c("red", "blue"))
axis(1)
axis(2)
}
}
par(mai = mai.old2)
par(mgp = mgp.old2)
par(mfrow = c(1, 1))
par(mai = mai.old1)
par(mgp = mgp.old1)
if(INTERACTIVE) x11(width = win.width, height = win.height)
matplot(KZ.fitted, depth.fitted[-1], type = "l",
lty = 1, ylim = rev(range(Depth)),
xlab = expression(K[z]~"("*E[d]*"z) (local K, i.e. derivative of Ed(z))"),
ylab = "Depth (m)",col = aop.cols)
grid(col = 1)
par(xpd = TRUE)
legend(par("usr")[1], par("usr")[4], legend = waves, xjust = 0, yjust = 0, lty = 1, lwd = 2, col = aop.cols, ncol = ceiling(length(waves) / 2), cex = 0.75)
par(xpd = FALSE)
if(INTERACTIVE) x11(width = win.width, height = win.height)
if (any(!is.na(K0.fitted))){
matplot(K0.fitted, depth.fitted[-1], type = "l",
lty = 1, ylim = rev(range(Depth)),
xlab = expression(K[0]~"("*E[d]*"z) (Depth-integrated from 0- to Z"),
ylab = "Depth (m)", col = aop.cols)
}
grid(col = 1)
par(xpd = TRUE)
legend(par("usr")[1], par("usr")[4], legend = waves, xjust = 0, yjust = 0, lty = 1, lwd = 2, col = aop.cols, ncol = ceiling(length(waves) / 2), cex = 0.75)
par(xpd = FALSE)
if (EXTRAPOLATION.0m) {
### Plot PAR
### The last depth is remove be cause it is often bad
par(mfrow = c(1, 2))
plot(PAR$PAR.z[1:(n.fitted-1)], depth.fitted[1:(n.fitted-1)],
type = "l", lwd=2,
lty = 1, ylim = rev(range(depth.fitted)), log="x",
xlab = expression(PAR~"("*mu*mol*.photon*.m^{-2}*.s^{-1}*")"),
ylab = "Depth (m)", main = "Downwelling Photosynthetic Active Radiation")
plot(PAR$PAR.z[1:(n.fitted-1)]/PAR$PAR.z[1]*100, depth.fitted[1:(n.fitted-1)],
type = "l", lwd=2,
lty = 1, ylim = rev(range(depth.fitted)),
xlab = expression(PAR~"(%)"), log="x",
ylab = "Depth (m)")
for (i in 1:5){
if (PAR$z.f.PAR[i,2]<max(depth.fitted)) {
points(PAR$z.f.PAR[i,], pch=19, col=2, cex=1.5)
posi=4
if (i == 1) posi = 2
text(PAR$z.f.PAR[i,], paste(PAR$z.f.PAR[i,1],
"% PAR at z=",
signif(PAR$z.f.PAR[i,2],3)),
pos=posi)
}
}
par(mfrow = c(1, 1))
}
if (!PLOT.LINEAR) Z.interval <-NA
return(list(
EdZ.fitted = aop.fitted,
KZ.EdZ.fitted = KZ.fitted,
K0.EdZ.fitted = K0.fitted,
EdZ.0m = aop.0,
K.EdZ.surf = K.surf,
EdZ.Z.interval = Z.interval,
EdZ.ix.Z.interval = ix.Z.interval,
EdZ.0m.linear = EdZ.0m.linear,
EdZ.linear.r2 = r2,
EdZ.KolmolSmirnov.p.value = KolmolSmirnov.p.value,
EdZ.detection.limit = EdZ.detect.lim,
PARd.z = PAR$PAR.z,
PARd.at.z = PAR$PAR.at.z
))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.