Code for preprocessing Rabatel (2023) data

Project.toml

@@ -9,6 +9,7 @@ CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 CodecZlib = "944b1d66-785c-5afd-91f1-9de20f533193"
 ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
 CoordRefSystems = "b46f11dc-f210-4604-bfba-323c1ec968cb"
+DateFormats = "44557152-fe0a-4de1-8405-416d90313ce6"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
@@ -17,6 +18,7 @@ Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
+NetCDF = "30363a11-5582-574a-97bb-aa9a979735b9"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Rasters = "a3a2b9e3-a471-40c9-b274-f788e487c689"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"

src/Sleipnir.jl

@@ -23,6 +23,8 @@ using Tar
 import NCDatasets
 using Unitful: m, rad, °
 using CoordRefSystems
+using Dates, DateFormats
+using NetCDF
 
 # ##############################################
 # ############    PARAMETERS     ###############
@@ -39,11 +41,13 @@ const global prepro_dir::String = joinpath(homedir(), ".ODINN", "ODINN_prepro")
 include("setup/config.jl")
 # All parameters needed for the models
 include("parameters/Parameters.jl")
-# Anything related to managing glacier topographical and climate data
+# Anything related to managing glacier topographical and climate variables
 include("glaciers/glacier/Glacier.jl")
+# Anything related to managing glacier data used for data assimilation
+include("glaciers/data/Data.jl")
 # All structures and functions related to ODINN models
 include("models/Model.jl")
 # Everything related to running simulations in ODINN
 include("simulations/Simulation.jl")
 
-end # module
+end # module

src/glaciers/data/Data.jl

@@ -0,0 +1,5 @@
+export AbstractData
+
+abstract type AbstractData end
+
+include("SurfaceVelocityData.jl")

src/glaciers/data/SurfaceVelocityData.jl

@@ -0,0 +1,75 @@
+export SurfaceVelocityData
+
+mutable struct SurfaceVelocityData{F <: AbstractFloat} <: AbstractData
+    x::Vector{F}
+    y::Vector{F}
+    lat::Vector{F}
+    lon::Vector{F}
+    vx::Array{F, 3}
+    vy::Array{F, 3}
+    vabs::Array{F, 3}
+    vx_error::Array{F, 1}
+    vy_error::Array{F, 1}
+    vabs_error::Array{F, 1}
+    date::Vector{DateTime}
+    date1::Vector{DateTime}
+    date2::Vector{DateTime}
+    date_error::Vector{Day}
+end
+
+"""
+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{Array{F, 3}, Nothing} = nothing,
+    vy::Union{Array{F, 3}, Nothing} = nothing,
+    vabs::Union{Array{F, 3}, Nothing} = nothing,
+    vx_error::Union{Array{F, 1}, Nothing} = nothing,
+    vy_error::Union{Array{F, 1}, Nothing} = nothing,
+    vabs_error::Union{Array{F, 1}, Nothing} = nothing,
+    date::Union{Vector{DateTime}, Nothing} = nothing,
+    date1::Union{Vector{DateTime}, Nothing} = nothing,
+    date2::Union{Vector{DateTime}, Nothing} = nothing,
+    date_error::Union{Vector{Day}, Nothing} = nothing,
+    ) where {F <: AbstractFloat}
+
+Constructor for ice surface velocity data based on Rabatel et. al (2023). 
+
+
+Important remarks:
+- Projections in longitude and latitude assume we are working in the north hemisphere. 
+  If working with south hemisphere glaciers, this needs to be changed.
+- 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).
+"""
+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{Array{F, 3}, Nothing} = nothing,
+    vy::Union{Array{F, 3}, Nothing} = nothing,
+    vabs::Union{Array{F, 3}, Nothing} = nothing,
+    vx_error::Union{Array{F, 1}, Nothing} = nothing,
+    vy_error::Union{Array{F, 1}, Nothing} = nothing,
+    vabs_error::Union{Array{F, 1}, Nothing} = nothing,
+    date::Union{Vector{DateTime}, Nothing} = nothing,
+    date1::Union{Vector{DateTime}, Nothing} = nothing,
+    date2::Union{Vector{DateTime}, Nothing} = nothing,
+    date_error::Union{Vector{Day}, Nothing} = nothing,
+    ) where {F <: AbstractFloat}
+
+    ft = typeof(x[begin])
+
+    return SurfaceVelocityData{ft}(x, y, lat, lon, vx, vy, vabs, vx_error, vy_error, vabs_error, date, date1, date2, date_error)
+end
+
+
+include("surfacevelocitydata_utils.jl")

src/glaciers/data/surfacevelocitydata_utils.jl

