Skip to content

Ice velocity integration #97

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 16 commits into from
Apr 25, 2025
Merged

Ice velocity integration #97

Show file tree
Hide file tree
Changes from 14 commits
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
ed4f0d8
fix error when keyword arguments are not provided in Results
albangossard Apr 9, 2025
b3b0015
update jld2 files
albangossard Apr 9, 2025
5633df2
test Results instantiation in unit tests
albangossard Apr 9, 2025
05692d4
trying to fix CI error
albangossard Apr 9, 2025
e2c4f9e
fix processVelocity type issue
albangossard Apr 10, 2025
c278657
put back plot tests
albangossard Apr 10, 2025
43e1f55
preprocessed data is redownloaded every day + add name field to Glaci…
albangossard Apr 15, 2025
084f982
add feature to grid surface velocity to the glacier grid + generate f…
albangossard Apr 16, 2025
4613425
fix tests + add docstrings
albangossard Apr 17, 2025
1f5230f
add thickness and velocity data to the glacier struct
albangossard Apr 17, 2025
1c326f6
fix unit test
albangossard Apr 17, 2025
c9b8221
improve plotting functions
albangossard Apr 18, 2025
9c1c119
finish updating the plotting functions
albangossard Apr 23, 2025
3a50573
fix unit tests
albangossard Apr 23, 2025
a270f5f
changes following PR reviews
albangossard Apr 24, 2025
c8a1ec4
bump version and update compat entries
albangossard Apr 25, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@ CoordRefSystems = "b46f11dc-f210-4604-bfba-323c1ec968cb"
DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
DateFormats = "44557152-fe0a-4de1-8405-416d90313ce6"
Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0"
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is already available from Rasters.jl. Did you add it just to have the necessary exported namespaces?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes this is just to access Dim in Sleipnir/test/surface_velocity.jl for example. If that bothers you we can probably do Rasters.Dim instead.

Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
Expand All @@ -23,7 +24,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"
Observables = "510215fc-4207-5dde-b226-833fc4488ee2"
Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
Rasters = "a3a2b9e3-a471-40c9-b274-f788e487c689"
Expand Down Expand Up @@ -51,7 +52,6 @@ Infiltrator = "1.8"
JLD2 = "0.4, 0.5"
JSON = "0.21.4"
NCDatasets = "0.14.6"
NetCDF = "0.12"
Rasters = "0.13.0, 0.14"
Revise = "3"
Statistics = "1"
Expand Down
17 changes: 14 additions & 3 deletions src/Sleipnir.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ using JLD2
using Distributed
using Statistics
using CairoMakie
using Observables
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is this for? I've seen it's some sort of improved Ref.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is used in lines like

Observables.connect!(cb.height, @lift CairoMakie.Fixed($(viewport(ax.scene)).widths[2]))

whose purpose is to make the colobar have the same height as the heatmap we are plotting.

using Downloads
using HDF5
using ComponentArrays
Expand All @@ -24,7 +25,6 @@ import NCDatasets
using Unitful: m, rad, °
using CoordRefSystems
using Dates, DateFormats
using NetCDF
using GR
using DataStructures
using ImageInTerminal
Expand All @@ -50,14 +50,25 @@ end
# ##############################################

include("setup/config.jl")
# All parameters needed for the models
include("parameters/Parameters.jl")

# Anything related to managing glacier data used for data assimilation
include("glaciers/data/Data.jl")

# Anything related to managing glacier topographical and climate variables
include("glaciers/glacier/Glacier.jl")

# All parameters needed for the models
include("parameters/Parameters.jl")

# The utils of surface velocity data, glaciers and climate need the struct to be already
# defined since they depend on each other. This is why we import them afterwards
include("glaciers/data/SurfaceVelocityData_utils.jl")
include("glaciers/glacier/glacier2D_utils.jl")
include("glaciers/climate/climate2D_utils.jl")

# All structures and functions related to ODINN models
include("models/Model.jl")

# Everything related to running simulations in ODINN
include("simulations/Simulation.jl")

Expand Down
8 changes: 7 additions & 1 deletion src/glaciers/data/Data.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,12 @@
export AbstractData

"""
AbstractData

Abstract type that represents data. Used to implement `ThicknessData` and `SurfaceVelocityData`.
"""
abstract type AbstractData end

