Sleipnir.jl API

AbstractData

Abstract type that represents data. Used to implement ThicknessData, SurfaceVelocityData and DhdtData.

AbstractGlacier

An abstract type representing a glacier. This serves as a base type for different glacier implementations in the Sleipnir package.

AbstractInput

Abstract type representing an input source for a Law. Concrete subtypes must implement:

• default_name(::ConcreteInput): returns the default field name (as a Symbol) under which this input will be accessed in the law's input NamedTuple.
• get_input(::ConcreteInput, simulation, glacier_idx, t): retrieves the actual input value given the simulation context, glacier index, and time t.

AbstractLaw

Abstract type representing a synthetic law. Currently it's only used for testing by making easier to create dumb laws, but in the future it may be cleaner to use different concrete type of laws (for example CallbackLaw, ContinuousLaw, or LearnableLaw)

Concrete subtypes must implement:

• apply_law!(::ConcreteLaw, state, simulation, glacier_idx, t, θ)
• init_cache(::ConcreteLaw, glacier, glacier_idx)
• law_VJP_input(::ConcreteLaw, cache, simulation, glacier_idx, t, θ)
• law_VJP_θ(::ConcreteLaw, cache, simulation, glacier_idx, t, θ)
• precompute_law_VJP(::ConcreteLaw, cache, vjpsPrepLaw, simulation, glacier_idx, t, θ)
• cache_type(::ConcreteLaw)
• is_callback_law(::ConcreteLaw)
• is_precomputable_law_VJP(::ConcreteLaw)
• callback_freq(::ConcreteLaw)
• inputs(::ConcreteLaw)
• inputs_defined(::ConcreteLaw)
• apply_law_in_model(::ConcreteLaw)

AbstractModel

An abstract type that serves as a base for all model types in ODINN.

AbstractParameters

An abstract type that serves as a base for all parameter-related types in the ODINN ecosystem.

AbstractPrepVJP

Abstract type representing the preparation of Vector-Jacobian Product (VJP) computations for the laws. Subtypes of AbstractPrepVJP are used to handle any precomputations or setup required before evaluating VJPs, such as configuring automatic differentiation backends or precompiling code. This type provides a flexible interface for implementing custom VJP preparation strategies for different laws. It is for internal use only and not exposed to the user.

Cache

Abstract type for defining a cache struct of a law.

Mandatory field:

• value: Store the result of a forward evaluation.

Optional fields:

• vjp_inp: Store the result of the evaluation of the vector-Jacobian product (VJP) with respect to the inputs.
• vjp_θ: Store the result of the evaluation of the vector-Jacobian product (VJP) with respect to the parameters θ.

Notes:

• If a concrete subtype does not implement vjp_inp and vjp_θ, the law should not be used for gradient computation, and therefore for inversions.

A mutable struct representing a 2D climate for a glacier with various buffers and datasets.

Climate2D{CLIMRAW <: RasterStack, CLIMRAWSTEP <: RasterStack, CLIMSTEP <: ClimateStep, CLIM2DSTEP <: Climate2Dstep, F <: AbstractFloat}

Fields

• raw_climate::CLIMRAW: Raw climate dataset for the whole simulation.

• climate_raw_step::CLIMRAWSTEP: Raw climate trimmed for the current step to avoid memory allocations.

• climate_step::ClimateStep: Climate data for the current step.

• climate_2D_step::Climate2Dstep: 2D climate data for the current step to feed to the mass balance (MB) model.

• longterm_temps::Vector{F}: Long-term temperatures for the ice rheology.

• avg_temps::F: Intermediate buffer for computing average temperatures.

• avg_gradients::F: Intermediate buffer for computing average gradients.

• ref_hgt::F: Reference elevation of the raw climate data.

  Climate2D( rgi_id, params::Parameters, S::Matrix{<: AbstractFloat}, Coords::Dict, )

Initialize the climate data given a RGI ID, a matrix of surface elevation and glacier coordinates.

Arguments

• rgi_id: The glacier RGI ID.
• params::Parameters: The parameters containing simulation settings and paths.
• S::Matrix{<: AbstractFloat}: Matrix of surface elevation used to initialize the downscaled climate data.
• Coords::Dict: Coordinates of the glacier.

Description

This function initializes the climate data for a glacier by:

1. Creating a dummy period based on the simulation time span and step.

2. Loading the raw climate data from a NetCDF file.

3. Calculating the cumulative climate data for the dummy period.

4. Downscaling the cumulative climate data to a 2D grid.

5. Retrieving long-term temperature data for the glacier.

6. Returning the climate data, including raw climate data, cumulative climate data, downscaled 2D climate data, long-term temperatures, average temperatures, and average gradients.

   Climate2D( rawclimate::RasterStack, climaterawstep::RasterStack, climatestep::ClimateStep, climate2Dstep::Climate2Dstep, longtermtemps::Vector{<: AbstractFloat}, avgtemps::AbstractFloat, avggradients::AbstractFloat, refhgt::AbstractFloat, )

Initialize the climate data with the fields provided as arguments. Refer to the list of fields for a complete description of the arguments.

Climate2Dstep{F <: AbstractFloat}

A mutable struct representing a 2D climate time step with various climate-related parameters.

Keyword arguments

• temp::Matrix{F}: Temperature matrix.

• PDD::Matrix{F}: Positive Degree Days matrix.

• snow::Matrix{F}: Snowfall matrix.

• rain::Matrix{F}: Rainfall matrix.

• elevation_diff::Matrix{F}: Elevation difference matrix.

• aspect::Matrix{F}: Surface aspect matrix in degrees.

• albedo::Matrix{F}: Albedo matrix.

• slhf::Matrix{F}: Surface latent heat flux matrix.

• slope::Matrix{F}: Surface slope matrix in degrees.

• sshf::Matrix{F}: Surface sensible heat flux matrix.

• ssrd::Matrix{F}: Surface shortwave radiation downwards matrix.

• str::Matrix{F}: Surface net thermal radiation matrix.

• gradient::F: Altitudinal gradient value.

• avg_gradient::F: Average gradient value.

• x::Vector{F}: X-coordinates vector.

• y::Vector{F}: Y-coordinates vector.

• ref_hgt::F: Reference height.

ClimateStep{F <: AbstractFloat}

Mutable struct that represents a climate step before downscaling.

Keyword arguments

• prcp::F: Cumulative precipitation for the given period.
• temp::F: Cumulative temperature at the reference elevation for the given period.
• gradient::F: Cumulative temperature gradient for the given period.
• avg_temp::F: Average temperature over the time step.
• avg_gradient::F: Average temperature gradient over the time step.
• ref_hgt::F: Reference elevation of the raw climate data.

ConstantLaw{CACHE_TYPE}(init_cache)

Creates a constant law of type ConstantLaw{CACHE_TYPE} that holds a fixed value for the entire simulation.

This is useful to inject glacier-specific or global constants into the simulation without modifying them over time. The update function is a no-op, and only the init_cache function matters.

Arguments

• init_cache::Function: A function init_cache(simulation, glacier_idx, θ)::T that provides the constant value.

Type Parameters

• CACHE_TYPE: The type of the cache. Must be specified manually and should match the return type of init_cache.

Examples

# Same value for all glaciers
n_law = ConstantLaw{Float64}(Returns(4.))

# Value depending on the glacier
n_law = ConstantLaw{Float64}((sim, i, θ) -> sim.glaciers[i].n)

# Learned value
n_law = ConstantLaw{Float64}((sim, i, θ) -> θ.n)

Container

Abstract type that defines a container to be used in the PDE solver. It is useful to retrieve the simulation object when applying callback laws.

CustomVJP <: VJPType

Indicates that a law uses a custom-defined function for VJP computation. This is used when the VJP is provided manually rather than computed automatically.

DIVJP <: VJPType

Indicates that a law uses the default VJP computation provided by DifferentiationInterface.jl. This is used when no custom VJP function is provided, and the Vector-Jacobian Product (VJP) is computed automatically using DifferentiationInterface.jl.

Simple snapshot of mean glacier surface elevation change. This represents the glacier-wide average dh/dt. The convention is that negative glacier-wide MB is represented as a negative dh/dt (we assume that the ice density is constant).

Note: if you use an external dataset please make sure that you use the correct unit (surface elevation change vs m.w.e.)

GenInputsAndApply{IN, F}

Given a tuple of AbstractInputs and a function f, returns a callable struct. This struct with_input_f can be evaluated as a function that generates the inputs and then applies the function's law with_input_f.f. It is for internal use only and it isn't exposed to the user.

A mutable struct representing a 2D glacier. Notice that all fields can be empty by providing nothing as the default value.

/!\ WARNING /!\ Glacier objects should not be constructed manually, but rather through the initialize_glaciers function.

Glacier2D{F <: AbstractFloat, I <: Integer, CLIM <: Climate2D, THICKDATA <: Union{<:ThicknessData, Nothing}, SURFVELDATA <: Union{<:SurfaceVelocityData, Nothing}, DHDTDATA <: Union{<:DhdtData, Nothing}} <: AbstractGlacier

Fields

• rgi_id::String: The RGI (Randolph Glacier Inventory) identifier for the glacier.
• name::String: The name of the glacier if available.
• climate::CLIM: The climate data associated with the glacier.
• H₀::Matrix{F}: Initial ice thickness matrix.
• H_glathida::Matrix{F}: Ice thickness matrix from the GLATHIDA dataset.
• S::Matrix{F}: Surface elevation matrix.
• B::Matrix{F}: Bedrock elevation matrix.
• V::Matrix{F}: Ice velocity magnitude matrix.
• Vx::Matrix{F}: Ice velocity in the x-direction matrix.
• Vy::Matrix{F}: Ice velocity in the y-direction matrix.
• A::F: Flow law parameter.
• C::F: Sliding law parameter.
• n::F: Flow law exponent.
• p::F: Power law exponent associated to Weertman sliding law (Power associated to basal drag).
• q::F: Power law exponent associated to Weertman sliding law (Power associated to normal pressure).
• slope::Matrix{F}: Surface slope matrix.
• dist_border::Matrix{F}: Distance to the glacier border matrix.
• mask::BitMatrix: Boolean matrix representing the glacier mask, where true values indicate regions constrained by the mask (i.e., no-ice zones)
• mask_loss::BitMatrix: Boolean matrix representing mask for inversion. Losses used for inversions will only be evaluated withing the mask.
• Coords::Dict{String, Vector{Float64}}: Coordinates dictionary with keys as coordinate names and values as vectors of coordinates.
• Δx::F: Grid spacing in the x-direction.
• Δy::F: Grid spacing in the y-direction.
• nx::I: Number of grid points in the x-direction.
• ny::I: Number of grid points in the y-direction.
• cenlon::F: Longitude of the glacier center.
• cenlat::F: Latitude of the glacier center.
• params_projection::Dict{String, Float64}: Projection parameters that allows mapping the regional grid to global WGS84 coordinates.
• thicknessData::THICKDATA: Thickness data structure that is used to store the reference values.
• velocityData::SURFVELDATA: Surface velocity data structure that is used to store the reference values.
• dhdtData::DHDTDATA: Structure that is used to store the reference values of the mean glacier surface elevation change. The uncertainty field on DhdtData carries the associated uncertainty (e.g. from Hugonnet et al. 2021).

Glacier2D(
    glacier::Glacier2D;
    thicknessData::Union{<: ThicknessData, Nothing} = nothing,
    velocityData::Union{<: SurfaceVelocityData, Nothing} = nothing,
    dhdtData::Union{<: DhdtData, Nothing} = nothing,
)

Copies a Glacier2D object and updates the thickness and/or surface velocity data.

Arguments

• glacier::Glacier2D: The original glacier struct.
• thicknessData::Union{<: ThicknessData, Nothing}: Thickness data structure that is used to store the reference values. Default is nothing which keeps the existing thickness data.
• velocityData::Union{<: SurfaceVelocityData, Nothing}: Surface velocity data structure that is used to store the reference values. Default is nothing which keeps the existing surface velocity data.
• dhdtData::Union{<: DhdtData, Nothing}: Structure that is used to store the reference values of the mean glacier surface elevation change. Default is nothing which keeps the existing mean glacier surface elevation change.

