[Help needed] Correlation length and lengraddepth

[ctroupin]

Request from user

I am working on regional 6y maps for DIN in the Baltic.

Last time I had fixed correlation lengths and I want to improve the result by for example limit the correlation length near the coast with the function lengraddepth. I don’t understand how to choose the parameters L and hmin?

For example: with L = 1000 there is a sharp border at the coast but with L = 10_000 there seems to also be limitations in the open sea:

I also tried fitvertcorrlen but it generated lenz varying between 2 and over 100 m and the result was very patchy. When I tried fithorzcorrlen – the code crashed.

I have tried a lenz dependent on the depth resolution but I get really strange gradients near the surface with high values at 0 m and a decrease at 5 m. I remember I got this result also last time and that is why I ended up with a fixed lenz=25m. In the Baltic in winter, the surface is well mixed down to the halocline and it would be nice to have this result also from the diva model.

An example: this is DIN-profiles from a position in the open sea for the winter. Lenz is twice the depth resolution and in the surface lenz is thus 10 m. There are very sharp gradients in the surface that are not present in observations, do you have an idea on what to tune to get a more realistic result? If I increase lenz in the surface it gets better but still there are these gradients. I also tried to do a run with lenz=0 but that was not good. (lenx=leny=60_000*RL)

Observations from this area (winter 2000-2006):

[jmbeckers]

For the first question, https://github.com/gher-uliege/DIVAnd.jl/blob/b750052d068f59c1ad5a7cb5eb85258f5fabd3c2/src/utils.jl#L468

provides the formula R_L = 1 / (1 + L |∇h| / max(h2,hmin))

In a nutshell h/|∇h| is a length scale but prone to noise. (specially when you have large topography gradients in shallow regions). That explains why you replace h by max(h2,hmin), where h2 could be a filtered version of h but by default is h, and hmin limits the depth. L is a reference value to which you compare h/|∇h| so as to have a relative length scale. So it is normal that you can have some variations in deeper layers too. If you only want to detect/include that you are close to the coast, I think there are some tools based on distance to coast that might be better.

For the maximum near the surface: when you say "I also tried to do a run with lenz=0 but that was not good. (lenx=leny=60_000*RL)" do you mean it had the same problem or just that the results were not "nice" for other reasons. If it is the former, then there is a problem with the data (bad data or very strange different spatial distribution in the surface and just below ?) not the vertical correlations (since they are deactivated when lenz=0).

If it is the vertical correlation that causes problems, you might try to use surfextend=true .

Another try could be to use option alphabc=0 (reverting to the old way boundary conditions were applied)

For the problems with the fitting, difficult to say without proper error identification

[ctroupin]

Yes, maybe 'crashed' is the wrong word, it works for some time and then stops. Maybe it is to few observations to perform the analyse with fithorzcorrlen.

Error log here:

timeindex = 3
[ Info: Time step 4 / 10
+ Warning: Observations equal to NaN: 290
+ @ DIVAnd C:\Users\a001819\.julia\packages\DIVAnd\MV3j9\src\DIVAnd_obs.jl:45
+ Warning: Too few data. Will use guesses (np = 4, RLz = 27501.283175834375, )
+ @ DIVAnd C:\Users\a001819\.julia\packages\DIVAnd\MV3j9\src\fit.jl:900
[ Info: Data points at z=0.0: 16501, horz. correlation length: 27501.283175834375 (preliminary)
+ Warning: Too few data. Will use guesses (np = 4, RLz = 27545.121879822887, )
+ @ DIVAnd C:\Users\a001819\.julia\packages\DIVAnd\MV3j9\src\fit.jl:900
[ Info: Data points at z=5.0: 17006, horz. correlation length: 27545.121879822887 (preliminary)
+ Warning: Too few data. Will use guesses (np = 4, RLz = 27561.501312525033, )
+ @ DIVAnd C:\Users\a001819\.julia\packages\DIVAnd\MV3j9\src\fit.jl:900
[ Info: Data points at z=10.0: 17050, horz. correlation length: 27561.501312525033 (preliminary)
+ Warning: Too few data. Will use guesses (np = 4, RLz = 27623.6639470855, )
+ @ DIVAnd C:\Users\a001819\.julia\packages\DIVAnd\MV3j9\src\fit.jl:900
[ Info: Data points at z=15.0: 17471, horz. correlation length: 27623.6639470855 (preliminary)
+ Warning: Too few data. Will use guesses (np = 4, RLz = 27655.63998846335, )
+ @ DIVAnd C:\Users\a001819\.julia\packages\DIVAnd\MV3j9\src\fit.jl:900
[ Info: Data points at z=20.0: 17505, horz. correlation length: 27655.63998846335 (preliminary)
+ Warning: Too few data. Will use guesses (np = 4, RLz = 27654.638446248297, )
+ @ DIVAnd C:\Users\a001819\.julia\packages\DIVAnd\MV3j9\src\fit.jl:900
[ Info: Data points at z=30.0: 17891, horz. correlation length: 27654.638446248297 (preliminary)
+ Warning: Too few data. Will use guesses (np = 4, RLz = 27720.887252623797, )
+ @ DIVAnd C:\Users\a001819\.julia\packages\DIVAnd\MV3j9\src\fit.jl:900
[ Info: Data points at z=40.0: 18302, horz. correlation length: 27720.887252623797 (preliminary)
+ Warning: Too few data. Will use guesses (np = 4, RLz = 28038.14847524304, )
+ @ DIVAnd C:\Users\a001819\.julia\packages\DIVAnd\MV3j9\src\fit.jl:900
[ Info: Data points at z=50.0: 17661, horz. correlation length: 28038.14847524304 (preliminary)
[ Info: Data points at z=60.0: 8539, horz. correlation length: 18781.950834440657 (preliminary)
[ Info: Data points at z=70.0: 5272, horz. correlation length: 134684.13053460122 (preliminary)
[ Info: Data points at z=80.0: 3837, horz. correlation length: 144931.0583987114 (preliminary)
[ Info: Data points at z=90.0: 2747, horz. correlation length: 67156.45478064583 (preliminary)
[ Info: Data points at z=100.0: 2092, horz. correlation length: 177597.39818494118 (preliminary)
[ Info: Data points at z=125.0: 1320, horz. correlation length: 68759.46545386632 (preliminary)
[ Info: Data points at z=150.0: 583, horz. correlation length: 24763.051291483564 (preliminary)
[ Info: Data points at z=175.0: 444, horz. correlation length: 15728.377900419542 (preliminary)
[ Info: Data points at z=200.0: 305, horz. correlation length: 3108.3430110719532 (preliminary)
[ Info: Data points at z=225.0: 216, horz. correlation length: 236246.32410298294 (preliminary)
[ Info: Data points at z=250.0: 101, horz. correlation length: 176581.87588485502 (preliminary)
[ Info: Data points at z=275.0: 113, horz. correlation length: 175791.3580190258 (preliminary)
TaskFailedException

    nested task error: InexactError: trunc(Int64, NaN)
    Stacktrace:
     [1] trunc
       @ .\float.jl:781 [inlined]
     [2] floor
       @ .\float.jl:357 [inlined]
     [3] fitlen(x::Tuple{Vector{Float64}, Vector{Float64}}, d::Vector{Float64}, weight::Vector{Float64}, nsamp::Int64, iter::DIVAnd.RandomCoupels{Random._GLOBAL_RNG}; distfun::typeof(DIVAnd.distfun_m), iseed::Int64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
       @ DIVAnd C:\Users\a001819\.julia\packages\DIVAnd\MV3j9\src\fit.jl:778
     [4] fitlen(x::Tuple{Vector{Float64}, Vector{Float64}}, d::Vector{Float64}, weight::Vector{Float64}, nsamp::Int64; rng::Random._GLOBAL_RNG, kwargs::Base.Pairs{Symbol, typeof(DIVAnd.distfun_m), Tuple{Symbol}, NamedTuple{(:distfun,), Tuple{typeof(DIVAnd.distfun_m)}}})
       @ DIVAnd C:\Users\a001819\.julia\packages\DIVAnd\MV3j9\src\fit.jl:683
     [5] macro expansion
       @ C:\Users\a001819\.julia\packages\DIVAnd\MV3j9\src\fit.jl:1074 [inlined]
     [6] (::DIVAnd.var"#4244#threadsfor_fun#585"{DIVAnd.var"#4244#threadsfor_fun#577#586"{Float64, typeof(DIVAnd.distfun_m), Bool, Random._GLOBAL_RNG, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Int64, Vector{Dict{Symbol, Any}}, Vector{Float64}, Vector{Float64}, UnitRange{Int64}}})(tid::Int64; onethread::Bool)
       @ DIVAnd .\threadingconstructs.jl:84
     [7] #4244#threadsfor_fun
       @ .\threadingconstructs.jl:51 [inlined]
     [8] (::Base.Threads.var"#1#2"{DIVAnd.var"#4244#threadsfor_fun#585"{DIVAnd.var"#4244#threadsfor_fun#577#586"{Float64, typeof(DIVAnd.distfun_m), Bool, Random._GLOBAL_RNG, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Int64, Vector{Dict{Symbol, Any}}, Vector{Float64}, Vector{Float64}, UnitRange{Int64}}}, Int64})()
       @ Base.Threads .\threadingconstructs.jl:30

Stacktrace:
  [1] wait
    @ .\task.jl:345 [inlined]
  [2] threading_run(fun::DIVAnd.var"#4244#threadsfor_fun#585"{DIVAnd.var"#4244#threadsfor_fun#577#586"{Float64, typeof(DIVAnd.distfun_m), Bool, Random._GLOBAL_RNG, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Int64, Vector{Dict{Symbol, Any}}, Vector{Float64}, Vector{Float64}, UnitRange{Int64}}}, static::Bool)
    @ Base.Threads .\threadingconstructs.jl:38
  [3] macro expansion
    @ .\threadingconstructs.jl:89 [inlined]
  [4] fithorzlen(x::Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}, value::Vector{Float64}, z::Vector{Float64}; tolrel::Float64, smoothz::Float64, smoothk::Int64, searchz::Float64, progress::Function, distfun::Function, limitfun::DIVAnd.var"#583#592", maxnsamp::Int64, limitlen::Bool, epsilon2::Vector{Float64}, min_rqual::Float64, rng::Random._GLOBAL_RNG)
    @ DIVAnd C:\Users\a001819\.julia\packages\DIVAnd\MV3j9\src\fit.jl:1060
  [5] (::DIVAnd.var"#447#463"{Vector{Float64}, Dict{Any, Any}, Dict{Any, Any}, typeof(DIVAnd.distfun_m), Nothing, Int64, Bool, Bool, Dict{Any, Any}, Dict{Any, Any}, Int64, Float64, Int64, Float64, Float64, Tuple{}, Bool, Vector{Tuple{String, Float64}}, typeof(DIVAndgo), Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:solver, :MEMTOFIT), Tuple{Symbol, Int64}}}, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{DateTime}}, Vector{Float64}, String, String, Bool, Matrix{DateTime}, Dict{Symbol, Any}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Vector{Float64}, Vector{Float64}, Int64, Tuple{Int64, Int64, Int64}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, DIVAnd.Anam.var"#invtrans#12"{DIVAnd.Anam.var"#invtrans#9#13"{Float64, Float64}}, DIVAnd.Anam.var"#trans#10"{DIVAnd.Anam.var"#trans#8#11"{Float64, Float64}}, Vector{DateTime}, TimeSelectorYearListMonthList{Vector{UnitRange{Int64}}, Vector{Vector{Int64}}}, Vector{Float64}, Int64})(ds::NCDataset{Nothing})
    @ DIVAnd C:\Users\a001819\.julia\packages\DIVAnd\MV3j9\src\diva.jl:551
  [6] NCDataset(::DIVAnd.var"#447#463"{Vector{Float64}, Dict{Any, Any}, Dict{Any, Any}, typeof(DIVAnd.distfun_m), Nothing, Int64, Bool, Bool, Dict{Any, Any}, Dict{Any, Any}, Int64, Float64, Int64, Float64, Float64, Tuple{}, Bool, Vector{Tuple{String, Float64}}, typeof(DIVAndgo), Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:solver, :MEMTOFIT), Tuple{Symbol, Int64}}}, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{DateTime}}, Vector{Float64}, String, String, Bool, Matrix{DateTime}, Dict{Symbol, Any}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Vector{Float64}, Vector{Float64}, Int64, Tuple{Int64, Int64, Int64}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, DIVAnd.Anam.var"#invtrans#12"{DIVAnd.Anam.var"#invtrans#9#13"{Float64, Float64}}, DIVAnd.Anam.var"#trans#10"{DIVAnd.Anam.var"#trans#8#11"{Float64, Float64}}, Vector{DateTime}, TimeSelectorYearListMonthList{Vector{UnitRange{Int64}}, Vector{Vector{Int64}}}, Vector{Float64}, Int64}, ::String, ::Vararg{String}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ NCDatasets C:\Users\a001819\.julia\packages\NCDatasets\sLdiM\src\dataset.jl:241
  [7] NCDataset(::Function, ::String, ::Vararg{String})
    @ NCDatasets C:\Users\a001819\.julia\packages\NCDatasets\sLdiM\src\dataset.jl:238
  [8] diva3d(xi::Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Vector{Float64}, TimeSelectorYearListMonthList{Vector{UnitRange{Int64}}, Vector{Vector{Int64}}}}, x::Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{DateTime}}, value::Vector{Float64}, len::Tuple{Array{Int64, 3}, Array{Int64, 3}, Array{Float64, 3}}, epsilon2::Float64, filename::String, varname::String; datadir::String, bathname::String, bathisglobal::Bool, plotres::var"#plotres#6"{Int64}, timeorigin::DateTime, moddim::Vector{Float64}, zlevel::Symbol, ncvarattrib::Dict{Any, Any}, ncglobalattrib::Dict{Any, Any}, transform::Tuple{DIVAnd.Anam.var"#trans#10"{DIVAnd.Anam.var"#trans#8#11"{Float64, Float64}}, DIVAnd.Anam.var"#invtrans#12"{DIVAnd.Anam.var"#invtrans#9#13"{Float64, Float64}}}, distfun::typeof(DIVAnd.distfun_m), mask::BitArray{3}, background::Nothing, background_epsilon2_factor::Nothing, background_lenz::Nothing, background_len::Nothing, background_lenz_factor::Int64, filterbackground::Int64, fitcorrlen::Bool, fithorzcorrlen::Bool, fitvertcorrlen::Bool, fithorz_param::Dict{Any, Any}, fitvert_param::Dict{Any, Any}, memtofit::Int64, overlapfactor::Float64, niter_e::Int64, minfield::Float64, maxfield::Float64, surfextend::Bool, velocity::Tuple{}, stat_per_timeslice::Bool, error_thresholds::Vector{Tuple{String, Float64}}, divamethod::typeof(DIVAndgo), kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:solver, :MEMTOFIT), Tuple{Symbol, Int64}}})
    @ DIVAnd C:\Users\a001819\.julia\packages\DIVAnd\MV3j9\src\diva.jl:375
  [9] macro expansion
    @ .\timing.jl:262 [inlined]
 [10] top-level scope
    @ .\In[37]:62

[jmbeckers]

Looks like multithreading is active ? Is that on purpose there ? is z=275.0 the last layer ?

[Alexander-Barth]

There are a lot of debug statement in fit.jl. Maybe we get some more information when we activate them with ENV["JULIA_DEBUG"] = "DIVAnd".

https://github.com/gher-uliege/DIVAnd.jl/blob/v2.7.9/src/fit.jl#L778

Maybe we have all observations at the same location around z=275.0 ? (in case, fitlen should handle even pathological situation more gracefully :-))

[jmbeckers]

For the problem at the top, can you also do a test with alphabc=0 as additional keyword argument ?

[kawess]

Thanks! Can you please guide me to the "tools based on distance to coast that might be better"? Perhaps that works better than lengraddepth.

The surface layer improved when I used a specific background so I definitelly need to use that. I am still not totally happy so I need to work more. I have not decided yet if I should use a background for the season and 6 years or the season and all years (more lika a climatology background). I also use surfextend=true and alphabc=0 as you suggested.

[jmbeckers]

Here is a piece of code (from Emodnet Chemistry), which could be used to create a correlation length lenfilled decreasing to a fifth of the normal value lenfwhen approaching a cost. The distance at which the change happens is defined by slen (here in meters, so metrics pmn are in /m). sz is size(mask)

The mask is 2D.

nearcoast = 1.0 * (1. .- mask)
nearcoastf = copy(nearcoast)

slen = (100e3,100e3)
nearcoastf = @time DIVAnd.diffusion(
    trues(sz),pmn,slen,nearcoast;
    boundary_condition! = x -> x[.!mask] .= 1
)

lenf2 = (1 .- nearcoastf) .* lenf + nearcoastf .* lenf/5

#@code_warntype DIVAnd.diffusion(mask,pmn,slen,len)
lenf2[.!mask] .= NaN

lenfilled = DIVAnd.ufill(lenf2,isfinite.(lenf2));

[GeoTuxMan]

Hello !

I use this code to reduce the CL near the coast: `x,y = DIVAnd.ndgrid(bx,by);

pm,pn = DIVAnd.DIVAnd_metric(x,y);

L = 11_100 # meters RL = DIVAnd.lengraddepth((pm,pn),b, L); RL[isnan.(RL)].=0.0 RL3D = repeat(RL, inner=(1,1,length(depthr))) #repeated for the different depth levels. @show size(RL3D) len=(RL3D,RL3D,RL3D);`

The resulted values are between 0.0 and 1.0; why are they so small ? are they in degrees ?

[jmbeckers]

The tool provides, as stated in the documentation

 relative correlation length-scale field `RL`

which should have values around 1. For the "real" correlation length, you have to multiply by a global correlation length in the units you work with.

The L=11_000 is just a parameter used in the calculation of RL, see documentation and explanation above.

[kawess]

If I have different normal values for lenf in x- and y-direction, should I do one lenfilledx and one lenfilledy?

[jmbeckers]

With a relative length scale you can still apply a different global scale in each direction if you have a specific reason for that. In any case, you need to fill in a field for each direction.

[kawess]

Ah ok, so resulting lenfilled (in the example above) is the actual correlation length... but... I can do: RL = lenfilled ./ lenf and get the relative length scale? I get a funny image:

[ctroupin]

If you plot by masking the land-values (with the mask created from the bathymetry), you would get rid of the weird parts.

[jmbeckers]

Yes, they are just filled to avoid any NaN in the L fields. Once you looked at the masked field, you can check if the width of the smaller L values is sufficient (otherwise increase the slen value in diffusion(....))

[jmbeckers]

Do we have any feedback on the debugging in fit.jl ? If not, should we simply change

nbmax = floor(Int, maxdist / ddist + 1)

into

nbmax=1
worktmp=maxdist / ddist 
if !isnan(worktmp)&&isfinite(worktmp)
nbmax = floor(Int, worktmp + 1)
end

and hope that later the tests of possible problems will catch that with one bin no fit should be made.

[jmbeckers]

No more feedback so considered closed