include("ThicknessData.jl")
include("SurfaceVelocityData.jl")
include("SurfaceVelocityMapping.jl")
include("SurfaceVelocityData.jl")
92 changes: 53 additions & 39 deletions src/glaciers/data/SurfaceVelocityData.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,36 +7,38 @@ 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{Array{F, 3}, Nothing}`: x component of surface velocity.
- `vy::Union{Array{F, 3}, Nothing}`: y component of surface velocity.
- `vabs::Union{Array{F, 3}, Nothing}`: Absolute ice surface velocity.
- `vx_error::Union{Array{F, 1}, Nothing}`: Error in `vx`
- `vy_error::Union{Array{F, 1}, Nothing}`: Error in `vy`
- `vabs_error::Union{Array{F, 1}, 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 adquisition.
- `date2::Union{Vector{DateTime}, Nothing}`: Second date of adquisition.
- `date_error::Union{Vector{Day}, Vector{Millisecond}, Nothing}`: Error in `date`.
- `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 component of surface velocity.
- `vy::Union{Vector{Matrix{F}}, Nothing}`: y component of surface velocity.
- `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 adquisition.
- `date2::Union{Vector{DateTime}, Nothing}`: Second date of adquisition.
- `date_error::Union{Vector{Day}, Vector{Millisecond}, Nothing}`: Error in `date`.
- `glacierGridded::Bool`: Whether the data have been gridded to the glacier grid or not.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's not entirely clear what you mean by this. Is there a case where the date is not gridded to the glacier grid? Like, do you use the same struct to store the original raw data and the final data in the glacier grid?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes we use the same struct to store the original raw data and the gridded ones. I don't think we need a specific struct to store the original raw data. Doing so will even come with many memory issues as here we load the data in a lazy mode only when we need them and many of the points in the original data aren't useful since they are outside of the glacier grid.
Most of the time, we probably won't use the raw data directly. However, it can be convenient in some cases for debugging or to visualize data of the netcdf files.

I changed the name to isGridGlacierAligned to make it clearer.

"""
mutable struct SurfaceVelocityData{F <: AbstractFloat} <: AbstractData
x::Union{Vector{F}, Nothing}
y::Union{Vector{F}, Nothing}
lat::Union{Vector{F}, Nothing}
lon::Union{Vector{F}, Nothing}
vx::Union{Array{F, 3}, Nothing}
vy::Union{Array{F, 3}, Nothing}
vabs::Union{Array{F, 3}, Nothing}
vx_error::Union{Array{F, 1}, Nothing}
vy_error::Union{Array{F, 1}, Nothing}
vabs_error::Union{Array{F, 1}, Nothing}
vx::Union{Vector{Matrix{F}}, Nothing}
vy::Union{Vector{Matrix{F}}, Nothing}
vabs::Union{Vector{Matrix{F}}, Nothing}
Comment on lines +31 to +33
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks good to me! Just curious: any reason why a vector of matrices is better than just a 3D tensor?

Copy link
Member

@JordiBolibar JordiBolibar Apr 24, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it's nicer to handle. Since two dimensions are spatial and one is temporal, I feel it's much nicer to iterate through time and get back spatial grids. Otherwise the user needs to know which of the 3 dimensions is time.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

At the origin it was because the plotting functions assume that spatio-temporal data are a vector of matrix. I don't have a strong opinion except that it should be uniform across the different structs.

vx_error::Union{Vector{F}, Nothing}
vy_error::Union{Vector{F}, Nothing}
vabs_error::Union{Vector{F}, Nothing}
date::Union{Vector{DateTime}, Nothing}
date1::Union{Vector{DateTime}, Nothing}
date2::Union{Vector{DateTime}, Nothing}
date_error::Union{Vector{Day}, Vector{Millisecond}, Nothing}
glacierGridded::Bool
end

"""
Expand All @@ -48,24 +50,23 @@ function SurfaceVelocityData(;
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,
vx::Union{Vector{Matrix{F}}, Nothing} = nothing,
vy::Union{Vector{Matrix{F}}, Nothing} = nothing,
vabs::Union{Vector{Matrix{F}}, Nothing} = nothing,
vx_error::Union{Vector{F}, Nothing} = nothing,
vy_error::Union{Vector{F}, Nothing} = nothing,
vabs_error::Union{Vector{F}, 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}, Vector{Millisecond}, Nothing} = nothing,
) where {F <: AbstractFloat}
glacierGridded::Bool = false,
) 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.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is this? Does this come from the paper(s)? And by oversestimated you mean that the error is larger than reported or lower than reported?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also, this agreement about the northern hemisphere is quite strong. At the very least I would give a warning somehow (e.g. in case there is a way to infer that we are not in the northern hemisphere) and clearly document this in the docs.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Concerning the errors, this is something that @facusapienza21 implemented, so I'm not 100% sure of the answer. I guess this is because for every point of the grid, we don't have the velocities for every time and we use an upper bound by taking for every grid point the maximum speed over time. In the end the error is overestimated. @facusapienza21 am I right?

