Fix the glathida processing part (#66)

* fix the glathida processing part

* bump version

Author: albangossard

Project.toml

@@ -1,7 +1,7 @@
 name = "Sleipnir"
 uuid = "f5e6c550-199f-11ee-3608-394420200519"
 authors = ["Jordi Bolibar <[email protected]>", "Facundo Sapienza <[email protected]>"]
-version = "0.7.0"
+version = "0.7.1"
 
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"

src/glaciers/glacier/glacier2D_utils.jl

@@ -36,7 +36,7 @@ function initialize_glaciers(rgi_ids::Vector{String}, params::Parameters; test=f
         pmap((rgi_id) -> generate_raw_climate_files(rgi_id, params.simulation), rgi_ids)
     end
 
-    glaciers = pmap((rgi_id) -> initialize_glacier(rgi_id, params; smoothing=false, test=test), rgi_ids)
+    glaciers::Vector{Glacier2D} = pmap((rgi_id) -> initialize_glacier(rgi_id, params; smoothing=false, test=test), rgi_ids)
 
     if params.simulation.use_glathida_data == true
 
@@ -86,11 +86,6 @@ function initialize_glacier(rgi_id::String, parameters::Parameters; smoothing=fa
     # Initialize glacier climate
     initialize_glacier_climate!(glacier, parameters)
 
-    if test
-        glacier.rgi_id = nothing # not sure of that line
-        glacier.S_coords = nothing
-    end
-
     return glacier
 end
 
@@ -193,23 +188,23 @@ function get_glathida!(glaciers::Vector{Glacier2D}, params::Parameters; force=fa
 end
 
 function get_glathida_glacier(glacier::Glacier2D, params::Parameters, force)
-    rgi_path = joinpath(prepro_dir, params.simulation.rgi_paths[rgi_id])
+    rgi_path = joinpath(prepro_dir, params.simulation.rgi_paths[glacier.rgi_id])
     gtd_path = joinpath(rgi_path, "glathida.h5")
     if isfile(gtd_path) && !force
         gtd_grid = h5read(gtd_path, "gtd_grid")
     else
-        glathida = CSV.File(joinpath(rgi_path, "glathida.csv"))
+        glathida = CSV.File(joinpath(rgi_path, "glathida_data.csv"))
         gtd_grid = zeros(size(glacier.H₀))
         count = zeros(size(glacier.H₀))
-        for (thick, i, j) in zip(glathida["elevation"], glathida["i_grid"], glathida["j_grid"])
+        for (thick, i, j) in zip(glathida["thickness"], glathida["i_grid"], glathida["j_grid"])
             count[i,j] += 1
             gtd_grid[i,j] += thick
         end
 
-        gtd_grid .= ifelse.(count > 0, gtd_grid ./ count, 0.0)
+        gtd_grid .= ifelse.(count .> 0, gtd_grid ./ count, 0.0)
 
         # Save file
-        h5open(joinpath(prepro_dir, params.simulation.rgi_paths[glacier.rgi_id], "glathida.h5"), "w") do file
+        h5open(joinpath(rgi_path, "glathida.h5"), "w") do file
             write(file, "gtd_grid", gtd_grid)
         end
     end

test/data/glaciers/glaciers2D_w_glathida.jld2

@@ -1,23 +1,31 @@
 
 
-function glaciers2D_constructor(; save_refs::Bool = false)
+function glaciers2D_constructor(; save_refs::Bool = false, use_glathida_data::Bool = false)
 
     rgi_paths = get_rgi_paths()
-    rgi_ids = ["RGI60-11.03638", "RGI60-11.01450"]
+    if use_glathida_data
+        rgi_ids = ["RGI60-07.00042", "RGI60-07.00065"] # Use glaciers that have glathida data
+        file_suffix = "w_glathida"
+    else
+        rgi_ids = ["RGI60-11.03638", "RGI60-11.01450"]
+        file_suffix = "wo_glathida"
+    end
+    # Filter out glaciers that are not used to avoid having references that depend on all the glaciers processed in Gungnir
+    rgi_paths = Dict(k => rgi_paths[k] for k in rgi_ids)
 
     params = Parameters(simulation=SimulationParameters(velocities=false,
-                                                        use_glathida_data=false,
+                                                        use_glathida_data=use_glathida_data,
                                                         working_dir=Sleipnir.root_dir,
                                                         test_mode=true,
                                                         rgi_paths=rgi_paths))
 
     glaciers = initialize_glaciers(rgi_ids, params; test=true)
 
     if save_refs
-        jldsave(joinpath(Sleipnir.root_dir, "test/data/glaciers/glaciers2D.jld2"); glaciers)
+        jldsave(joinpath(Sleipnir.root_dir, string("test/data/glaciers/glaciers2D_", file_suffix, ".jld2")); glaciers)
     end
 
-    glaciers_ref = load(joinpath(Sleipnir.root_dir,"test/data/glaciers/glaciers2D.jld2"))["glaciers"]
+    glaciers_ref = load(joinpath(Sleipnir.root_dir, string("test/data/glaciers/glaciers2D_", file_suffix, ".jld2")))["glaciers"]
 
     @test all(glaciers == glaciers_ref)
 

test/data/glaciers/glaciers2D.jld2 test/data/glaciers/glaciers2D_wo_glathida.jld2

@@ -1,7 +1,10 @@
 
 function params_constructor_specified(; save_refs::Bool = false)
 
+    rgi_id = "RGI60-11.03638"
     rgi_paths = get_rgi_paths()
+    # Filter out glaciers that are not used to avoid having references that depend on all the glaciers processed in Gungnir
+    rgi_paths = Dict(rgi_id => rgi_paths[rgi_id])
 
     physical_params = PhysicalParameters(ρ = 900.0,
                                         g = 9.81,

test/data/params/params_specified.jld2

@@ -21,6 +21,8 @@ ENV["GKSwstype"]="nul"
 
 @testset "Parameters constructors by default" params_constructor_default()
 
-@testset "Glaciers 2D constructors" glaciers2D_constructor()
+@testset "Glaciers 2D constructors w/o glathida data" glaciers2D_constructor(use_glathida_data=false)
+
+@testset "Glaciers 2D constructors w/ glathida data" glaciers2D_constructor(use_glathida_data=true)
 
 #@testset "Glaciers 2D plots" glaciers2D_plots()