fix various bugs

src/Sleipnir.jl

@@ -19,6 +19,7 @@ using CSV
 using JSON
 using CodecZlib
 using Tar
+import NCDatasets
 
 # ##############################################
 # ############    PARAMETERS     ###############

src/glaciers/climate/Climate1D.jl

@@ -1,23 +1,33 @@
 
-@kwdef mutable struct Climate1Dstep{F <: AbstractFloat} 
+@kwdef mutable struct Climate1Dstep{F <: AbstractFloat}
     temp::Vector{F}
     PDD::Vector{F}
     snow::Vector{F}
     rain::Vector{F}
-    gradient::Ref{F}
-    avg_gradient::Ref{F}
+    gradient::F
+    avg_gradient::F
     x::Vector{F}
     y::Vector{F}
-    ref_hgt::Ref{F}
+    ref_hgt::F
 end
 
-@kwdef mutable struct Climate1D{F <: AbstractFloat} 
+Base.:(==)(a::Climate1Dstep, b::Climate1Dstep) = a.temp == b.temp && a.PDD == b.PDD &&
+                                      a.snow == b.snow && a.rain == b.rain &&
+                                      a.gradient == b.gradient && a.avg_gradient == b.avg_gradient &&
+                                      a.x == b.x && a.y == b.y && a.ref_hgt == b.ref_hgt
+
+@kwdef mutable struct Climate1D{F <: AbstractFloat}
     raw_climate::RasterStack # Raw climate dataset for the whole simulation
     # Buffers to avoid memory allocations
-    climate_raw_step::Ref{RasterStack} # Raw climate trimmed for the current step
-    climate_step::Ref{RasterStack} # Climate data for the current step
+    climate_raw_step::RasterStack # Raw climate trimmed for the current step
+    climate_step::Dict # Climate data for the current step
     climate_2D_step::Climate2Dstep # 2D climate data for the current step to feed to the MB model
     longterm_temps::Vector{F} # Longterm temperatures for the ice rheology
-    avg_temps::Ref{RasterStack} # Intermediate buffer for computing average temperatures
-    avg_gradients::Ref{RasterStack} # Intermediate buffer for computing average gradients
+    avg_temps::F # Intermediate buffer for computing average temperatures
+    avg_gradients::F # Intermediate buffer for computing average gradients
 end
+
+Base.:(==)(a::Climate1D, b::Climate1D) = a.raw_climate == b.raw_climate && a.climate_raw_step == b.climate_raw_step &&
+                                      a.climate_step == b.climate_step && a.climate_2D_step == b.climate_2D_step &&
+                                      a.longterm_temps == b.longterm_temps && a.avg_temps == b.avg_temps &&
+                                      a.avg_gradients == b.avg_gradients

src/glaciers/climate/Climate2D.jl

@@ -1,25 +1,35 @@
 
 export Climate2Dstep, Climate2D
 
-@kwdef mutable struct Climate2Dstep{F <: AbstractFloat} 
+@kwdef mutable struct Climate2Dstep{F <: AbstractFloat}
     temp::Matrix{F}
     PDD::Matrix{F}
     snow::Matrix{F}
     rain::Matrix{F}
-    gradient::Ref{F}
-    avg_gradient::Ref{F}
+    gradient::F
+    avg_gradient::F
     x::Vector{F}
     y::Vector{F}
-    ref_hgt::Ref{F}
+    ref_hgt::F
 end
 
-@kwdef mutable struct Climate2D{F <: AbstractFloat} 
+Base.:(==)(a::Climate2Dstep, b::Climate2Dstep) = a.temp == b.temp && a.PDD == b.PDD &&
+                                      a.snow == b.snow && a.rain == b.rain &&
+                                      a.gradient == b.gradient && a.avg_gradient == b.avg_gradient &&
+                                      a.x == b.x && a.y == b.y && a.ref_hgt == b.ref_hgt
+
+@kwdef mutable struct Climate2D{F <: AbstractFloat}
     raw_climate::RasterStack # Raw climate dataset for the whole simulation
     # Buffers to avoid memory allocations