As for the northern hemisphere restriction, actually this is something that I removed in this PR and the code now supports data in the southern hemisphere too.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't remember exactly what was the problem, but as long as you can specify if the projection is done in the north or south hemisphere it should be fine!

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think right now we can just do north hemisphere unless we change this code:

  # Spatial preprocessing
  proj_zone = ncgetatt(file, "mapping", "utm_zone_number")
  transform(X,Y) = Sleipnir.UTMercator(X, Y; zone=proj_zone, hemisphere=:north)

Right now this is done for north hemisphere. I think this was also something from the dataset: since the dataset for now is just for the north hemisphere, it is not as clear what will be the argument in the future from the dataset to determine if the datacube is north or south hemisphere. This is a good question for Romain/Antoine, but for now I will keep the comment until we are sure the preprocessing will work with South data.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, I suggest we open an issue to keep track of this.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In fact, I already changed that line in the PR. The code is now:

hemisphere = (y[end]+y[begin])/2 >= 0 ? :north : :south
transform(X,Y) = Sleipnir.UTMercator(X, Y; zone=Int(params_projection["zone"]), hemisphere=hemisphere)

Since in the UTM coordinate system, y corresponds to the northing, then y=0 corresponds to the equator (I double checked this numerically).

If you think there can be bugs, we can add an assert, put back the remark in the documentation and create an issue for that, but I think that this code will work with dataset in the south hemisphere.


Expand All @@ -80,26 +81,39 @@ function SurfaceVelocityData(;
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,
vx::Union{Vector{Matrix{F}}, Nothing} = nothing,
vy::Union{Vector{Matrix{F}}, Nothing} = nothing,
vabs::Union{Vector{Matrix{F}}, Nothing} = nothing,
vx_error::Union{Vector{F}, Nothing} = nothing,
vy_error::Union{Vector{F}, Nothing} = nothing,
vabs_error::Union{Vector{F}, 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}, Vector{Millisecond}, Nothing} = nothing,
) where {F <: AbstractFloat}
glacierGridded::Bool = false,
) where {F <: AbstractFloat}

ft = Sleipnir.Float

return SurfaceVelocityData{ft}(
x, y, lat, lon,
vx, vy, vabs, vx_error, vy_error, vabs_error,
date, date1, date2, date_error
date, date1, date2, date_error, glacierGridded
)
end


include("surfacevelocitydata_utils.jl")
Base.:(==)(a::SurfaceVelocityData, b::SurfaceVelocityData) =
a.x == b.x && a.y == b.y && a.lat == b.lat && a.lon == b.lon &&
a.vx == b.vx && a.vy == b.vy && a.vabs == b.vabs &&
a.vx_error == b.vx_error && a.vy_error == b.vy_error && a.vabs_error == b.vabs_error &&
a.date == b.date && a.date1 == b.date1 && a.date2 == b.date2 && a.date_error == b.date_error &&
a.glacierGridded == b.glacierGridded


Base.:(≈)(a::SurfaceVelocityData, b::SurfaceVelocityData) =
safe_approx(a.x, b.x) && safe_approx(a.y, b.y) && safe_approx(a.lat, b.lat) && safe_approx(a.lon, b.lon) &&
safe_approx(a.vx, b.vx) && safe_approx(a.vy, b.vy) && safe_approx(a.vabs, b.vabs) &&
safe_approx(a.vx_error, b.vx_error) && safe_approx(a.vy_error, b.vy_error) && safe_approx(a.vabs_error, b.vabs_error) &&
safe_approx(a.date, b.date) && safe_approx(a.date1, b.date1) && safe_approx(a.date2, b.date2) && safe_approx(a.date_error, b.date_error) &&
a.glacierGridded == b.glacierGridded
Loading
Loading