GLEON / rLakeAnalyzer

An R version of Lake Analyzer
43 stars 26 forks source link

Thermocline depth calculation #118

Open vikifc opened 9 months ago

vikifc commented 9 months ago

HI! I am finding a problem when calculating the thermocline depth of a temperature profile in R. Here it goes:

wtr = c(13.9, 13.8, 13.7, 13.7, 13.5, 13.3, 13.1, 13.0, 12.9, 12.8) water.density(wtr) [1] 999.2871 999.3007 999.3143 999.3143 999.3411 999.3674 999.3931 999.4059 999.4184 999.4309 depths = c(0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5) thermo.depth(wtr, depths, Smin = 0.5, seasonal = TRUE, mixed.cutoff = 1) [1] 1.991029

I don't understand why it finds a thermocline at 1.99 m with the Smin=0.5 condition. I also tried with really high and low Smin numbers and the result is the same, do someone know how to fix that?

I would really appreciate the help, thanks!

hdugan commented 9 months ago

Smin is a bit confusing, because it is easily overridden in the code if it's set too high or too low.

What are you trying to accomplish? Would the mixed.cutoff parameter be better suited to your goal?

vikifc commented 9 months ago

I think mixed.cutoff considers the maximum difference between de highest temperature and the lowest temperature in the whole profile, right? I am trying to calculate the depth of the thermocline, where both conditions are met: 1) the density gradient is at least 0.5 kg/m3/m, and 2) the temperature gradient is at least 1ºC/m. At 1.99 m none of that is true.

hdugan commented 9 months ago

Yeah I seem the problem. Without recoding the current function, you could apply your own density cutoff. Something like the following:

library(LakeMetabolizer)
library(tidyverse)
wtr = c(13.9, 13.8, 13.7, 13.7, 13.5, 13.3, 13.1, 13.0, 12.9, 12.8)
depths = c(0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5)

# Derivative function
drv <- function(x, y) c(NA, (diff(y) /diff(x))) 

# Create your own density cutoff 
data.frame(wtr = wtr, depths = depths) |> 
  mutate(dens = water.density(wtr)) |> 
  mutate(dens.der = drv(depths, dens)) |> 
  mutate(thermo.depth = thermo.depth(wtr,depths)) |> 
  mutate(thermo.depth = if_else(max(dens.der, na.rm = T) < 0.5, NA, thermo.depth[1]))

Using the data.frame structure probably makes more sense if you have multiple profiles and group them accordingly by date. But hopefully this is useful.

vikifc commented 9 months ago

Yes, I think I'll go for a solution like that. Thank you for your help!