Returns

• A Glacier2D object that is a copy of the original one with the thickness, surface velocity and/or elevation change data updated.

Glacier2D(
    rgi_id::String,
    params::Parameters;
    masking::Union{Int, Nothing, BitMatrix} = 2,
    smoothing=false
)

Build glacier object for a given RGI ID and parameters.

Arguments

• rgi_id::String: The RGI ID of the glacier.

• params::Parameters: A Parameters object containing simulation parameters.

• masking::Union{Int, Nothing, BitMatrix}: Type of mask applied to the glacier to determine regions with no ice.

  • When masking is an Int, the mask is based on the initial ice thickness H₀ and it is set to true for pixels outside at a distance of the glacier borders greater than the value of masking.
  • When masking is set to nothing, the mask is set to a BitMatrix full of falses.
  • When masking is a BitMatrix, this matrix is used for the mask. Defaults to 2.

• smoothing::Bool=false: Optional; whether to apply smoothing to the initial ice thickness. Default is false.

• test::Bool=false: Optional; test flag. Default is false.

Returns

• glacier::Glacier2D: A Glacier2D object initialized with the glacier data.

Description

This function loads and initializes the glacier data for a given RGI ID. It retrieves the initial ice thickness conditions based on the specified source in the parameters, applies optional smoothing, and initializes the glacier's topographical and velocity data. The function also handles Mercator projection for the glacier coordinates and filters glacier borders in high elevations to avoid overflow problems.

Notes

• The function reverses the matrices for ice thickness, bedrock, and other data to match the required orientation.
• If the Mercator projection includes latitudes larger than 80°, a warning is issued.
• If the glacier data is missing, the function updates a list of missing glaciers and issues a warning.

Glacier2D(;
    rgi_id::String = "",
    name::String = "",
    climate::Union{Climate2D, Nothing} = nothing,
    H₀::Matrix{F} = Matrix{Sleipnir.Float}([;;]),
    H_glathida::Matrix{F} = Matrix{Sleipnir.Float}([;;]),
    S::Matrix{F} = Matrix{Sleipnir.Float}([;;]),
    B::Matrix{F} = Matrix{Sleipnir.Float}([;;]),
    V::Matrix{F} = Matrix{Sleipnir.Float}([;;]),
    Vx::Matrix{F} = Matrix{Sleipnir.Float}([;;]),
    Vy::Matrix{F} = Matrix{Sleipnir.Float}([;;]),
    A::F = 0.0,
    C::F = 0.0,
    n::F = 0.0,
    p::F = 0.0,
    q::F = 0.0,
    slope::Matrix{F} = Matrix{Sleipnir.Float}([;;]),
    dist_border::Matrix{F} = Matrix{Sleipnir.Float}([;;]),
    mask::BitMatrix = BitMatrix([;;]),
    mask_loss::BitMatrix = BitMatrix([;;]),
    Coords::Dict{String, Vector{Float64}} = Dict{String, Vector{Float64}}("lon" => [], "lat" => []),
    Δx::F = 0.0,
    Δy::F = 0.0,
    nx::I = 0,
    ny::I = 0,
    cenlon::F = NaN,
    cenlat::F = NaN,
    params_projection::Dict{String, Float64} = Dict{String, Float64}(),
    thicknessData::THICKDATA = nothing,
    velocityData::SURFVELDATA = nothing,
    dhdtData::DHDTDATA = nothing,
) where {
    F <: AbstractFloat,
    I <: Integer,
    THICKDATA <: Union{<: ThicknessData, Nothing},
    SURFVELDATA <: Union{<: SurfaceVelocityData, Nothing},
    DHDTDATA <: Union{<:DhdtData, Nothing},
}

Constructs a Glacier2D object with the given parameters, including default ones.

Arguments

• rgi_id::String: The RGI identifier for the glacier.
• name::String: The name of the glacier if available.
• climate::Union{Climate2D, Nothing}: The climate data associated with the glacier. Defaults to nothing which falls back to DummyClimate2D().
• H₀::Matrix{F}: Initial ice thickness matrix.
• H_glathida::Matrix{F}: Ice thickness matrix from GLATHIDA.
• S::Matrix{F}: Surface elevation matrix.
• B::Matrix{F}: Bed elevation matrix.
• V::Matrix{F}: Ice velocity magnitude matrix.
• Vx::Matrix{F}: Ice velocity in the x-direction matrix.
• Vy::Matrix{F}: Ice velocity in the y-direction matrix.
• A::F: Flow law parameter.
• C::F: Sliding law parameter.
• n::F: Flow law exponent.
• p::F: Power law exponent associated to Weertman sliding law (Power associated to basal drag).
• q::F: Power law exponent associated to Weertman sliding law (Power associated to normal pressure).
• slope::Matrix{F}: Slope matrix.
• dist_border::Matrix{F}: Distance to border matrix.
• mask::BitMatrix: Boolean matrix representing the glacier mask, where true values indicate regions constrained by the mask (i.e., no-ice zones)
• mask_loss::BitMatrix: Boolean matrix representing mask for inversion. Losses used for inversions will only be evaluated withing the mask.
• Coords::Dict{String, Vector{Float64}}: Coordinates dictionary with keys "lon" and "lat".
• Δx::F: Grid spacing in the x-direction.
• Δy::F: Grid spacing in the y-direction.
• nx::I: Number of grid points in the x-direction.
• ny::I: Number of grid points in the y-direction.
• cenlon::F: Central longitude of the glacier.
• cenlat::F: Central latitude of the glacier.
• params_projection::Dict{String, Float64}: Projection parameters that allows mapping the regional grid to global WGS84 coordinates.
• thicknessData::THICKDATA: Thickness data structure that is used to store the reference values.
• velocityData::SURFVELDATA: Surface velocity data structure that is used to store the reference values.
• dhdtData::DHDTDATA: Structure that is used to store the reference values of the mean glacier surface elevation change.

Returns

• A Glacier2D object with the specified parameters.

IntegratedTrajectoryMapping <: VelocityMapping

Integrated trajectory mapping. This mapping is closer to reality as it consists in integrating over time the instantaneous ice surface velocities along ice flow trajectories in a Lagrangian way. This integrated velocity is then compared to the velocity of the datacube. It has not been implemented yet but its computational cost will likely be expensive.

Fields

• spatialInterp::Symbol: The spatial interpolation to use to map the ice surface velocity grid to the glacier grid. For the moment only :nearest is supported.

Law{T}(;
    inputs = nothing,
    f!, f_VJP_input! = nothing, f_VJP_θ! = nothing,
    init_cache, callback_freq = nothing,
    p_VJP! = nothing,
    max_value = NaN, min_value = NaN, name = :unknown,
) where{T}

Defines a physical or empirical law applied to a glacier model that mutates an internal state T at each simulation time step.

# ❌ Will not work: Float64 is immutable, so cache .= ... has no effect
Law{Float64}(;
    f! = (cache, _, _, t, θ) -> cache = θ.scale * sin(2π * t + θ.shift),
    init_cache = (_, _, _) -> 0.0,
)

# ✅ Correct: using a 0-dimensional array allows in-place mutation
Law{Array{Float64, 0}}(;
    f! = (cache, _, _, t, θ) -> cache .= θ.scale * sin(2π * t + θ.shift) + θ.bias,
    init_cache = (_, _, _) -> zeros(),
)

Arguments

• f!::Function: A function with signature f!(cache::T, simulation, glacier_idx, t, θ) that updates the internal state. If inputs are provided, the function instead takes the form f!(cache::T, inputs, θ).
• init_cache::Function: A function init_cache(simulation, glacier_idx, θ)::T that initializes the internal state for a given glacier.
• callback_freq::Union{Nothing, Real, Period, Month}: Optional. If provided, the law is treated as a callback law and is only applied every callback_freq time units. If callback_freq is set to zero, then the law is applied only once at the beginning of the simulation. If callback_freq is set to nothing (default), then the law is applied at every iteration. If callback_freq is provided as a Period or Month, it is converted to a float value in a yearly basis.
• f_VJP_input!: A function with signature (cache::T, simulation, glacier_idx, t, θ) that updates cache.vjp_inp which is the VJP with respect to the inputs.
• f_VJP_θ!: A function with signature (cache::T, simulation, glacier_idx, t, θ) that updates cache.vjp_θ which is the VJP with respect to the parameters θ.
• p_VJP!: A function with signature (cache::T, vjpsPrepLaw, simulation, glacier_idx, t, θ) that performs the precomputation of the VJPs.
• inputs::Union{Nothing, Tuple{<:AbstractInput}}: Optional. Provides automatically generated inputs passed to f! at runtime.
• max_value::Float64: Optional. The maximum value that the law can take, used for plotting and capping the function output.
• min_value::Float64: Optional. The minimum value that the law can take, used for plotting and capping the function output.
• name::Symbol: A name for the law, used for identification and plotting.

Type Parameters

• T: The type of the internal state. Must be specified manually and should match the return type of init_cache.

Notes

• Refer to the tutorials in the documentation for a complete description of the VJP options.

Examples

# A law applied at every timestep, storing a scalar value
Law{Array{Float64, 0}}(;
    f! = (cache, _, _, t, θ) -> cache .= θ.scale * sin(2π * t + θ.shift) + θ.bias,
    init_cache = (_, _, _) -> zeros(),
)

# A callback law applied once per month (assuming time in years)
Law{Array{Float64, 0}}(;
    f! = (cache, _, _, t, θ) -> cache .= θ.scale * sin(2π * t + θ.shift) + θ.bias,
    init_cache = (_, _, _) -> zeros(),
    callback_freq = 1 / 12,
)

MatrixCache <: Cache

A cache structure for storing a two-dimensional array of Float64 values along with their associated vector-Jacobian products (VJP). This is typically used for spatially varying laws. Fields:

• value::Array{Float64, 2}: The cached matrix.
• vjp_inp::Array{Float64, 2}: VJP with respect to inputs.
• vjp_θ::Vector{Float64}: VJP with respect to parameters.

MatrixCacheNoVJP <: Cache

A mutable cache structure for storing a two-dimensional array of Float64 values. This is typically used for spatially varying laws. This struct is intended for use cases where the law is not differentiated, and hence the vector-Jacobian products (VJP) are not required. Fields:

• value::Array{Float64, 2}: The cached matrix.

MeanDateVelocityMapping <: VelocityMapping

Mean date velocity mapping. It is the most simple mapping one can build and it consists in taking the 2D vector field of ice velocity associated to a given mean date and compare it to the instantaneous ice surface velocity obtained from the ice flow model. It is valid only for ice surface velocities estimated from short time windows since the velocity can vary within this time window.

Fields

• spatialInterp::Symbol: The spatial interpolation to use to map the ice surface velocity grid to the glacier grid. For the moment only :nearest is supported.

Model{IFM <: AbstractEmptyModel, MBM <: Union{<:AbstractEmptyModel, Vector{<:AbstractEmptyModel}}, TC <: AbstractEmptyModel}

A mutable struct that represents a model with three components: iceflow, mass balance, and machine learning.

Model(
    iceflow::IFM,
    mass_balance::MBM,
    trainable_components::TC,
) where {IFM <: AbstractEmptyModel, MBM <: Union{<:AbstractEmptyModel, Vector{<:AbstractEmptyModel}}, TC <: AbstractEmptyModel}

Model(;iceflow, mass_balance) = Model(iceflow, mass_balance, nothing)

Initialize Model (no machine learning model).

Keyword arguments

• iceflow::IFM: Represents the iceflow component, which is an instance of IFM.
• mass_balance::MBM: Represents the mass balance component. Either a single MBM instance or a Vector{MBM} of per-glacier models.
• trainable_components::TC: Represents the trainable components, which is an instance of TC.

Type Parameters

• IFM: A subtype of AbstractEmptyModel representing the type of the iceflow model.
• MBM: Either a subtype of AbstractEmptyModel (single shared model) or a Vector of such subtypes (per-glacier models). The field stores the exact runtime type — no Union in the field — which is required for AD compatibility.
• TC: A subtype of AbstractEmptyModel representing the type of the trainable components.