-    climate_raw_step::Ref{RasterStack} # Raw climate trimmed for the current step
-    climate_step::Ref{RasterStack} # Climate data for the current step
+    climate_raw_step::RasterStack # Raw climate trimmed for the current step
+    climate_step::Dict # Climate data for the current step
     climate_2D_step::Climate2Dstep # 2D climate data for the current step to feed to the MB model
     longterm_temps::Vector{F} # Longterm temperatures for the ice rheology
-    avg_temps::Ref{RasterStack} # Intermediate buffer for computing average temperatures
-    avg_gradients::Ref{RasterStack} # Intermediate buffer for computing average gradients
+    avg_temps::F # Intermediate buffer for computing average temperatures
+    avg_gradients::F # Intermediate buffer for computing average gradients
 end
+
+Base.:(==)(a::Climate2D, b::Climate2D) = a.raw_climate == b.raw_climate && a.climate_raw_step == b.climate_raw_step &&
+                                      a.climate_step == b.climate_step && a.climate_2D_step == b.climate_2D_step &&
+                                      a.longterm_temps == b.longterm_temps && a.avg_temps == b.avg_temps &&
+                                      a.avg_gradients == b.avg_gradients

src/glaciers/climate/climate2D_utils.jl

@@ -72,12 +72,16 @@ function get_cumulative_climate!(climate, period, gradient_bounds=[-0.009, -0.00
 end
 
 function get_cumulative_climate(climate, gradient_bounds=[-0.009, -0.003], default_grad=-0.0065)
-    climate.temp.data .= max.(climate.temp.data, 0.0) # get PDDs
-    climate.gradient.data .= clamp.(climate.gradient.data, gradient_bounds[1], gradient_bounds[2]) # Clip gradients within plausible values
-    climate_sum = Dict("temp" => sum(climate.temp),
+    avg_temp = mean(climate.temp)
+    avg_gradient = mean(climate.gradient)
+    copy_climate = deepcopy(climate)
+    copy_climate.temp.data .= max.(copy_climate.temp.data, 0.0) # get PDDs
+    copy_climate.gradient.data .= clamp.(copy_climate.gradient.data, gradient_bounds[1], gradient_bounds[2]) # Clip gradients within plausible values
+    climate_sum = Dict("temp" => sum(copy_climate.temp),
                        "prcp" => sum(climate.prcp),
-                       "avg_temp" => mean(climate.temp),
-                       "avg_gradient" => mean(climate.gradient),
+                       "gradient" => sum(copy_climate.gradient),
+                       "avg_temp" => avg_temp,
+                       "avg_gradient" => avg_gradient,
                        "ref_hgt" => metadata(climate)["ref_hgt"])
     return climate_sum
 end
@@ -111,7 +115,7 @@ Applies temperature gradients to the glacier 2D climate data based on a DEM.
 """
 function apply_t_grad!(climate::RasterStack, dem::Raster)
     # We apply the gradients to the temperature
-    climate.temp.data .= climate.temp.data .+ climate.gradient.data .* (mean(dem.data[:]) .- climate.ref_hgt) 
+    climate.temp.data .= climate.temp.data .+ climate.gradient.data .* (mean(dem.data[:]) .- metadata(climate)["ref_hgt"])
 end
 
 """
@@ -147,11 +151,11 @@ function downscale_2D_climate(climate_step::Dict, glacier::Glacier2D)
                        PDD=PDD_2D,
                        snow=snow_2D,
                        rain=rain_2D,
-                       gradient=climate_step["gradient"],
-                       avg_gradient=climate_step["avg_gradient"],
-                       x=glacier.S_coords.x.data,
-                       y=glacier.S_coords.y.data,
-                       ref_hgt=climate_step["ref_hgt"])
+                       gradient=Float64(climate_step["gradient"]),
+                       avg_gradient=Float64(climate_step["avg_gradient"]),
+                       x=glacier.S_coords["x"],
+                       y=glacier.S_coords["y"],
+                       ref_hgt=Float64(climate_step["ref_hgt"]))
 
     # Apply temperature gradients and compute snow/rain fraction for the selected period
     apply_t_cumul_grad!(climate_2D_step, reshape(glacier.S, size(glacier.S))) # Reproject current S with xarray structure

src/glaciers/glacier/Glacier1D.jl

@@ -19,7 +19,7 @@ mutable struct Glacier1D{F <: AbstractFloat, I <: Integer} <: AbstractGlacier
     λ::Union{Vector{F}, Nothing}
     slope::Union{Vector{F}, Nothing}
     dist_border::Union{Vector{F}, Nothing}
-    S_coords::Union{Dict{String, Integer}, Nothing}
+    S_coords::Union{Dict{String, Vector{Float64}}, Nothing}
     Δx::Union{F, Nothing}
     Δy::Union{F, Nothing}
     nx::Union{I, Nothing}
@@ -41,7 +41,7 @@ function Glacier1D(;
     λ::Union{Vector{F}, Nothing} = nothing,
     slope::Union{Vector{F}, Nothing} = nothing,
     dist_border::Union{Vector{F}, Nothing} = nothing,
-    S_coords::Union{Dict{String, Integer}, Nothing} = nothing,
+    S_coords::Union{Dict{String, Vector{Float64}}, Nothing} = nothing,
     Δx::Union{F, Nothing} = nothing,
     Δy::Union{F, Nothing} = nothing,
     nx::Union{I, Nothing} = nothing,
@@ -64,7 +64,7 @@ function Glacier1D(;
     λ::Union{Vector{F}, Nothing} = nothing,
     slope::Union{Vector{F}, Nothing} = nothing,
     dist_border::Union{Vector{F}, Nothing} = nothing,
-    S_coords::Union{Dict{String, Integer}, Nothing} = nothing,
+    S_coords::Union{Dict{String, Vector{Float64}}, Nothing} = nothing,
     Δx::Union{F, Nothing} = nothing,
     Δy::Union{F, Nothing} = nothing,
     nx::Union{I, Nothing} = nothing,

src/glaciers/glacier/Glacier2D.jl

@@ -20,7 +20,7 @@ mutable struct Glacier2D{F <: AbstractFloat, I <: Integer} <: AbstractGlacier
     n::Union{F, Nothing}
     slope::Union{Matrix{F}, Nothing}
     dist_border::Union{Matrix{F}, Nothing}
-    S_coords::Union{Dict{String, Integer}, Nothing}
+    S_coords::Union{Dict{String, Vector{Float64}}, Nothing}
     Δx::Union{F, Nothing}
     Δy::Union{F, Nothing}
     nx::Union{I, Nothing}
@@ -45,7 +45,7 @@ function Glacier2D(;
     n::Union{F, Nothing} = nothing,
     slope::Union{Matrix{F}, Nothing} = nothing,
     dist_border::Union{Matrix{F}, Nothing} = nothing,
-    S_coords::Union{Dict{String, Integer}, Nothing} = nothing,
+    S_coords::Union{Dict{String, Vector{Float64}}, Nothing} = nothing,
     Δx::Union{F, Nothing} = nothing,
     Δy::Union{F, Nothing} = nothing,
     nx::Union{I, Nothing} = nothing,
@@ -71,7 +71,7 @@ function Glacier2D(;
     n::Union{F, Nothing} = nothing,
     slope::Union{Matrix{F}, Nothing} = nothing,
     dist_border::Union{Matrix{F}, Nothing} = nothing,
-    S_coords::Union{Dict{String, Integer}, Nothing} = nothing,
+    S_coords::Union{Dict{String, Vector{Float64}}, Nothing} = nothing,
     Δx::Union{F, Nothing} = nothing,
     Δy::Union{F, Nothing} = nothing,
     nx::Union{I, Nothing} = nothing,
@@ -105,7 +105,7 @@ Base.:(≈)(a::Glacier2D, b::Glacier2D) = a.rgi_id == b.rgi_id && a.climate == b
                                         safe_approx(a.S, b.S) && safe_approx(a.B, b.B) && safe_approx(a.V, b.V) &&
                                         safe_approx(a.A, b.A) && safe_approx(a.C, b.C) && safe_approx(a.n, b.n) &&
                                         isapprox(a.slope, b.slope; rtol=1e-3) && safe_approx(a.dist_border, b.dist_border) &&
-                                        a.S_coords == b.S_coords && safe_approx(a.Δx, b.Δx) && safe_approx(a.Δy, b.Δy) &&
+                                        safe_approx(a.S_coords, b.S_coords) && safe_approx(a.Δx, b.Δx) && safe_approx(a.Δy, b.Δy) &&
                                         safe_approx(a.nx, b.nx) && safe_approx(a.ny, b.ny) &&
                                         safe_approx(a.cenlon, b.cenlon) && safe_approx(a.cenlat, b.cenlat)
 

src/glaciers/glacier/glacier2D_utils.jl

@@ -122,15 +122,15 @@ function initialize_glacier_data(rgi_id::String, params::Parameters; smoothing=f
 
     try
         # We filter glacier borders in high elevations to avoid overflow problems
-        dist_border = glacier_gd.dis_from_border.data
+        dist_border::Matrix{Float64} = glacier_gd.dis_from_border.data
 
             # H_mask = (dist_border .< 20.0) .&& (S .> maximum(S)*0.7)
             # H₀[H_mask] .= 0.0
 
         B = glacier_gd.topo.data .- H₀ # bedrock
 
-        S_coords = Dict{"x"=> glacier_gd.topo.dims[1], "y"=> glacier_gd.topo.dims[2]}
-        S = glacier_gd.topo.data
+        S_coords = Dict{String,Vector{Float64}}("x"=> dims(glacier_gd, 1).val, "y"=> dims(glacier_gd, 2).val)
+        S::Matrix{Float64} = glacier_gd.topo.data
         #smooth!(S)
 
         if params.simulation.velocities
@@ -149,7 +149,7 @@ function initialize_glacier_data(rgi_id::String, params::Parameters; smoothing=f
         ny = glacier_grid["nxny"][2] 
         Δx = abs.(glacier_grid["dxdy"][1])
         Δy = abs.(glacier_grid["dxdy"][2])
-        slope = glacier_gd.slope.data
+        slope::Matrix{Float64} = glacier_gd.slope.data
 
         # We initialize the Glacier with all the initial topographical 
         glacier = Glacier2D(rgi_id = rgi_id, 

src/setup/config.jl

@@ -54,7 +54,7 @@ end
 
 function get_rgi_paths()
     rgi_paths = JSON.parsefile(joinpath(prepro_dir, "rgi_paths.json"))
-    rgi_paths = Dict(k => string(v) for (k,v) in pairs(rgi_paths)) # Convert Dict{String, Any} to Dict{String, String}
+    rgi_paths = Dict(k => joinpath(prepro_dir, string(v)) for (k,v) in pairs(rgi_paths)) # Convert Dict{String, Any} to Dict{String, String}
     return rgi_paths
 end
 

test/glaciers_construction.jl

@@ -19,7 +19,7 @@ function glaciers2D_constructor(; save_refs::Bool = false)
 
     glaciers_ref = load(joinpath(Sleipnir.root_dir,"test/data/glaciers/glaciers2D.jld2"))["glaciers"]
 
-    @test all(glaciers .≈ glaciers_ref)
+    @test all(glaciers == glaciers_ref)
 
 
 end

test/params_construction.jl

@@ -28,6 +28,7 @@ function params_constructor_specified(; save_refs::Bool = false)
                                             working_dir = "",
                                             rgi_paths = rgi_paths)
 
+    # TODO: check that the tests pass if script is launched from a different computer with different values inside rgi_paths
 
     params = Parameters(physical=physical_params,
                         simulation=simulation_params)