Code for preprocessing Rabatel (2023) data #71
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,5 @@ | ||
| export AbstractData | ||
|
|
||
| abstract type AbstractData end | ||
|
|
||
| include("SurfaceVelocityData.jl") |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,75 @@ | ||
| export SurfaceVelocityData | ||
|
|
||
| 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} | ||
| date::Union{Vector{DateTime}, Nothing} | ||
| date1::Union{Vector{DateTime}, Nothing} | ||
| date2::Union{Vector{DateTime}, Nothing} | ||
| date_error::Union{Vector{Day}, Vector{Millisecond}, Nothing} | ||
| 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}, Vector{Millisecond}, 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}, Vector{Millisecond}, 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") | ||
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,109 @@ | ||
| 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, compute_vabs_error::Bool=true) | ||
|
|
||
| # 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") | ||
| # Extract string from metadata of when days start counting | ||
| # Default format of attribute "units" is "days since 2015-08-14 00:00:00" | ||
| date_mean_since = ncgetatt(file, "mid_date", "units")[12:21] # e.g., "2015-08-14" | ||
| date1_offset_since = ncgetatt(file, "date1", "units")[12:21] # e.g., "2015-07-30" | ||
| date2_offset_since = ncgetatt(file, "date2", "units")[12:21] # e.g., "2015-08-29" | ||
| # Convertion to Julia datetime | ||
| date_mean_offset = datetime2julian(DateTime(date_mean_since)) - 2400000.5 | ||
|
There was a problem hiding this comment. Why do we still need to convert to Julian time? This was used before when we had to interact with Python dates, but this should no longer be the case, right? There was a problem hiding this comment. I am not sure what you are referring @JordiBolibar , but this transformation is required since that is the original format of the dataset. Days are count from customize starting dates, so this is necessary to set all in the same time reference. Anyways... this is just for the intepolated dataset, I would not pay much attention to this for now. |
||
| date1_offset = datetime2julian(DateTime(date1_offset_since)) - 2400000.5 | ||
| date2_offset = datetime2julian(DateTime(date_mean_offset)) - 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 | ||
| if compute_vabs_error | ||
| vx_ratio_max = map(i -> max_or_empty(abs.(vx[:,:,i][vabs[:,:,i] .> 0.0]) ./ vabs[:,:,i][vabs[:,:,i] .> 0.0]), 1:size(vx)[3]) | ||
| vy_ratio_max = map(i -> max_or_empty(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 | ||
| vabs_error = nothing | ||
| end | ||
| 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 | ||
albangossard marked this conversation as resolved.
Show resolved
Hide resolved
|
||
|
|
||
| # 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_empty(A::Array) | ||
| if length(A) == 0 | ||
| return 0.0 | ||
| else | ||
| return maximum(A) | ||
| end | ||
| end | ||
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Maybe we could add an assert on the value of
latto make sure that the user doesn't use data from the southern hemisphere, don't you think?There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
mmm my understanding is that the projection is defined with the zone and the hemisphere information. A priori, you don't have the latitude, you need to compute it with the projection knowing the hemisphere. Maybe I am wrong, but I don't think there is big danger in leaving it as it is.