Model(;
    iceflow::Union{<: AbstractModel, Nothing} = nothing,
    mass_balance::Union{<: AbstractModel, Nothing} = nothing,
    regressors::Union{NamedTuple, Nothing} = nothing
)

Creates a new model instance using the provided iceflow, mass balance, and machine learning components.

Arguments

• iceflow::Union{IFM, Nothing}: The iceflow model to be used. Can be a single model or nothing.
• mass_balance::Union{MBM, Nothing}: The mass balance model to be used. Can be a single model or nothing.
• regressors::Union{NamedTuple, Nothing}: The regressors to be used in the laws.

Returns

• model: A new instance of Sleipnir.Model initialized with the provided components.

Note

Since Model is a keyword argument function, it is defined in Sleipnir. Otherwise we would have to define it multiple times, making the ODINN precompilation impossible because of method overwriting.

ModelCache{IFC, MBC}

Cache struct that holds the internal state or memory buffers for the components of a Model.

Typically used to store per-glacier preallocated buffers or intermediate results that persist across time steps during simulation.

Fields

• iceflow::IFC: Cache associated with the iceflow model.
• mass_balance::MBC: Cache associated with the mass balance model.

Type Parameters

• IFC: Cache type for the iceflow model.
• MBC: Cache type for the mass balance model.

NullLaw <: AbstractLaw

This struct represents a law that is not used in the iceflow model.

mutable struct Parameters{PPHY <: AbstractEmptyParams, PSIM <: AbstractEmptyParams, PHY <: AbstractEmptyParams,
    PSOL <: AbstractEmptyParams, PUDE <: AbstractEmptyParams}

A mutable struct that holds various parameter sets for different aspects of a simulation or model.

Fields

• physical::PPHY: Physical parameters.
• simulation::PSIM: Simulation parameters.
• hyper::PHY: Hyperparameters.
• solver::PSOL: Solver parameters.
• UDE::PUDE: Universal Differential Equation (UDE) parameters.

Type Parameters

• PPHY: Type of the physical parameters, must be a subtype of AbstractEmptyParams.
• PSIM: Type of the simulation parameters, must be a subtype of AbstractEmptyParams.
• PHY: Type of the hyperparameters, must be a subtype of AbstractEmptyParams.
• PSOL: Type of the solver parameters, must be a subtype of AbstractEmptyParams.
• PUDE: Type of the UDE parameters, must be a subtype of AbstractEmptyParams.

Parameters(; physical::PhysicalParameters = PhysicalParameters(), simulation::SimulationParameters = SimulationParameters())

Constructs a Parameters object with the given physical and simulation parameters.

Arguments

• physical::PhysicalParameters: An instance of PhysicalParameters (default: PhysicalParameters()).
• simulation::SimulationParameters: An instance of SimulationParameters (default: SimulationParameters()).

Returns

• A Parameters object initialized with the provided physical and simulation parameters.

Notes

• If simulation.multiprocessing is enabled, multiprocessing is configured with the specified number of workers.

A structure representing physical parameters used in simulations.

PhysicalParameters{F <: AbstractFloat}

Fields

• ρ::F: Density of ice.
• g::F: Gravitational acceleration.
• ϵ::F: Regularization used in the square root of norms for AD numerical stability.
• η₀::F: Initial viscosity.
• maxA::F: Maximum A.
• minA::F: Minimum A.
• maxC::F: Maximum C.
• minC::F: Minimum C.
• maxTlaw::F: Maximum temperature according to some law.
• minTlaw::F: Minimum temperature according to some law.
• noise_A_magnitude::F: Magnitude of noise in A.
• ρ_w::F: Density of water (kg m⁻³), used for ice-to-water-equivalent conversions.
• DDF_min::F: Minimum degree-day factor for TI model calibration (m w.e. °C⁻¹ d⁻¹).
• DDF_max::F: Maximum degree-day factor for TI model calibration (m w.e. °C⁻¹ d⁻¹).
• prcp_fac_min::F: Minimum precipitation correction factor for TI model calibration.
• prcp_fac_max::F: Maximum precipitation correction factor for TI model calibration.
• temp_bias_min::F: Minimum temperature bias (°C) for TI model calibration. Default: -10.0.
• temp_bias_max::F: Maximum temperature bias (°C) for TI model calibration. Default: 10.0.

Initialize the physical parameters of a model.

PhysicalParameters(;
    ρ::Float64 = 900.0,
    g::Float64 = 9.81,
    ϵ::Float64 = 1e-10,
    η₀::F = 1.0,
    maxA::Float64 = 8e-17,
    minA::Float64 = 8.5e-20,
    maxC::Float64 = 8e-17, # TODO: to be revised
    minC::Float64 = 8.5e-20,
    maxTlaw::Float64 = 1.0,
    minTlaw::Float64 = -25.0,
    noise_A_magnitude::Float64 = 5e-18
    )

Keyword arguments

