[WIP] continue the migration, tests of params_construction.jl are pas…

…sing

src/Sleipnir.jl

@@ -14,6 +14,7 @@ using Statistics
 using CairoMakie
 using Downloads
 using HDF5
+using Rasters
 
 # ##############################################
 # ############    PARAMETERS     ###############

src/glaciers/climate/climate2D_utils.jl

@@ -17,10 +17,10 @@ Initializes the `Climate` data structure for a given `Glacier``
 function initialize_glacier_climate!(glacier::AbstractGlacier, params::Parameters)
 
     dummy_period = partial_year(Day, params.simulation.tspan[1]):Day(1):partial_year(Day, params.simulation.tspan[1] + params.simulation.step)
-    raw_climate = RasterStack(joinpath(params.simulation.rgi_path, "raw_climate_$(params.simulation.tspan).nc"))
+    raw_climate = RasterStack(joinpath(params.simulation.rgi_paths[glacier.rgi_id], "raw_climate_$(params.simulation.tspan).nc"))
     climate_step = get_cumulative_climate(raw_climate[At(dummy_period)])
     climate_2D_step = downscale_2D_climate(climate_step, glacier)
-    longterm_temps = get_longterm_temps(params, raw_climate)
+    longterm_temps = get_longterm_temps(glacier.rgi_id, params, raw_climate)
     glacier.climate = Climate2D(raw_climate = raw_climate,
                             climate_raw_step = raw_climate[At(dummy_period)],
                             #climate_cum_step = raw_climate.sel(time=dummy_period).sum(),
@@ -190,21 +190,22 @@ function partial_year(period::Type{<:Period}, float)
     partial = period(round(Dates.value(year) * Δ))
     year_start + partial
 end
-partial_year(float) = partial_year(Day, float) 
+partial_year(float) = partial_year(Day, float)
 
 
-function get_longterm_temps(params::Parameters, tspan)
-    glacier_gd = RasterStack(joinpath(params.simulation.rgi_path, "gridded_data.nc"))
-    climate = RasterStack(joinpath(gdir.dir, "raw_climate_$tspan.nc"))
+function get_longterm_temps(rgi_id::String, params::Parameters)
+    rgi_path = params.simulation.rgi_paths[rgi_id]
+    glacier_gd = RasterStack(joinpath(rgi_path, "gridded_data.nc"))
+    climate = get_raw_climate_data(rgi_path)
     apply_t_grad!(climate, glacier_gd.topo)
-    longterm_temps = climate.groupby("time.year").mean().temp.data # TODO: change to Rasters.jl syntax
+    longterm_temps = mean.(groupby(climate.temp, Ti=>year)).data
     return longterm_temps
 end
 
-function get_longterm_temps(params::Parameters, climate::RasterStack)
-    glacier_gd = RasterStack(joinpath(params.simulation.rgi_path, "gridded_data.nc"))
+function get_longterm_temps(rgi_id::String, params::Parameters, climate::RasterStack)
+    glacier_gd = RasterStack(joinpath(params.simulation.rgi_paths[rgi_id], "gridded_data.nc"))
     apply_t_grad!(climate, glacier_gd.topo)
-    longterm_temps = climate.groupby("time.year").mean().temp.data # TODO: change to Rasters.jl syntax
+    longterm_temps = mean.(groupby(climate.temp, Ti=>year)).data
     return longterm_temps
 end
 

src/glaciers/glacier/Glacier.jl

@@ -1,4 +1,4 @@
 
 include("Glacier2D.jl")
 
-include("Glacier1D.jl")
+# include("Glacier1D.jl") # TODO: uncomment this line

src/glaciers/glacier/Glacier2D.jl

@@ -1,7 +1,7 @@
 
 export Glacier2D, Climate2D
 
-abstract type AbstractGlacier end 
+abstract type AbstractGlacier end
 
 include("../climate/Climate2D.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{Py, Nothing}
+    S_coords::Union{Dict{String, Integer}, Nothing}
     Δx::Union{F, Nothing}
     Δy::Union{F, Nothing}
     nx::Union{I, Nothing}
@@ -32,27 +32,32 @@ end
 """
 function Glacier2D(;
     rgi_id::Union{String, Nothing} = nothing,
-    gdir::Union{Py, Nothing} = nothing,
     climate::Union{Climate2D, Nothing} = nothing,
     H₀::Union{Matrix{F}, Nothing} = nothing,
-    H_glathida::Union{Matrix{F}, Nothing},
+    H_glathida::Union{Matrix{F}, Nothing} = nothing,
     S::Union{Matrix{F}, Nothing} = nothing,
     B::Union{Matrix{F}, Nothing} = nothing,
     V::Union{Matrix{F}, Nothing}= nothing,
+    Vx::Union{Matrix{F}, Nothing}= nothing,
+    Vy::Union{Matrix{F}, Nothing}= nothing,
+    A::Union{F, Nothing} = nothing,
+    C::Union{F, Nothing} = nothing,
+    n::Union{F, Nothing} = nothing,
     slope::Union{Matrix{F}, Nothing} = nothing,
     dist_border::Union{Matrix{F}, Nothing} = nothing,
-    S_coords::Union{Py, Nothing} = nothing,
+    S_coords::Union{Dict{String, Integer}, Nothing} = nothing,
     Δx::Union{F, Nothing} = nothing,
     Δy::Union{F, Nothing} = nothing,
     nx::Union{I, Nothing} = nothing,
-    ny::Union{I, Nothing} = nothing
-    ) where {F <: AbstractFloat, I <: Integer} 
+    ny::Union{I, Nothing} = nothing,
+    cenlon::Union{F, Nothing} = nothing,
+    cenlat::Union{F, Nothing} = nothing
+    ) where {F <: AbstractFloat, I <: Integer}
 
 Constructor for empty 2D Glacier object.
 """
 function Glacier2D(;
     rgi_id::Union{String, Nothing} = nothing,
-    gdir::Union{Py, Nothing} = nothing,
     climate::Union{Climate2D, Nothing} = nothing,
     H₀::Union{Matrix{F}, Nothing} = nothing,
     H_glathida::Union{Matrix{F}, Nothing} = nothing,
@@ -66,42 +71,42 @@ function Glacier2D(;
     n::Union{F, Nothing} = nothing,
     slope::Union{Matrix{F}, Nothing} = nothing,
     dist_border::Union{Matrix{F}, Nothing} = nothing,
-    S_coords::Union{Py, Nothing} = nothing,
+    S_coords::Union{Dict{String, Integer}, Nothing} = nothing,
     Δx::Union{F, Nothing} = nothing,
     Δy::Union{F, Nothing} = nothing,
     nx::Union{I, Nothing} = nothing,
     ny::Union{I, Nothing} = nothing,
     cenlon::Union{F, Nothing} = nothing,
     cenlat::Union{F, Nothing} = nothing
-    ) where {F <: AbstractFloat, I <: Integer} 
+    ) where {F <: AbstractFloat, I <: Integer}
 
     # Define default float and integer type for constructor
     ft = typeof(Δx)
     it = typeof(nx)
-  
-    return Glacier2D{ft,it}(rgi_id, gdir, climate, H₀, H_glathida, S, B, V, Vx, Vy, A, C, n, slope, dist_border, S_coords, Δx, Δy, nx, ny, cenlon, cenlat)
+
+    return Glacier2D{ft,it}(rgi_id, climate, H₀, H_glathida, S, B, V, Vx, Vy, A, C, n, slope, dist_border, S_coords, Δx, Δy, nx, ny, cenlon, cenlat)
 end
 
 ###############################################
 ################### UTILS #####################
 ###############################################
 
 
-Base.:(==)(a::Glacier2D, b::Glacier2D) = a.rgi_id == b.rgi_id && a.gdir == b.gdir && a.climate == b.climate && 
+Base.:(==)(a::Glacier2D, b::Glacier2D) = a.rgi_id == b.rgi_id && a.climate == b.climate &&
                                       a.H₀ == b.H₀ && a.H_glathida == b.H_glathida && a.S == b.S && a.B == b.B && a.V == b.V &&
-                                      a.A == b.A && a.C == b.C && a.n == b.n && 
-                                      a.slope == b.slope && a.dist_border == b.dist_border && 
+                                      a.A == b.A && a.C == b.C && a.n == b.n &&
+                                      a.slope == b.slope && a.dist_border == b.dist_border &&
                                       a.S_coords == b.S_coords && a.Δx == b.Δx && a.Δy == b.Δy && a.nx == b.nx && a.ny == b.ny &&
                                       a.cenlon == b.cenlon && a.cenlat == b.cenlat
 
 
-Base.:(≈)(a::Glacier2D, b::Glacier2D) = a.rgi_id == b.rgi_id && a.gdir == b.gdir && a.climate == b.climate &&
-                                        safe_approx(a.H₀, b.H₀) && safe_approx(a.H_glathida, b.H_glathida) && 
+Base.:(≈)(a::Glacier2D, b::Glacier2D) = a.rgi_id == b.rgi_id && a.climate == b.climate &&
+                                        safe_approx(a.H₀, b.H₀) && safe_approx(a.H_glathida, b.H_glathida) &&
                                         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.nx, b.nx) && safe_approx(a.ny, b.ny) && 
+                                        safe_approx(a.nx, b.nx) && safe_approx(a.ny, b.ny) &&
                                         safe_approx(a.cenlon, b.cenlon) && safe_approx(a.cenlat, b.cenlat)
 
 include("glacier2D_utils.jl")

src/glaciers/glacier/glacier2D_utils.jl

@@ -69,15 +69,14 @@ end
 
 
 """
-    initialize_glacier(gdir::Py, tspan, step; smoothing=false, velocities=true)
+    initialize_glacier(rgi_id::String, parameters::Parameters; smoothing=false, velocities=true)
+
+Initialize a single `Glacier`s, including its `Climate`, based on a `rgi_id` and timestepping arguments.
 
-Initialize a single `Glacier`s, including its `Climate`, based on a `gdir` and timestepping arguments.
-    
 Keyword arguments
 =================
-    - `gdir`: Glacier directory
-    - `tspan`: Tuple specifying the initial and final year of the simulation
-    - `step`: Step in years for the surface mass balance processing
+    - `rgi_id`: Glacier RGI ID
+    - `parameters`: Parameters including the physical and simulation ones
     - `smoothing` Flag determining if smoothing needs to be applied to the surface elevation and ice thickness.
     - `velocities` Flag determining if the ice surface velocities need to be retrieved.
 """
@@ -89,15 +88,15 @@ function initialize_glacier(rgi_id::String, parameters::Parameters; smoothing=fa
     initialize_glacier_climate!(glacier, parameters)
 
     if test
-        glacier.gdir = nothing
+        glacier.rgi_id = nothing # not sure of that line
         glacier.S_coords = nothing
     end
 
     return glacier
 end
 
 """
-    initialize_glacier(gdir::Py; smoothing=false, velocities=true)
+    initialize_glacier(rgi_id::String, params::Parameters; smoothing=false, velocities=true)
 
 Retrieves the initial glacier geometry (bedrock + ice thickness) for a glacier with other necessary data (e.g. grid size and ice surface velocities).
 """
@@ -192,20 +191,21 @@ function get_glathida!(glaciers::Vector{Glacier2D}, params::Parameters; force=fa
     end
     jldsave(joinpath(params.simulation.working_dir, "data/missing_glaciers.jld2"); missing_glaciers)
 
-   # Apply deletion to both gtd_grids and gdirs using the same set of indices
+   # Apply deletion to both gtd_grids and glaciers using the same set of indices
     indices_to_remove = findall(x -> length(x[x .!= 0.0]) == 0, gtd_grids)
     deleteat!(gtd_grids, indices_to_remove)
     deleteat!(glaciers, indices_to_remove)
-   
+
     return gtd_grids, glaciers
 end
 
 function get_glathida_glacier(glacier::Glacier2D, params::Parameters, force)
-    gtd_path = joinpath(gdir.dir, "glathida.h5")
+    rgi_path = params.simulation.rgi_paths[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(params.simulation.rgi_path, "glathida.csv")) 
+        glathida = CSV.File(joinpath(rgi_path, "glathida.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"])
@@ -214,10 +214,10 @@ function get_glathida_glacier(glacier::Glacier2D, params::Parameters, force)
         end
 
         gtd_grid .= ifelse.(count > 0, gtd_grid ./ count, 0.0)
-        
-        # Save file 
-        h5open(joinpath(gdir.dir, "glathida.h5"), "w") do file
-            write(file, "gtd_grid", gtd_grid)  
+
+        # Save file
+        h5open(joinpath(params.simulation.rgi_paths[glacier.rgi_id], "glathida.h5"), "w") do file
+            write(file, "gtd_grid", gtd_grid)
         end
     end
     return gtd_grid
@@ -277,6 +277,7 @@ function filter_missing_glaciers!(rgi_ids::Vector{String}, params::Parameters) #
 
     # Check which glaciers we can actually process
     rgi_stats = pd[].read_csv(utils[].file_downloader("https://cluster.klima.uni-bremen.de/~oggm/rgi/rgi62_stats.csv"), index_col=0)
+    TODO: update here
     # rgi_stats = rgi_stats.loc[rgi_ids]
 
     # if any(rgi_stats.Connect .== 2)

src/parameters/Parameters.jl

@@ -4,23 +4,23 @@ abstract type AbstractParameters end
 
 const AbstractEmptyParams = Union{AbstractParameters,Nothing}
 
-struct Parameters{PPHY <: AbstractEmptyParams, PSIM <: AbstractEmptyParams, PHY <: AbstractEmptyParams, 
-        PSOL <: AbstractEmptyParams, PUDE <: AbstractEmptyParams, POGGM <: AbstractEmptyParams, PINV <: AbstractEmptyParams}  
+struct Parameters{PPHY <: AbstractEmptyParams, PSIM <: AbstractEmptyParams, PHY <: AbstractEmptyParams,
+        PSOL <: AbstractEmptyParams, PUDE <: AbstractEmptyParams, PINV <: AbstractEmptyParams}
         physical::PPHY
         simulation::PSIM
         hyper::PHY
         solver::PSOL
         UDE::PUDE
-        inversion::PINV  
+        inversion::PINV
 end
 
 include("PhysicalParameters.jl")
 include("SimulationParameters.jl")
 
 """
 Parameters(;
+        physical::PhysicalParameters = PhysicalParameters(),
         simulation::SimulationParameters = SimulationParameters()
-        physical::PhysicalParameters = PhysicalParameters()
         )
 Initialize ODINN parameters
 
@@ -35,7 +35,7 @@ function Parameters(;
 
     # Build the parameters based on all the subtypes of parameters
     parameters = Parameters(physical, simulation,
-                            nothing, nothing, nothing, nothing)  
+                            nothing, nothing, nothing, nothing)
 
     if parameters.simulation.multiprocessing
             enable_multiprocessing(parameters.simulation.workers)

src/parameters/SimulationParameters.jl

@@ -4,17 +4,18 @@ struct SimulationParameters{I <: Integer, F <: AbstractFloat} <: AbstractParamet
     use_MB::Bool
     use_iceflow::Bool
     plots::Bool
-    velocities::Bool 
+    velocities::Bool
     overwrite_climate::Bool
     use_glathida_data::Bool
     float_type::DataType
     int_type::DataType
     tspan::Tuple{F, F}
     step::F
     multiprocessing::Bool
-    workers::I 
+    workers::I
     working_dir::String
     test_mode::Bool
+    rgi_paths::Dict
 end
 
 
@@ -29,8 +30,12 @@ end
                         float_type::DataType = Float64,
                         int_type::DataType = Int64,
                         tspan::Tuple{F, F} = (2010.0,2015.0),
+                        step::F = 1/12,
                         multiprocessing::Bool = true,
-                        workers::I = 4
+                        workers::I = 4,
+                        working_dir::String = "",
+                        test_mode::Bool = false,
+                        rgi_paths::Dict{String, String} = Dict()
         )
 Initialize the parameters for a simulation.
 Keyword arguments
@@ -55,7 +60,7 @@ function SimulationParameters(;
             workers::I = 4,
             working_dir::String = "",
             test_mode::Bool = false,
-            rgi_paths::Dict{String, String} = Dict()
+            rgi_paths::Dict = Dict() # TODO: once the mismatch in the prototype is fixed, revert to rgi_paths::Dict{String, String} = Dict()
             ) where {I <: Integer, F <: AbstractFloat}
 
     simulation_parameters = SimulationParameters(use_MB, use_iceflow, plots, velocities,

src/setup/helper_utilities.jl

@@ -11,21 +11,20 @@ function safe_approx(a, b)
     end
 end
 
-# Function for Python objects
-function safe_getproperty(obj::Py, prop_name::Symbol)
-    if PyCall.hasproperty(obj, prop_name)
-        return PyCall.getproperty(obj, prop_name)
-    else
-        return 0.0
-    end
-end
+# # Function for Python objects
+# function safe_getproperty(obj::Py, prop_name::Symbol)
+#     if PyCall.hasproperty(obj, prop_name)
+#         return PyCall.getproperty(obj, prop_name)
+#     else
+#         return 0.0
+#     end
+# end
 
-# Function for Julia objects
-function safe_getproperty(obj, prop_name::Symbol)
-    if hasproperty(obj, prop_name)
-        return getproperty(obj, prop_name)
-    else
-        return 0.0
-    end
-end
-
+# # Function for Julia objects
+# function safe_getproperty(obj, prop_name::Symbol)
+#     if hasproperty(obj, prop_name)
+#         return getproperty(obj, prop_name)
+#     else
+#         return 0.0
+#     end
+# end