@@ -0,0 +1,99 @@
+export initialize_surfacevelocitydata
+export plot_timeseries, plot_count
+
+"""
+    initialize_surfacevelocitydata(file::String; interp=false)
+
+Initialize SurfaceVelocityData from Rabatel et. al (2023). 
+
+Arguments
+=================
+    - `file`: name of netcdf file with data
+    - `interp`: boolean variable indicating if we use the inporpolated data or not
+"""
+function initialize_surfacevelocitydata(file::String; interp=false)
+
+    # Date of first adquisition
+    date1 = ncread(file, "date1")
+    # Date of second adquisition
+    date2 = ncread(file, "date2")
+
+    # Middle date
+    if !interp 
+        date_mean = 0.5 .* date1 .+ 0.5 .* date2
+        date_mean_offset = date1_offset = date2_offset = 0
+    else 
+        date_mean = ncread(file, "mid_date")
+        date_mean_offset = datetime2julian(DateTime("2015-08-14")) - 2400000.5
+        date1_offset = datetime2julian(DateTime("2015-07-30")) - 2400000.5
+        date2_offset = datetime2julian(DateTime("2015-08-29")) - 2400000.5
+    end
+
+    # Convert dates from Modified Julian Days to datetipes
+    date1 = mjd.(date1 .+ date1_offset)
+    date2 = mjd.(date2 .+ date2_offset)
+    date_mean = mjd.(date_mean .+ date_mean_offset)
+    date_error = date2 .- date1
+
+    # We convert from Date to DateTime
+    date_mean = DateTime.(date_mean)
+    date1 = DateTime.(date1)
+    date2 = DateTime.(date2)
+
+    # Read data from netcdf file
+    x = ncread(file, "x")
+    y = ncread(file, "y")
+
+    # Velocity in the x direction (m/yr)
+    vx = ncread(file, "vx")              
+    vy = ncread(file, "vy")    
+
+    # Compute absolute velocity
+    vabs = (vx.^2 .+ vy.^2).^0.5  
+    # The sqrt operation in Julia promotes Float32 to Float64. We convert manually 
+    # to keep type consistency     
+    vabs = convert(typeof(vx), vabs)   
+
+    # Error is reported once per timespan, so upper bounds are given by absolute error
+    if !interp
+        vx_error = ncread(file, "error_vx")
+        vy_error = ncread(file, "error_vy")
+        # Absolute error uncertanty using propagation of uncertanties 
+        vx_ratio_max = map(i -> max_or_empyt(abs.(vx[:,:,i][vabs[:,:,i] .> 0.0]) ./ vabs[:,:,i][vabs[:,:,i] .> 0.0]), 1:size(vx)[3])
+        vy_ratio_max = map(i -> max_or_empyt(abs.(vy[:,:,i][vabs[:,:,i] .> 0.0]) ./ vabs[:,:,i][vabs[:,:,i] .> 0.0]), 1:size(vy)[3])
+        vabs_error = ((vx_ratio_max .* vx_error).^2 .+ (vy_ratio_max .* vy_error).^2).^0.5
+        vabs_error = convert(typeof(vx_error), vabs_error)  
+    else         
+        vx_error = nothing
+        vy_error = nothing
+        vabs_error = nothing
+    end
+
+    # Run some basic tests 
+    nx, ny, ntimes = size(vx)
+    @assert length(date1) == length(date2) == ntimes
+    @assert nx == ny == 250
+
+    # Spatial preprocessing
+    proj_zone = ncgetatt(file, "mapping", "utm_zone_number")
+    transform(X,Y) = Sleipnir.UTMercator(X, Y; zone=proj_zone, hemisphere=:north) 
+    Coordinates = transform.(x, y)
+    latitudes = map(x -> x.lat.val, Coordinates)
+    longitudes = map(x -> x.lon.val, Coordinates)
+
+    return SurfaceVelocityData(x=x, y=y, lat=latitudes, lon=longitudes, vx=vx, vy=vy, vabs=vabs, vx_error=vx_error, vy_error=vy_error, vabs_error=vabs_error, date=date_mean, date1=date1, date2=date2, date_error=date_error)
+end
+
+"""
+    max_or_empyt(A::Array)
+
+Return maximum value for non-empty arrays. 
+This is just required to compute the error in the absolute velocity.
+"""
+function max_or_empyt(A::Array)
+    if length(A) == 0
+        return 0.0
+    else
+        return maximum(A)
+    end
+end

src/parameters/PhysicalParameters.jl

@@ -13,28 +13,29 @@ end
 
 """
     PhysicalParameters(;
-        ρ::Float64 = 900.0,
-        g::Float64 = 9.81,
-        ϵ::Float64 = 1e-3,
-        η₀::F = 1.0, 
-        maxA::Float64 = 8e-17,
-        minA::Float64 = 8.5e-20,
-        maxTlaw::Float64 = 1.0,
-        minTlaw::Float64 = -25.0,
-        noise_A_magnitude::Float64 = 5e-18
-        )
+            ρ::F = 900.0,
+            g::F = 9.81,
+            ϵ::F = 1e-3,
+            η₀::F = 1.0, 
+            maxA::F = 8e-17,
+            minA::F = 8.5e-20,
+            maxTlaw::F = 1.0,
+            minTlaw::F = -25.0,
+            noise_A_magnitude::F = 5e-18
+            )
+
 Initialize the physical parameters of a model.
 Keyword arguments
 =================
     - `ρ`: Ice density
     - `g`: Gravitational constant
-    - `n`: Glen's exponent
-    - `A`: Glen's coefficient
     - `ϵ`: Small number
-    - `C`: Sliding coefficient
     - `η₀`:  
     - `maxA`: Maximum value for `A` (Glen's coefficient)
     - `minA`: Minimum value for `A` (Glen's coefficient)
+    - `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
 """
 function PhysicalParameters(;
             ρ::F = 900.0,