- `ρ`: Ice density
- `g`: Gravitational acceleration.
- `ϵ`: Regularization used in the square root of norms for AD numerical stability.
- `η₀`: Factor to cap surface elevation differences with the upstream ice thickness to impose boundary condition in the iceflow equation
- `maxA`: Maximum value for `A` (Glen's coefficient)
- `minA`: Minimum value for `A` (Glen's coefficient)
- `maxC`: Maximum value of sliding coefficient `C`
- `minC`: Minimum value of sliding coefficient `C`
- `maxTlaw`: Maximum value of Temperature used in simulations on fake law
- `minTlaw`: Minimum value of Temperature used in simulations on fake law
- `noise_A_magnitude`: Magnitude of noise added to A
- `ρ_w`: Water density (kg m⁻³). Default: 1000.0.
- `DDF_min`: Minimum degree-day factor for TI model calibration (m w.e. °C⁻¹ d⁻¹). Default: 0.5×10⁻³.
- `DDF_max`: Maximum degree-day factor for TI model calibration (m w.e. °C⁻¹ d⁻¹). Default: 20.0×10⁻³.
- `prcp_fac_min`: Minimum precipitation correction factor for TI model calibration. Default: 0.1.
- `prcp_fac_max`: Maximum precipitation correction factor for TI model calibration. Default: 10.0.
- `temp_bias_min`: Minimum temperature bias (°C) for TI model calibration. Default: -10.0.
- `temp_bias_max`: Maximum temperature bias (°C) for TI model calibration. Default: 10.0.

mutable struct Results{F <: AbstractFloat, I <: Integer}

A mutable struct to store the results of simulations.

Fields

• rgi_id::String: Identifier for the RGI (Randolph Glacier Inventory).
• H::Vector{Matrix{F}}: Vector of matrices representing glacier ice thickness H over time.
• H_glathida::Matrix{F}: Optional matrix for Glathida ice thicknesses.
• H_ref::Vector{Matrix{F}}: Reference data for ice thickness.
• S::Matrix{F}: Glacier surface altimetry.
• B::Matrix{F}: Glacier bedrock.
• V::Matrix{F}: Glacier ice surface velocities.
• Vx::Matrix{F}: x-component of the glacier ice surface velocity V.
• Vy::Matrix{F}: y-component of the glacier ice surface velocity V.
• V_ref::Matrix{F}: Reference data for glacier ice surface velocities V.
• Vx_ref::Matrix{F}: Reference data for the x-component of the glacier ice surface velocity Vx.
• Vy_ref::Matrix{F}: Reference data for the y-component of the glacier ice surface velocity Vy.
• date_Vref::Vector{F}: Date of velocity observation (mean of date1 and date2).
• date1_Vref::Vector{F}: First date of velocity acquisition.
• date2_Vref::Vector{F}: Second date of velocity acquisition.
• t_dhdt::Tuple{F, F}: Time window of the mean surface elevation change.
• dhdt_ref::F: Mean surface elevation change.
• Δx::F: Grid spacing in the x-direction.
• Δy::F: Grid spacing in the y-direction.
• lon::F: Longitude of the glacier grid center.
• lat::F: Latitude of the glacier grid center.
• nx::I: Number of grid points in the x-direction.
• ny::I: Number of grid points in the y-direction.
• tspan::Vector{F}: Time span of the simulation.

Results(glacier::G, ifm::IF;
    rgi_id::String = glacier.rgi_id,
    H::Vector{Matrix{F}} = Vector{Matrix{Sleipnir.Float}}([[;;]]),
    H_glathida::Matrix{F} = glacier.H_glathida,
    H_ref::Vector{Matrix{F}} = Vector{Matrix{Sleipnir.Float}}([[;;]]),
    S::Matrix{F} = zeros(Sleipnir.Float, size(ifm.S)),
    B::Matrix{F} = zeros(Sleipnir.Float, size(glacier.B)),
    V::Vector{Matrix{F}} = Vector{Matrix{Sleipnir.Float}}([[;;]]),
    Vx::Vector{Matrix{F}} = Vector{Matrix{Sleipnir.Float}}([[;;]]),
    Vy::Vector{Matrix{F}} = Vector{Matrix{Sleipnir.Float}}([[;;]]),
    V_ref::Vector{Matrix{F}} = Vector{Matrix{Sleipnir.Float}}([[;;]]),
    Vx_ref::Vector{Matrix{F}} = Vector{Matrix{Sleipnir.Float}}([[;;]]),
    Vy_ref::Vector{Matrix{F}} = Vector{Matrix{Sleipnir.Float}}([[;;]]),
    date_Vref::Vector{F} = Vector{Sleipnir.Float}([]),
    date1_Vref::Vector{F} = Vector{Sleipnir.Float}([]),
    date2_Vref::Vector{F} = Vector{Sleipnir.Float}([]),
    t_dhdt::Union{Tuple{F, F}, Nothing} = nothing,
    dhdt_ref::Union{F, Nothing} = nothing,
    Δx::F = glacier.Δx,
    Δy::F = glacier.Δy,
    lon::F = glacier.cenlon,
    lat::F = glacier.cenlat,
    nx::I = glacier.nx,
    ny::I = glacier.ny,
    t::Vector{F} = Vector{Sleipnir.Float}([]),
    tspan::Tuple{F, F} = (NaN, NaN),
) where {G <: AbstractGlacier, F <: AbstractFloat, IF <: AbstractModel, I <: Integer}

Construct a Results object for a glacier simulation.

Arguments

• glacier::G: The glacier object, subtype of AbstractGlacier.
• ifm::IF: The model object, subtype of AbstractModel.
• rgi_id::String: The RGI identifier for the glacier. Defaults to glacier.rgi_id.
• H::Vector{Matrix{F}}: Ice thickness matrices. Defaults to an empty vector.
• H_glathida::Matrix{F}: Ice thickness from GlaThiDa. Defaults to glacier.H_glathida.
• H_ref::Vector{Matrix{F}}: Reference ice thickness. Defaults to an empty vector.
• S::Matrix{F}: Surface elevation matrix. Defaults to a zero matrix of the same size as ifm.S.
• B::Matrix{F}: Bed elevation matrix. Defaults to a zero matrix of the same size as glacier.B.
• V::Vector{Matrix{F}}: Velocity magnitude matrix. Defaults to an empty vector.
• Vx::Vector{Matrix{F}}: Velocity in the x-direction matrix. Defaults to an empty vector.
• Vy::Vector{Matrix{F}}: Velocity in the y-direction matrix. Defaults to an empty vector.
• V_ref::Vector{Matrix{F}}: Reference velocity magnitude matrix. Defaults to an empty vector.
• Vx_ref::Vector{Matrix{F}}: Reference velocity in the x-direction matrix. Defaults to an empty vector.
• Vy_ref::Vector{Matrix{F}}: Reference velocity in the y-direction matrix. Defaults to an empty vector.
• date_Vref::Vector{F}: Date of velocity observation (mean of date1 and date2). Defaults to an empty vector.
• date1_Vref::Vector{F}: First date of velocity acquisition. Defaults to an empty vector.
• date2_Vref::Vector{F}: Second date of velocity acquisition. Defaults to an empty vector.
• t_dhdt::Union{Tuple{F, F}, Nothing}: Time window of the mean surface elevation change. Defaults to nothing in which case if glacier.dhdtData exists, glacier.dhdtData.t is used instead.
• dhdt_ref::Union{F, Nothing}: Mean surface elevation change. Defaults to nothing in which case glacier.dhdtData.dhdt is used when glacier.dhdtData is available.
• Δx::F: Grid spacing in the x-direction. Defaults to glacier.Δx.
• Δy::F: Grid spacing in the y-direction. Defaults to glacier.Δy.
• lon::F: Longitude of the glacier grid center. Defaults to glacier.cenlon.
• lat::F: Latitude of the glacier grid center. Defaults to glacier.cenlat.
• nx::I: Number of grid points in the x-direction. Defaults to glacier.nx.
• ny::I: Number of grid points in the y-direction. Defaults to glacier.ny.
• tspan::Tuple(F, F): Timespan of the simulation.
• θ::Union{Nothing, ComponentArray{F}}: Model parameters. Defaults to nothing.
• loss::Union{Nothing, Vector{F}}: Loss values. Defaults to nothing.

Returns

• results::Results: A Results object containing the simulation results.

ScalarCache <: Cache

A cache structure for storing a scalar value as a zero-dimensional array of Float64 along with their associated vector-Jacobian products (VJP). This is typically used for constant per glacier laws. Fields:

• value::Array{Float64, 0}: The cached scalar value.
• vjp_inp::Array{Float64, 0}: VJP with respect to inputs.
• vjp_θ::Vector{Float64}: VJP with respect to parameters.

ScalarCacheNoVJP <: Cache

A mutable cache structure for storing a scalar value as a zero-dimensional array of Float64. This is typically used for constant per glacier laws. This struct is intended for use cases where the law is not differentiated, and hence the vector-Jacobian products (VJP) are not required. Fields:

• value::Array{Float64, 0}: The cached scalar value.

Simulation

An abstract type representing a generic simulation. This type is intended to be subclassed by specific simulation types to provide a common interface and shared functionality for all simulations.

A structure to hold simulation parameters for a simulation in ODINN.

struct SimulationParameters{I <: Integer, F <: AbstractFloat, VM <: VelocityMapping} <: AbstractParameters

Fields

• use_MB::Bool: Flag to indicate whether mass balance should be used.
• calibrate_MB::Bool: Flag to indicate whether the mass balance model should be calibrated per glacier against geodetic observations when building a simulation. The calibration routine dispatches on the mass balance model type; model types without a calibration method are left unchanged.
• use_iceflow::Bool: Flag to indicate whether ice flow should be used.
• plots::Bool: Flag to indicate whether plots should be generated.
• use_velocities::Bool: Flag to indicate whether velocities should be calculated.
• f_surface_velocity_factor::F: Numerical factor representing the ratio between depth integrated ice velocity and surface velocity.
• overwrite_climate::Bool: Flag to indicate whether to overwrite climate data.
• use_glathida_data::Bool: Flag to indicate whether to use GLATHIDA data.
• tspan::Tuple{F, F}: Time span for the simulation.
• step_MB::F: Time step for the MB simulation.
• multiprocessing::Bool: Flag to indicate whether multiprocessing should be used.
• workers::I: Number of workers for multiprocessing.
• working_dir::String: Directory for working files.
• test_mode::Bool: Flag to indicate whether to run in test mode.
• rgi_paths::Dict{String, String}: Dictionary of RGI paths.
• ice_thickness_source::Symbol: Source of ice thickness data.
• velocity_product::Symbol: Source of velocity product data.
• mapping::VM: Mapping to use in order to grid the data from the coordinates of the velocity product datacube to the glacier grid.
• gridScalingFactor::I: Grid downscaling factor, used to speed-up the tests. Default value is 1 which means no downscaling is applied.
• catch_errors::Bool: Whether to catch errors during glacier initialization and write the glaciers to discard in a missing_glaciers.jld2.

Constructor for SimulationParameters type, including default values.

SimulationParameters(;
    use_MB::Bool = true,
    calibrate_MB::Bool = true,
    use_iceflow::Bool = true,
    plots::Bool = true,
    use_velocities::Bool = true,
    f_surface_velocity_factor::F = 1.0,
    overwrite_climate::Bool = false,
    use_glathida_data::Bool = false,
    tspan::Tuple{F, F} = (2010.0, 2015.0),
    step_MB::F = 1/12,
    multiprocessing::Bool = true,
    workers::I = 4,
    working_dir::String = "",
    test_mode::Bool = false,
    rgi_paths::Dict{String, String} = Dict{String, String}(),
    ice_thickness_source::Symbol = :Farinotti19,
    velocity_product::Symbol = :Millan22,
    climate_data_source::Symbol = :W5E5,
    mapping::VM = MeanDateVelocityMapping(),
    gridScalingFactor::I = 1,
    catch_errors::Bool = false
) where {I <: Integer, F <: AbstractFloat, VM <: VelocityMapping}

Keyword arguments

• use_MB::Bool: Whether to use mass balance (default: true).
• calibrate_MB::Bool: Whether to calibrate the mass balance model per glacier against geodetic observations when building a simulation (default: true).
• use_iceflow::Bool: Whether to use ice flow (default: true).
• plots::Bool: Whether to generate plots (default: true).
• use_velocities::Bool: Whether to calculate velocities (default: true).
• f_surface_velocity_factor::F: Numerical factor representing the ratio between depth integrated ice velocity and surface velocity (default: 1.0).
• overwrite_climate::Bool: Whether to overwrite climate data (default: false).
• use_glathida_data::Bool: Whether to use GLATHIDA data (default: false).
• velocity_product::Symbol: Source of velocity product data (default: :Millan22).
• float_type::DataType: Data type for floating point numbers (default: Float64).
• int_type::DataType: Data type for integers (default: Int64).
• tspan::Tuple{F, F}: Time span for the simulation (default: (2010.0, 2015.0)).
• step_MB::F: Time step for the MB simulation (default: 1/12).
• multiprocessing::Bool: Whether to use multiprocessing (default: true).
• workers::I: Number of workers for multiprocessing (default: 4).
• working_dir::String: Working directory for the simulation (default: "").
• test_mode::Bool: Whether to run in test mode (default: false).
• rgi_paths::Dict{String, String}: Dictionary of RGI paths (default: Dict{String, String}()).
• ice_thickness_source::Symbol: Source of ice thickness data, either :Millan22 or :Farinotti19 (default: :Farinotti19).
• climate_data_source::Symbol: Source of climate data (default: :W5E5).
• mapping::VM: Mapping to use in order to grid the data from the coordinates of the velocity product datacube to the glacier grid.
• gridScalingFactor::I: Grid downscaling factor, used to speed-up the tests. Default value is 1 which means no downscaling is applied.
• catch_errors::Bool: Whether to catch errors during glacier initialization and write the glaciers to discard in a missing_glaciers.jld2.

Returns

• simulation_parameters: A new SimulationParameters object.

Throws

• AssertionError: If ice_thickness_source is not :Millan22 or :Farinotti19.

Notes

• If the global variable ODINN_OVERWRITE_MULTI is set to true, multiprocessing is is enabled in any case and the number of workers specified in the simulation parameters must correspond to the number of processes with which Julia has been started. This is to allow the documentation to build successfully in ODINN as we cannot change the number of process in the CI.

A mutable struct representing a surface velocity data. Notice that all fields can be empty by providing nothing as the default value.

SurfaceVelocityData{F <: AbstractFloat} <: AbstractData

Fields

• x::Union{Vector{F}, Nothing}: Easting of observation.
• y::Union{Vector{F}, Nothing}: Northing of observation.
• lat::Union{Vector{F}, Nothing}: Latitude of observation.
• lon::Union{Vector{F}, Nothing}: Longitude of observation.
• vx::Union{Vector{Matrix{F}}, Nothing}: x / longitudinal component of surface velocity. Positive velocities correspond to the direction of increasing index of the glacier x-coordinate.
• vy::Union{Vector{Matrix{F}}, Nothing}: y / latitudinal component of surface velocity. Positive velocities correspond to the direction of increasing index of the glacier y-coordinate.
• vabs::Union{Vector{Matrix{F}}, Nothing}: Absolute ice surface velocity.
• vx_error::Union{Vector{F}, Nothing}: Error in vx
• vy_error::Union{Vector{F}, Nothing}: Error in vy
• vabs_error::Union{Vector{F}, Nothing}: Error in vabs.
• date::Union{Vector{DateTime}, Nothing}: Date of observation (mean of date1 and date2)
• date1::Union{Vector{DateTime}, Nothing}: First date of acquisition.
• date2::Union{Vector{DateTime}, Nothing}: Second date of acquisition.
• date_error::Union{Vector{Day}, Vector{Millisecond}, Nothing}: Error in date.
• flag::Union{BitMatrix, Nothing}: Flag indicating whether a pixel is considered as reliable or not.
• isGridGlacierAligned::Bool: Whether the data have been gridded to the glacier grid or not.

Constructs SurfaceVelocityData using data from Rabatel et. al (2023) with the given parameters, including default ones.

function SurfaceVelocityData(; x::Union{Vector{F}, Nothing} = nothing, y::Union{Vector{F}, Nothing} = nothing, lat::Union{Vector{F}, Nothing} = nothing, lon::Union{Vector{F}, Nothing} = nothing, vx::Union{Vector{Matrix{F}}, Nothing} = nothing, vy::Union{Vector{Matrix{F}}, Nothing} = nothing, vabs::Union{Vector{Matrix{F}}, Nothing} = nothing, vxerror::Union{Vector{F}, Nothing} = nothing, vyerror::Union{Vector{F}, Nothing} = nothing, vabserror::Union{Vector{F}, Nothing} = nothing, date::Union{Vector{DateTime}, Nothing} = nothing, date1::Union{Vector{DateTime}, Nothing} = nothing, date2::Union{Vector{DateTime}, Nothing} = nothing, dateerror::Union{Vector{Day}, Vector{Millisecond}, Nothing} = nothing, flag::Union{BitMatrix, Nothing} = nothing, isGridGlacierAligned::Bool = false, ) where {F <: AbstractFloat}

Constructor for ice surface velocity data based on Rabatel et. al (2023).

Important remarks:

• Velocities values are reported in m/yr. Positive velocities correspond to the direction of increasing index of the glacier. When the glacier is oriented in east-west and south-north (see latitude and coordinate ordering), positive velocities of the ice surface velocity correspond to positive east-west and south-north velocity component.
• The error in velocity is unique per timestamp, rather than being pixel distributed.
• The error in the absolute velocities vabs_error is overestimated.

References:

• Rabatel, A., Ducasse, E., Millan, R. & Mouginot, J. Satellite-Derived Annual Glacier Surface Flow Velocity Products for the European Alps, 2015–2021. Data 8, 66 (2023).

Simple time series of ice thickness data to test transient inversion

VJPType

Abstract type representing the mode of Vector-Jacobian Product (VJP) computation used in a law. Subtypes of VJPType define how the VJP is evaluated (e.g., custom or using DifferentiationInterface.jl).

VelocityMapping

Abstract type representing the mapping to use in order to map the ice velocity products onto the glacier grid. It contains all needed information to build both the spatial projection, and how to interpolate the data in time.

iTopoAspect{F<:AbstractFloat} <: AbstractInput

Input representing dynamic smoothed surface aspect (degrees, [0, 360)) computed from the current glacier surface topography.

iTopoSlope{F<:AbstractFloat} <: AbstractInput

Input representing dynamic surface slope (degrees) computed from the current glacier surface topography.

DummyClimate2D(;
    longterm_temps_scalar::Vector{F} = Vector{Sleipnir.Float}([]),
    longterm_temps_gridded::Matrix{F} = Matrix{Sleipnir.Float}(zeros(0, 0))
) where {F <: AbstractFloat}

Dummy climate initialization for very specific use cases where we don't have climate data and we need to build a minimalistic climate with only a few data. For the moment it supports only the initialization of the long term temperatures. It returns a minimalistic Climate2D instance.

Arguments:

• longterm_temps_scalar::Vector{F}: Scalar long term temperatures.
• longterm_temps_gridded::Matrix{F}: Distributed long term temperatures (matrix).

ReverseUTMercator(x::F, y::F; k=0.9996, cenlon=0.0, cenlat=0.0, x0=0.0, y0=0.0, zone::Union{Nothing, Int}=nothing, hemisphere=nothing) where {F <: AbstractFloat}

Transverse Mercator Projection. This function reprojects latitude/longitude into northing/easting coordinates.

Keyword arguments

- `k`: scale factor of the projection
- `cenlon`: Central longitude used in the projection
- `cenlat`: Central latitude used in the projection
- `x0`: Shift in easting
- `y0`: Shift in northing
- `zone` : Zone of the projection
- `hemisphere`: Either :north or :south

UTMercator(x::F, y::F; k=0.9996, cenlon=0.0, cenlat=0.0, x0=0.0, y0=0.0, zone::Union{Nothing, Int}=nothing, hemisphere=nothing) where {F <: AbstractFloat}

Transverse Mercator Projection. This function reprojects northing/easting coordinates into latitude/longitude.

Keyword arguments

- `k`: scale factor of the projection
- `cenlon`: Central longitude used in the projection
- `cenlat`: Central latitude used in the projection
- `x0`: Shift in easting
- `y0`: Shift in northing
- `zone` : Zone of the projection
- `hemisphere`: Either :north or :south

_decorate_geo_axis!(ax, nx, ny, lon, lat, Δx; scale_text_size, num_vars)

Add geographic labels (central lon/lat tick), axis styling, and a Δx-based scale bar to a glacier heatmap axis. Shared by heatmap, difference-evolution, and gridded-data plotting functions.

Return true when at least one ERA5-specific field in a ClimateStep/Climate2Dstep has a non-zero value — used by standalone show methods that lack parent context.

_overlay_contour!(ax, ctr)

Draw glacier-boundary contour lines on an existing axis. No-op when ctr is nothing.

_resolve_temporal_snapshot(data, timeIdx) -> Matrix

If data is a vector of matrices, pick element timeIdx (or the last one when timeIdx is nothing). Otherwise return data as-is.

_results_plot_metadata(results; plotContour) -> NamedTuple

Extract common metadata (lon, lat, x, y, rgi_id, Δx, mask, contour) from a Results object. The glacier mask is derived from the initial thickness field.

accumulate_gridded_data(
    gridded_data::Vector{Matrix{F}};
    weights::Union{Nothing,AbstractVector{<:Real}}=nothing,
) where {F <: AbstractFloat}

Accumulate a time series of gridded matrices into a single matrix.

Arguments

• gridded_data::Vector{Matrix{F}}: Sequence of gridded fields to accumulate.
• weights::Union{Nothing,AbstractVector{<:Real}}: Optional per-step weights. If provided, weighted accumulation is performed as sum(weights[i] * gridded_data[i]).

Returns

• Matrix{F}: The accumulated matrix.

apply_all_callback_laws!(model::AbstractModel, cache, simulation, glacier_idx, t, θ)

This function is a placeholder and must be implemented for your custom model type.

It is intended to apply all callback laws in the simulation to update the model cache for a given glacier at time t with parameters θ. By default, calling this function will throw an error to indicate that the user should provide their own implementation tailored to their model.

Arguments

• model::AbstractModel: The model instance.
• cache: The cache object storing state variables.
• simulation: The simulation context.
• glacier_idx: Index identifying the glacier.
• t: The current simulation time.
• θ: The parameter vector.

Throws

• Always throws an error: "This function should not be called. Implement apply_all_callback_laws! for your own model."

apply_all_non_callback_laws!(model::AbstractModel, cache, simulation, glacier_idx, t, θ)

This function is a placeholder and must be implemented for your custom model type.

It is intended to apply all non-callback laws in the simulation to update the model cache for a given glacier at time t with parameters θ. By default, calling this function will throw an error to indicate that the user should provide their own implementation tailored to their model.

Arguments

• model::AbstractModel: The model instance.
• cache: The cache object storing state variables.
• simulation: The simulation context.
• glacier_idx: Index identifying the glacier.
• t: The current simulation time.
• θ: The parameter vector.

Throws

• Always throws an error: "This function should not be called. Implement apply_all_non_callback_laws! for your own model."

apply_t_grad!(climate::RasterStack, dem::Raster)

Apply temperature gradients to the climate data based on the digital elevation model (DEM).

Arguments

• climate::RasterStack: A RasterStack object containing climate data, including temperature and gradient information.
• dem::Raster: A Raster object representing the digital elevation model (DEM) data.

Description

This function adjusts the temperature data in the climate object by applying the temperature gradients. The adjustment is based on the difference between the mean elevation from the DEM data and a reference height specified in the metadata of the climate object.

block_average(mat::Matrix{F}, n::Int) where {F <: AbstractFloat}

Downsamples a matrix by averaging non-overlapping n x n blocks. Returns a matrix of the block-averaged values with size (div(X, n), div(Y, n)) where (X, Y) = size(mat).

Arguments

• mat::Matrix{F}: Input 2D matrix.
• n::Int: Block size for downsampling. Both matrix dimensions must be divisible by n.

block_average_pad_edge(mat::Matrix{F}, n::Int) where {F <: AbstractFloat}

Downsamples a matrix by averaging n x n blocks, using edge-replication padding when the matrix dimensions are not divisible by n. Edge padding replicates the last row/column values to expand the matrix so that both dimensions are divisible by n. Returns a matrix of averaged values with size (ceil(Int, X/n), ceil(Int, Y/n)).

Arguments

• mat::Matrix{F}: Input 2D matrix.
• n::Int: Block size for downsampling.

block_average_pad_edge_masked(
    mat::Matrix{F},
    mask::BitMatrix,
    n::Int;
    empty_value::F = F(NaN),
) where {F <: AbstractFloat}

Downsamples a matrix by averaging n x n blocks intersecting with a mask, using edge-replication padding when the matrix dimensions are not divisible by n. Edge padding replicates the last row/column values to expand the matrix so that both dimensions are divisible by n. Returns a matrix of averaged values with size (ceil(Int, X/n), ceil(Int, Y/n)). The average discards values where mask is false. If the mask is full of falses for a given block, the average is replaced by the prescribed empty value.

Arguments

• mat::Matrix{F}: Input 2D matrix.
• mask::BitMatrix: Mask of valid data to average.
• n::Int: Block size for downsampling.
• empty_value::F: Fallback value for blocks that do not have valid values. Defaults to NaN.

build_affect(law::AbstractLaw, cache, glacier_idx, θ)

Return a !-style function suitable for use in a callback, which applies the given law to update the cache for a specific glacier and parameters θ, using the simulation time.

check_concrete_types(x, indent=0, show=true)

Check recursively that the fields of a nested struct are all concrete. Print information for each field and return false if any of the fields is not concrete.

check_field_types(x, indent=0; show=true)

Check recursively that the fields of type are all concrete. Print information for each field and return false if any of the fields is not concrete.

combine_velocity_data(refVelocities; merge=false)

Combine multiple ice surface velocity datasets into a single SurfaceVelocityData object.

Arguments

• refVelocities::Vector{SurfaceVelocityData}: A vector of ice surface velocity datasets to combine. Each element must have the same grid alignment (isGridGlacierAligned must be true for all).
• merge::Bool=false: If true, velocities with the same date are averaged, and corresponding date ranges (date1, date2) are reduced to their min/max. If false, data are simply concatenated.

Returns

• SurfaceVelocityData: A single object containing the combined velocity data, including vx, vy, vabs and their associated errors, as well as coordinate (x, y, lat, lon) and date information. The isGridGlacierAligned field reflects whether all input datasets were aligned.

Notes

• The function asserts that all input datasets are aligned on the same grid.
• When merge=true, velocities and errors are averaged over datasets sharing the same date.
• Uses nanmean for averaging to handle missing data.
• date_error is set to nothing when merging. # Check all surfaces are on the same grid as the glacier

compute_surface_topography(
    S::Matrix{<: AbstractFloat},
    Δx::Sleipnir.Float,
    Δy::Sleipnir.Float;
    window_m::Sleipnir.Float = Sleipnir.Float(200.0),
)

Compute dynamic slope and aspect fields from the current glacier surface S. The surface is spatially smoothed with a square moving window (window_m) before computing finite-difference gradients to reduce pixel-scale noise.

create_results(
    simulation::SIM,
    glacier_idx::I,
    solution,
    tstops::Vector{F};
    processVelocity::Union{Nothing, Function} = nothing,
) where {SIM <: Simulation, I <: Integer}

Create a Results object from a given simulation and solution.

Arguments

• simulation::SIM: The simulation object of type Simulation.
• glacier_idx::I: The index of the glacier within the simulation.
• solution: The solution object containing all the steps including intermediate ones.
• tstops::Vector{F}: The list of time steps to use to construct the results.
• processVelocity::Union{Nothing, Function}: Post processing function to map the ice thickness to the surface velocity. It is called before creating the results. It takes as inputs simulation, ice thickness (matrix) and the associated time and returns 3 variables Vx, Vy, V which are all matrix. Defaults is nothing which means no post processing is applied.

Returns

• results: A Results object containing the processed simulation data.

Details

The function processes the solution to select the last value for each time step. It then constructs a Results object containing various attributes from the simulation and the iceflow model.

downscale_2D_climate!(glacier::Glacier2D)

Update the 2D climate structure for a given glacier by downscaling climate data.

Arguments

• glacier::Glacier2D: The glacier object containing the climate data to be downscaled.

Description

This function updates the 2D climate structure of the given glacier by:

1. Updating the temperature, PDD (Positive Degree Days), snow, and rain fields in the 2D climate step with the corresponding values from the climate step.
2. Updating the gradients and average gradients in the 2D climate step.
3. Applying temperature gradients and computing the snow/rain fraction for the selected period by reprojecting the current S with the RasterStack structure.

Notes

# Update 2D climate structure

• The function modifies the glacier object in place.

downscale_2D_climate(climate_step::ClimateStep, S::Matrix{<: AbstractFloat}, Coords::Dict)

Downscales climate data to a 2D grid based on the provided matrix of surface elevation and coordinates.

Arguments

• climate_step::ClimateStep: A struct containing climate data for a specific time step. Expected fields are:

  • "avg_temp": Average temperature.
  • "temp": Temperature.
  • "prcp": Precipitation.
  • "gradient": Temperature gradient.
  • "avg_gradient": Average temperature gradient.
  • "ref_hgt": Reference height.

• S::Matrix{<: AbstractFloat}: Surface elevation data.

• Coords::Dict: A dictionary with keys "lon" and "lat" for longitude and latitude coordinates.

Returns

• Climate2Dstep: A Climate2Dstep object containing the downscaled climate data with fields:

  • temp: 2D array of temperature.
  • PDD: 2D array of positive degree days.
  • snow: 2D array of snow precipitation.
  • rain: 2D array of rain precipitation.
  • gradient: Temperature gradient.
  • avg_gradient: Average temperature gradient.
  • x: Longitude coordinates.
  • y: Latitude coordinates.
  • ref_hgt: Reference height.

Description # Create dummy 2D arrays to have a base to apply gradients afterwards

This function creates dummy 2D arrays based on the provided surface elevation data and applies the climate step data to these arrays. It then constructs a Climate2Dstep object with the downscaled climate data and applies temperature gradients to compute the snow/rain fraction for the selected period.

emptyPrepVJP(cache, vjpPrep, simulation, glacier_idx, t, θ)
emptyPrepVJPWithInputs(cache, vjpPrep, inputs, θ)

Function that defines an empty !-style function for the preparation of the VJPs. The two methods define the two possible signatures for the function that updates the cache in-place. Trying to apply this function will yield an error. It is for internal use only and it isn't exposed to the user.

emptyVJP(cache, simulation, glacier_idx, t, θ)
emptyVJPWithInputs(cache, inputs, θ)

Function that defines an empty !-style function for the VJPs. The two methods define the two possible signatures for the function that updates the cache in-place. It is for internal use only and it isn't exposed to the user.

fake_interpolated_datacube()

Create a fake datacube of ice surface velocity time series. It corresponds to the interpolated data.

fake_multi_datacube()

Create a fake datacube of ice surface velocity time series. It corresponds to the filtered multi source data.

fillNaN!(A::AbstractArray, fill::Number=zero(eltype(A)))

Replace all NaN values in the array A with the specified fill value.

Arguments

• A::AbstractArray: The array in which NaN values will be replaced.
• fill::Number: The value to replace NaN with. Defaults to zero(eltype(A)).

fillNaN(A::AbstractArray, fill::Number=zero(eltype(A)))

Replace all NaN values in the array A with the specified fill value. If no fill value is provided, it defaults to the zero value of the element type of A.

Arguments

• A::AbstractArray: The input array that may contain NaN values.
• fill::Number: The value to replace NaNs with. Defaults to zero(eltype(A)).

Returns

• An array of the same type and shape as A, with all NaN values replaced by fill.

fillZeros!(A::AbstractArray, fill::Number=NaN)

Replace all zero elements in the array A with the specified fill value.

Arguments

• A::AbstractArray: The array in which to replace zero elements.
• fill::Number: The value to replace zero elements with. Defaults to NaN.

fillZeros(A::AbstractArray, fill::Number=NaN) -> AbstractArray

Replace all zero elements in the array A with the specified fill value.

Arguments

• A::AbstractArray: The input array in which zero elements are to be replaced.
• fill::Number: The value to replace zero elements with. Defaults to NaN.

Returns

• AbstractArray: A new array with zero elements replaced by the fill value.

filter_missing_glaciers!(rgi_ids::Vector{String}, params::Parameters)

Filter out glaciers that cannot be processed from the given list of RGI IDs.

Arguments

• rgi_ids::Vector{String}: A vector of RGI IDs representing glaciers.
• params::Parameters: A Parameters object containing simulation parameters.

Description

This function filters out glaciers from the provided rgi_ids list based on two criteria:

1. Glaciers that are marked as level 2 in the RGI statistics CSV file.
2. Glaciers listed in the missing_glaciers.jld2 file located in the params.simulation.working_dir directory.

Notes

TODO: see if this is necessary, otherwise remove

• The RGI statistics CSV file is downloaded from a remote server.
• If the missing_glaciers.jld2 file is not available, a warning is logged and the function skips this filtering step. # Check which glaciers we can actually process # TODO: see if this is necessary, otherwise remove

generate_raw_climate_files(rgi_id::String, simparams::SimulationParameters)

Generate raw climate files for a given RGI (Randolph Glacier Inventory) ID and simulation parameters.

Arguments

• rgi_id::String: The RGI ID for which to generate raw climate files.
• simparams::SimulationParameters: The simulation parameters containing the time span and RGI paths.

Description

This function generates raw climate files for a specified RGI ID if they do not already exist. It retrieves raw climate data, ensures the desired period is covered, crops the data to the desired time period, and saves the raw climate data to disk.

Details

1. Constructs the path to the RGI directory using the provided rgi_id and simparams.

2. Checks if the raw climate file for the specified time span already exists.

3. If the file does not exist:

   • Retrieves the raw climate data.
   • Ensures the desired period is covered by the climate data.
   • Crops the climate data to the desired time period.
   • Saves the cropped climate data to disk. # Initialize RGI path to be accessible outside the try block
   • Triggers garbage collection to free up memory.

get_cumulative_climate!(climate, t::Sleipnir.Float, step::Sleipnir.Float, gradient_bounds=[-0.009, -0.003])
get_cumulative_climate!(climate, period::StepRange{Date, Day}, gradient_bounds=[-0.009, -0.003])

Calculate and update the cumulative climate data for a given period. The user can choose between providing a specific time t and a time step step, or a time period defined by period.

Keyword arguments

• climate::Climate: The climate object containing raw climate data.
• gradient_bounds::Vector{Float64}: Optional. The bounds within which to clamp the gradient values. Default is [-0.009, -0.003]. Optional parameters to specify the time period:
• t::Sleipnir.Float: Time at which the cumulative climate data should be computed.
• step::Sleipnir.Float: Time step used to compute the cumulative climate data. Together with t they define a time period. or
• period::StepRange{Date, Day}: The time period for which to compute the cumulative climate data.

Updates

• climate.climate_raw_step: The raw climate data for the given period.
• climate.avg_temps: The average temperature for the given period.
• climate.avg_gradients: The average gradient for the given period.
• climate.climate_step.prcp: The cumulative precipitation for the given period.
• climate.climate_step.temp: The cumulative temperature for the given period.
• climate.climate_step.gradient: The cumulative gradient for the given period.
• climate.climate_step.avg_temp: The average temperature for the given period.
• climate.climate_step.avg_gradient: The average gradient for the given period.
• climate.climate_step.ref_hgt: The reference height from the raw climate data.

get_cumulative_climate(
    climate::RasterStack,
    gradient_bounds::Vector{Float64}=[-0.009, -0.003],
)

Calculate cumulative climate statistics from the given climate data.

Keyword arguments

• climate::RasterStack: A RasterStack object containing temperature, precipitation, and gradient data.
• gradient_bounds::Vector{Float64}: A two-element vector specifying the lower and upper bounds for the gradient values. Defaults to [-0.009, -0.003].

Returns

• climate_sum::ClimateStep: A struct containing the following fields:

  • "temp": The sum of positive degree days (PDDs) from the temperature data.
  • "prcp": The sum of precipitation data.
  • "gradient": The sum of gradient data, clipped within the specified bounds.
  • "avg_temp": The average temperature.
  • "avg_gradient": The average gradient.
  • "ref_hgt": The reference height from the climate metadata.

Notes

• The temperature data is modified to only include positive degree-day values (PDDs).
• The gradient data is clipped within the specified bounds to ensure plausible values.

get_glathida!(glaciers::Vector{G}, params::Parameters; force=false) where {G <: Glacier2D}

Retrieve and process glacier thickness data for a vector of Glacier2D objects.

Arguments

• glaciers::Vector{Glacier2D}: A vector of Glacier2D objects for which the glacier thickness data is to be retrieved.
• params::Parameters: A Parameters object containing simulation parameters.
• force::Bool=false: A boolean flag indicating whether to force the retrieval of glacier thickness data.

Returns

• gtd_grids::Vector: A vector of glacier thickness data grids.
• glaciers::Vector{Glacier2D}: The updated vector of Glacier2D objects after removing glaciers with no data.

Description

This function retrieves glacier thickness data for each glacier in the input vector using parallel processing. It updates a list of missing glaciers if any glacier has all data points equal to zero. The function then removes glaciers with no data from both the gtd_grids and glaciers vectors and returns the updated vectors.

Notes

• The function uses pmap for parallel processing of glaciers.
• The list of missing glaciers is stored in a JLD2 file located at params.simulation.working_dir/data/missing_glaciers.jld2.
• Glaciers with no data are identified and removed based on the condition that all data points in their thickness grid are zero.

get_glathida_glacier(glacier::Glacier2D, params::Parameters, force)

Retrieve or generate the glathida glacier grid for a given glacier.

Arguments

• glacier::Glacier2D: The glacier object for which the glathida grid is to be retrieved or generated.
• params::Parameters: The parameters object containing simulation settings.
• force: A boolean flag indicating whether to force regeneration of the glathida grid even if it already exists.

Returns

• gtd_grid: A 2D array representing the glathida glacier grid.

Description

This function checks if the glathida glacier grid file (glathida.h5) exists in the specified path. If the file exists and force is false, it reads the grid from the file. Otherwise, it reads the glacier thickness data from a CSV file (glathida_data.csv), computes the average thickness for each grid cell, and saves the resulting grid to an HDF5 file (glathida.h5).

get_longterm_temps(rgi_id::String, params::Parameters, climate::RasterStack) -> Array{Float64}

Calculate the long-term average temperatures for a given glacier.

Arguments

• rgi_id::String: The RGI (Randolph Glacier Inventory) identifier for the glacier.
• params::Parameters: A struct containing simulation parameters, including paths to RGI data.
• climate::RasterStack: A RasterStack object containing climate data.

Returns

• Array{Float64}: An array of long-term average temperatures.

Description

This function retrieves the gridded data for the specified glacier using its RGI identifier. It then applies a temperature gradient to the climate data based on the glacier's topography. Finally, it calculates the long-term average temperatures by grouping the temperature data by year and computing the mean for each group.

get_raw_climate_data(rgi_path::String) -> RasterStack

Load raw climate data from a specified path.

Arguments

• rgi_path::String: The file path to the directory containing the climate data file.

Returns

• RasterStack: A RasterStack object containing the climate data from the specified file.

get_result_id_from_rgi(glacier_id::I, simulation::SIM) where {I <: Integer, SIM <: Simulation}

Extract results of specific simulation from the Simulation object.

Arguments

• glacier_id::I: Numerical ID of glacier used to generate simulation.
• simulation::SIM`: The simulation object containing the parameters and results.

get_winter_prcp_factor(glacier::AbstractGlacier, params::Parameters; prcp_fac_bounds) -> Float

Compute a glacier-specific precipitation correction factor from mean winter precipitation, following OGGM's decide_winter_precip_factor.

Arguments

• glacier::AbstractGlacier: Glacier providing rgi_id, cenlat (hemisphere) and climate.climate_data_source (:W5E5 or :ERA5).
• params::Parameters: Simulation parameters, used to locate the climate file.
• prcp_fac_bounds: Lower/upper clip bounds (default (0.1, 10.0), matching OGGM).

Description

Reads the full daily climate record (1979–2019, the window OGGM's coefficients were fit on), averages winter daily precipitation (NH: Oct–Apr, SH: Apr–Oct) in kg/m²/day, and maps it to a precipitation factor: log fit for W5E5, linear fit for ERA5. If the daily file is missing, falls back to 2.5. Note: we average daily values directly (OGGM averages monthly daily-rates; equivalent within noise).

glacierName(rgi_id::String)
glacierName(rgi_ids::Vector{String})

Returns the name(s) of one or multiple glaciers based the given RGI ID(s). It uses the rgi62_stats.csv file from OGGM.

grid(
    glacier::G,
    latitudes::Vector{F},
    longitudes::Vector{F},
    vx::Union{FileArray, Array{Union{Missing, F}, 3}},
    vy::Union{FileArray, Array{Union{Missing, F}, 3}},
    mapping::VM
) where {
    G <: AbstractGlacier,
    F <: AbstractFloat,
    VM <: VelocityMapping,
    FileArray <: Rasters.FileArray
}

Grid velocity data onto the glacier grid following the prescribed mapping. This function maps the 3 dimensional surface velocities (x, y and t) to the glacier grid. The provided surface velocities can be a Rasters.FileArray which happens when the RasterStack is instantiated in lazy mode. In this situation, only the smallest cube that contains all the needed data to construct the mapping is read from disk. The returned velocity variables have shape (nTimes, nx, ny) where nTimes is the number of time steps and (nx, ny) is the size of the glacier grid.

Arguments:

• glacier::G: Glacier instance which determines the glacier on which the velocities are projected onto.
• latitudes::Vector{F}: Vector of latitude values of the original surface velocity grid.
• longitudes::Vector{F}: Vector of longitude values of the original surface velocity grid.
• vx::Union{FileArray, Array{Union{Missing, F}, 3}}: X component of the original surface velocities. It can be either a Rasters.FileArray if the datacube is read in lazy mode, or a plain 3 dimensional array.
• vy::Union{FileArray, Array{Union{Missing, F}, 3}}: Y component of the original surface velocities. It can be either a Rasters.FileArray if the datacube is read in lazy mode, or a plain 3 dimensional array.
• mapping::VM: Mapping to use.

Returns:

• xG: A vector that gives the x coordinates of the glacier grid.
• yG: A vector that gives the y coordinates of the glacier grid.
• vxG: A 3 dimensional array of the x component of the velocity gridded onto the glacier grid.
• vyG: A 3 dimensional array of the y component of the velocity gridded onto the glacier grid.

indFromT(tspan, tstops, t)

Find indices in time vector t corresponding to target times in tstops.

Arguments

• tspan: Tuple of simulation start and end times
• tstops: Vector of target times for which to find corresponding indices
• t: Solution time vector containing all time steps

Returns

Vector of indices such that t[indices] are the times in t closest to each value in tstops.

Details

• For the initial time step (tspan[1]), returns the first matching index to avoid out-of-range interpolation
• For subsequent time steps, returns the last matching index to capture state after potential callbacks
• Uses approximate equality with relative tolerance rtol=1e-10 to handle numerical rounding

initialize_glacier(rgi_id::String, parameters::Parameters; smoothing=false)

Initialize a glacier with the given RGI ID and parameters.

Arguments

• rgi_id::String: The RGI (Randolph Glacier Inventory) ID of the glacier.
• parameters::Parameters: A struct containing various parameters required for initialization.
• smoothing::Bool: Optional. If true, apply smoothing to the initial topography. Default is false.
• masking::Union{Int, Nothing, Matrix}: Type of mask applied to the glacier to determine regions with no ice.
• velocityDatacubes::Union{Dict{String, String}, Dict{String, RasterStack}}: A dictionary that provides for each RGI ID either the path to the datacube or the RasterStack with velocity data.

Returns

• glacier: An initialized glacier object containing the initial topography and climate data.

initialize_glaciers(
    rgi_ids::Vector{String},
    params::Parameters;
    velocityDatacubes::Union{Dict{String, String}, Dict{String, RasterStack}}=Dict(),
)

Initialize glaciers based on provided RGI IDs and parameters.

Arguments

• rgi_ids::Vector{String}: A vector of RGI IDs representing the glaciers to be initialized.
• params::Parameters: A Parameters object containing simulation parameters.
• velocityDatacubes::Union{Dict{String, String}, Dict{String, RasterStack}}: A dictionary that provides for each RGI ID either the path to the datacube or the RasterStack with velocity data.

Returns

• glaciers::Vector{Glacier2D}: A vector of initialized Glacier2D objects.

Description

This function performs the following steps:

1. Generates a file for missing glaciers if it does not already exist.
2. Filters out missing glaciers from the provided RGI IDs.
3. Generates raw climate data for the glaciers if necessary.
4. Initializes the glaciers using the provided RGI IDs and parameters.
5. If use_glathida_data is enabled in the simulation parameters, assigns GlaThiDa data to the glaciers.

Errors

• Throws an error if none of the provided RGI IDs have GlaThiDa data.

Warnings

• Issues a warning if not all glaciers have GlaThiDa data available.

Example

# We declare a list of glaciers to be initialized with their RGI IDs
rgi_ids = ["RGI60-11.03638", "RGI60-11.01450", "RGI60-11.02346", "RGI60-08.00203"]
# We initialize those glaciers based on the RGI IDs and the parameters we previously specified
glaciers = initialize_glaciers(rgi_ids, params)

initialize_surfacevelocitydata(
    raster::Union{String, <: RasterStack};
    glacier::Union{G, Nothing}=nothing,
    mapping::VM=MeanDateVelocityMapping(),
    compute_vabs_error::Bool=true,
    flag::Union{String, <: RasterStack, Nothing} = nothing
) where {G <: AbstractGlacier, VM <: VelocityMapping}

Initialize SurfaceVelocityData from Rabatel et. al (2023).

Arguments:

• raster::Union{String, RasterStack}: RasterStack or path of the netCDF file with surface velocity data.
• glacier::Union{G, Nothing}: Glacier associated to the surface velocity datacube. When provided, the surface velocity data are gridded on the glacier grid using the mapping.
• mapping::VM: Mapping to use in order to grid the data from the coordinates of the velocity product datacube to the glacier grid.
• compute_vabs_error::Bool: Whether to compute the absolute error uncertainty.
• flag::Union{String, <: RasterStack, Nothing}: Option to provide a RasterStack containing a fmask raster and which provides an indicator of whether each pixel of the ice surface velocity data is considered as reliable or not.

initialize_surfacevelocitydata_mask(data::RasterStack, flag::Union{String, <:RasterStack, Nothing}=nothing)

Create a mask for a glacier surface velocity dataset by subsetting and aligning a flag raster to the spatial domain of data.

Arguments

• data::RasterStack: The glacier datacube or surface velocity dataset.
• flag::Union{String, RasterStack, Nothing}: Type of flag to be applied

Description

This function extracts the portion of the flag raster that spatially overlaps with the glacier domain defined by data. Because the flag raster is centered on pixel centers, its bounding box is shifted by half a grid step when subsetting.

is_in_glacier(A::Matrix{F}, distance::I) where {I <: Integer, F <: AbstractFloat}

Return a matrix with booleans indicating if a given pixel is at distance at least distance in the set of non zero values of the matrix. This usually allows discarding the border pixels of a glacier. A positive value of distance indicates a measurement from inside the glacier, while a negative distance indicates one from outside.

Arguments:

• A::Matrix{F}: Matrix from which to compute the matrix of booleans.
• distance::I: Distance to the border, computed as the number of pixels we need to move from within the glacier to find a pixel with value zero.

Compute distance between one point in the format of a TransverseMercator point and a set of points defined through the coordinates x and y in meters.

Generate a thickness-evolution animation. The output format (e.g. MP4 or GIF) is inferred from the extension of pathVideo, since Makie's record selects the encoder from the file extension.

mapVelocity(
    velocityMapping::MeanDateVelocityMapping,
    velocityData::SurfaceVelocityData,
    t::AbstractFloat,
)

Retrieve the reference ice surface velocity for a given time step. This mapping uses the nearest snap shot available within a time window whose length is controlled by velocityMapping.thresDate. If no snapshot is found in the time window of length 2*velocityMapping.thresDate, the returned ice velocity components are empty matrices and the returned boolean flag useVel is set to false.

Arguments:

• velocityMapping::MeanDateVelocityMapping: Mapping to map the reference ice velocity to a target time step t.
• velocityData::SurfaceVelocityData: Surface velocity data. This is usually an attribute of a glacier.
• t::AbstractFloat: Current time step.

Returns

• Vx_ref: Matrix of the x-component of the ice surface velocity.
• Vy_ref: Matrix of the y-component of the ice surface velocity.
• V_ref: Matrix of the ice surface velocity magnitude.
• useVel: Boolean indicating whether the returned ice surface velocity can be used or not. The value of this boolean depends on the success of the ice surface velocity mapping at the current time step t.

max_or_empty(A::Array)

Return maximum value for non-empty arrays. This is just required to compute the error in the absolute velocity.

parse_proj(proj::String)

Parses the string containing the information of the projection to filter for important information "+proj=tmerc +lat0=0 +lon0=6.985 +k=0.9996 +x0=0 +y0=0 +datum=WGS84 +units=m +no_defs"

partial_year(float::Sleipnir.Float) -> Sleipnir.Float

Calculate the partial year value based on the given floating-point number.

Arguments

• float::Sleipnir.Float: A floating-point number representing the fraction of the year.

Returns

• Sleipnir.Float: The calculated partial year value.

partial_year(period::Type{<:Period}, float)

Calculate a partial year date based on a floating-point year value.

Arguments

• period::Type{<:Period}: The type of period to use (e.g., Month, Day).
• float::Sleipnir.Float: The floating-point year value.

Returns

• Date: The calculated date corresponding to the partial year.

plot_bias(
    results,
    variables;
    treshold = [0, 0],
    figsize::Union{Nothing, Tuple{Int64, Int64}}=nothing,
)

Plot the bias of the glacier integrated volume over the specified time span.

Arguments

• results::Results: The results object containing the data to be plotted.
• variables::Vector{Symbol}: The variables to be plotted.
• title_mapping::Dict{Symbol, String}: A dictionary mapping variable names to their titles.
• tspan::Tuple{Float64, Float64}: A tuple representing the start and end time for the simulation.
• figsize::Union{Nothing, Tuple{Int64, Int64}}: Size of the figure.

Returns

• A plot of the glacier integrated volume bias.

plot_cumulative_gridded_data(
    gridded_data::Vector{Matrix{F}},
    results::Results;
    weights::Union{Nothing,AbstractVector{<:Real}}=nothing,
    kwargs...
) where {F <: AbstractFloat}

Plot the cumulative field of a time series of gridded matrices using plot_gridded_data. This is a thin utility wrapper to avoid duplicating plotting code.

Arguments

• gridded_data::Vector{Matrix{F}}: Sequence of gridded fields to accumulate.
• results::Results: Results object with glacier metadata for plotting.
• weights::Union{Nothing,AbstractVector{<:Real}}: Optional per-step weights.
• kwargs...: Additional keyword arguments forwarded to plot_gridded_data.

Returns

• Figure: Cumulative field figure.

plot_cumulative_mb(results::Results; kwargs...)

Plot a spatial map of the cumulative surface mass balance accumulated over the simulation period from the per-callback MB fields stored in results.MB.

Each entry of results.MB is the gridded mass-balance increment produced at one MB callback (one per MB time step). The function sums these increments cell-by-cell (via accumulate_gridded_data) to obtain the total mass balance over the period results.tspan = (t0, t1).

Modes (annual_MB)

• annual_MB = false (default): plots the raw cumulative MB over the whole period, i.e. Σᵢ MBᵢ, in m w.e..
• annual_MB = true: divides every increment by the period length T = t1 - t0 (in years) before summing, i.e. Σᵢ (MBᵢ / T) = (Σᵢ MBᵢ) / T. Since Σᵢ MBᵢ is the total over T years, this is the mean annual mass-balance rate (the "annually-averaged equivalent"), in m w.e. yr⁻¹.

Returns nothing (with a warning) when results.MB is empty — e.g. when the simulation was run without use_MB = true.

Arguments

• results::Results: results carrying the per-callback MB maps in results.MB.

Keyword arguments

• title::String: figure title prefix (the period (t0–t1) is appended automatically).
• colorbar_label::Union{Nothing,String}: colorbar label; defaults to m w.e. (or m w.e. yr⁻¹ when annual_MB = true).
• annual_MB::Bool = false: switch between cumulative total and mean annual rate.
• colormap: diverging colormap (red→white→blue by default).
• kwargs...: forwarded to plot_gridded_data.

Returns

• Figure, or nothing if there is no MB history.

plot_glacier(results, plot_type, variables; kwargs...) -> Figure

High-level entry point that dispatches to specific glacier plotting functions based on plot_type:

• "heatmaps" → plot_glacier_heatmaps
• "quivers" → plot_glacier_quivers
• "evolution difference" → plot_glacier_difference_evolution
• "evolution statistics" → plot_glacier_statistics_evolution
• "integrated volume" → plot_glacier_integrated_volume
• "bias" → plot_bias
• "dem" → plot_glacier_dem

plot_glacier_dem(glacier_or_results; kwargs...)

Plot the glacier DEM (surface elevation field S) with a terrain colormap, geographic coordinate labels, colorbar, and glacier contour overlay.

Arguments

- `glacier_or_results`: Either `Results` or `Glacier2D`.

Keyword Arguments

- `title::String`: Figure title prefix.
- `colorbar_label::String`: Label for the colorbar.
- `plotContour::Bool`: Whether to overlay glacier contour (default `true`).
- `colormap`: Makie colormap symbol (default `:terrain`).
- `kwargs...`: Additional keyword arguments forwarded to `plot_gridded_data`.

Returns

- `Figure`: DEM figure.

plot_glacier_difference_evolution(
    results::Results,
    variables::Vector{Symbol},
    title_mapping;
    tspan::Tuple{F,F}=results.tspan,
    metrics::Vector{String}="difference",
    figsize::Union{Nothing, Tuple{Int64, Int64}}=nothing,
) where {F<:AbstractFloat}

Plot the evolution of the difference in a glacier variable over time.

Arguments

• results::Results: The simulation results object containing the data to be plotted.
• variables::Vector{Symbol}: The variable to be plotted.
• title_mapping: A dictionary mapping variable names to their titles.
• tspan::Tuple{F,F}: A tuple representing the start and end time for the simulation.
• metrics::Vector{String}: Metrics to visualize, e.g., ["difference"].
• figsize::Union{Nothing, Tuple{Int64, Int64}}: Size of the figure.

Returns

• A plot of the glacier difference evolution.

plot_glacier_heatmaps(
    results::Results,
    variables::Vector{Symbol},
    title_mapping::Dict;
    scale_text_size::Union{Nothing,Float64}=nothing,
    timeIdx::Union{Nothing,Int64}=nothing,
    figsize::Union{Nothing, Tuple{Int64, Int64}} = nothing,
    plotContour::Bool=false,
) -> Figure

Plot heatmaps for glacier variables.

Arguments

• results::Results: The results object containing the data to be plotted.
• variables::Vector{Symbol}: A list of variables to be plotted.
• title_mapping::Dict: A dictionary mapping variable names to their titles and colormaps.
• scale_text_size::Union{Nothing,Float64}: Optional argument to scale the text size.
• timeIdx::Union{Nothing,Int64}:: Optional argument to select the index at which data should be plotted when dealing with vector of matrix. Default is nothing which selects the last element available.
• figsize::Union{Nothing, Tuple{Int64, Int64}}: Size of the figure.
• plotContour::Bool: Whether to add a contour plot representing the glacier borders at the beginning of the simulation on top of each of the figures. Default is false.

Returns

• A plot of the glacier heatmaps.

plot_glacier_integrated_volume(
    results,
    variables,
    title_mapping;
    tspan,
    figsize::Union{Nothing, Tuple{Int64, Int64}}=nothing,
)

Plot the integrated volume of a glacier variable over time.

Arguments

• results::Results: The results object containing the data to be plotted.
• variables::Vector{Symbol}: The variable to be plotted.
• title_mapping: A dictionary mapping variable names to their titles.
• tspan: A tuple representing the start and end time for the simulation.
• figsize::Union{Nothing, Tuple{Int64, Int64}}: Size of the figure.

Returns

• A plot of the glacier integrated volume.

plot_glacier_quivers(
    results::Results,
    variables::Vector{Symbol},
    title_mapping::Dict;
    timeIdx::Union{Nothing,Int64} = nothing,
    figsize::Union{Nothing, Tuple{Int64, Int64}} = nothing,
    lengthscale::Float64 = Float64(0.00001),
    tiplength::Float64 = Float64(0.5),
) -> Figure

Plot quivers for glacier variables.

Arguments

• results::Results: The results object containing the data to be plotted.
• variables::Vector{Symbol}: A list of variables to be plotted.
• title_mapping::Dict: A dictionary mapping variable names to their titles and colormaps.
• timeIdx::Union{Nothing,Int64}:: Optional argument to select the index at which data should be plotted when dealing with vector of matrix. Default is nothing which selects the last element available.
• figsize::Union{Nothing, Tuple{Int64, Int64}}: Size of the figure.
• lengthscale::Float64: Lengthscale of the arrows in the quiver plot.
• tiplength::Float64: Length of the arrow in the quiver plot.

Returns

• A plot of the glacier quivers.

plot_glacier_statistics_evolution(
    results::Results,
    variables::Vector{Symbol},
    title_mapping;
    metrics="median",
    tspan,
    threshold=0.5,
    figsize::Union{Nothing, Tuple{Int64, Int64}}=nothing,
)

Plot the evolution of statistics for multiple glacier variables over time.

Arguments

• results::Results: The simulation results object containing the data to be plotted.
• variables::Vector{Symbol}: A list of variables to be plotted.
• title_mapping: A dictionary mapping variable names to their titles.
• metrics: Metrics to visualize, e.g., "average", "median", "min", "max", and "std". Default is "median".
• tspan: A tuple representing the start and end time for the simulation.
• threshold: A threshold value to filter the data. Default is 0.5.
• figsize::Union{Nothing, Tuple{Int64, Int64}}: Size of the figure.

Returns

• A plot of the glacier statistics evolution.

plot_glacier_vid(
    plot_type::String,
    results::Results,
    glacier::Glacier2D,
    tspan,
    step,
    pathVideo::String;
    framerate::Int=24,
    baseTitle::String=""
)

Generate various types of videos for glacier data. For now only the evolution of the glacier ice thickness is supported. More types of visualizations will be added in the future.

Arguments

• plot_type: Type of plot to generate. Options are:

  • "thickness": Heatmap of the glacier thickness.

• results: A result object containing the simulation results including ice thickness over time.

• glacier: A glacier instance.

• tspan: The simulation time span.

• step: Time step to use to retrieve the results and generate the video.

• pathVideo: Path of the output animation. The format is inferred from the file extension — e.g. .mp4 for a video or .gif for an animated GIF.

Optional Keyword Arguments

• framerate: The framerate to use for the video generation.
• baseTitle: The prefix to use in the title of the frames. In each frame it is concatenated with the value of the year in the form " (t=XXXX)".

plot_gridded_data(
    gridded_data::Union{Vector{Matrix{F}}, Matrix{F}},
    results::Results;
    scale_text_size::Union{Nothing,Float64}=nothing,
    timeIdx::Union{Nothing,Int64}=nothing,
    figsize::Union{Nothing, Tuple{Int64, Int64}} = nothing,
    plotContour::Bool=false,
    colormap = :cool,
    logPlot = false,
) where {F <: AbstractFloat}

Plot a gridded matrix (or a time series of matrices) as a heatmap using metadata from results.

Arguments

• gridded_data::Union{Vector{Matrix{F}}, Matrix{F}}: Single snapshot or time series (defaults to last timestep).
• results::Results: Supplies lon, lat, x, y, rgi_id, Δx and H (mask).
• scale_text_size, figsize, colormap: Optional plotting params.
• timeIdx::Union{Nothing,Int64}: Select timestep when gridded_data is a vector.
• plotContour::Bool: overlay glacier-mask contour from results.H.
• logPlot::Bool: Use log10 colorscale (positive non-NaN values determine range).

Behavior

• Masks out cells where results.H[begin] .<= 0 (set to NaN).
• Adds colorbar, central lon/lat tick, and a Δx-based scale bar in km.
• If plotContour, draws mask boundary lines.
• Returns a CairoMakie.Figure.

Errors

• Asserts gridded_data is non-empty and timeIdx (if provided) is in range.

precompute_all_VJPs_laws!(model::AbstractModel, cache, simulation, glacier_idx, t, θ)

This function is a placeholder and must be implemented for your custom model type.

It is intended to precompute the VJPs for all the laws that are used in a model. By default, calling this function will throw an error to indicate that the user should provide their own implementation tailored to their model.

Arguments

• model::AbstractModel: The model instance.
• cache: The cache object storing state variables.
• simulation: The simulation context.
• glacier_idx: Index identifying the glacier.
• t: The current simulation time.
• θ: The parameter vector.

Throws

• Always throws an error: "This function should not be called. Implement precompute_all_VJPs_laws! for your own model."

prepare_vjp_law(simulation, law::AbstractLaw, law_cache, θ, glacier_idx)

Function used to prepare the VJPs at the initialization of the model cache. It is used for example to compile VJPs of the laws to be differentiated using DifferentiationInterface.jl.

random_spatially_coherent_mask(h::Integer, w::Integer; sigma::Real=1.0, threshold::Real=0.0) -> BitMatrix
random_spatially_coherent_mask(mask::BitMatrix; sigma::Real=1.0, threshold::Real=0.0) -> BitMatrix

Generate a random binary mask with spatially correlated patches rather than pixel-wise independent noise. This is done by drawing white noise, applying a Gaussian low-pass filter in the frequency domain, and thresholding the result.

Arguments

• h::Integer, w::Integer: Height and width of the mask.
• mask::BitMatrix: An existing binary mask. The generated spatially coherent mask will be applied elementwise (.&) to this mask.
• sigma::Real=1.0: Controls the spatial correlation length. Larger values produce smoother, larger patches.
• threshold::Real=0.0: Threshold applied to the filtered noise. Higher values result in sparser masks. Statistically, setting the threshold to zero results in a mask with half pixels to true.

Returns

A BitMatrix of size (h, w) containing true in patchy regions and false elsewhere.

Examples

# Generate a new 256×256 patchy mask
mask = random_spatially_coherent_mask(256, 256; sigma = 8.0, threshold = 0.0)

# Apply patchy masking to an existing mask
base = trues(128, 128)
patchy = random_spatially_coherent_mask(base; sigma = 5.0, threshold = 0.3)    # 1) white noise

ratio_max(v, vabs)

Compute the maximum ratio between v and vabs at points where the value of vabs is not a NaN.

retrieve_simulation(p)
retrieve_simulation(p::Container)

Function that retrieves the simulation object from integrator.p when called from a callback. If p is a subtype of Container, then p.simulation is returned, otherwise it returns p. It is for internal use only and it isn't exposed to the user.

reverseForHeatmap(inp, x, y) -> Matrix

Out-of-place reverse of a matrix so that the heatmap is displayed in the correct geographic orientation. Flips along each axis whose coordinate vector is descending.

Arguments

• inp::Matrix{F}: The matrix to reverse.
• x::Vector{F}: Values of the x axis.
• y::Vector{F}: Values of the y axis.

Returns

• Out-of-place copy of inp that has been reversed if needed.

save_figure(fig, path)

Save fig to path, creating any missing parent directories automatically. Returns path.

save_results_file!(results_list::Vector{Results{F, I}}, simulation::SIM; path::Union{String,Nothing}=nothing) where {F <: AbstractFloat, I <: Integer, SIM <: Simulation}

Save the results of a simulation to a file.

Arguments

• results_list::Vector{Results{F, I}}: A vector containing the results of the simulation.
• simulation::SIM: The simulation object containing the parameters and results.
• path::Union{String,Nothing}: Optional. The path where the results file will be saved. If not provided, a default path will be used.

Description

This function saves the results of a simulation to a file in JLD2 format. If the path argument is not provided, the function will create a default path based on the current project directory. The results are saved in a file named prediction_<nglaciers>glaciers_<tspan>.jld2, where <nglaciers> is the number of glaciers in the simulation and <tspan> is the simulation time span.

smooth!(A)

Smooths the interior of a 2D array A using a simple averaging method. The function modifies the array A in place.

Arguments

• A::AbstractMatrix: A 2D array to be smoothed.

Details

The function updates the interior elements of A (excluding the boundary elements) by adding a weighted average of the second differences along both dimensions. The boundary elements are then set to the values of their nearest interior neighbors to maintain the boundary conditions.

stop_condition_tstops(u, t, integrator, tstops)

Check if the current time t is in the list of stop times tstops.

Arguments

• u: The current state of the system (not used in this function).
• t::AbstractFloat: The current time.
• integrator: The integrator object (not used in this function).
• tstops::Vector{<: AbstractFloat}: A collection of times at which the integration should stop.

Returns

• Bool: true if t is in tstops, otherwise false.

tdata(data::Nothing)
tdata(data::ThicknessData)
tdata(data::DhdtData)
tdata(data::Nothing, mapping::MeanDateVelocityMapping)
tdata(data::SurfaceVelocityData, mapping::MeanDateVelocityMapping)

Retrieve the time steps at which data is available for ice thickness and surface velocity data. If the provided data is nothing, returns an empty vector.

trim_period(period, climate)

Adjusts the given period to fit within the bounds of the climate data, ensuring it aligns with hydrological years.

Arguments

• period::UnitRange{Date}: The initial date range to be trimmed.
• climate::AbstractArray: The climate data array, which should have a time dimension Ti.

Returns

• UnitRange{Date}: The adjusted date range that fits within the climate data's time bounds.

Details

• If the start of the climate data is later than the start of the period, the period is adjusted to start from October 1st of the year of the climate data's start.
• If the end of the climate data is earlier than the end of the period, the period is adjusted to end on September 30th of the year of the climate data's end.

∂law∂inp!()

This function serves as a placeholder and should be replaced by other implementations in ODINN. This implementation throws an error. It is for internal use only and is not exposed to the user.

∂law∂θ!()

This function serves as a placeholder and should be replaced by other implementations in ODINN. This implementation throws an error. It is for internal use only and is not exposed to the user.