Calculation of Transversal Mercator Projection for glacier coordinates

Project.toml

@@ -8,6 +8,7 @@ CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 CodecZlib = "944b1d66-785c-5afd-91f1-9de20f533193"
 ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
+CoordRefSystems = "b46f11dc-f210-4604-bfba-323c1ec968cb"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
@@ -22,6 +23,7 @@ Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Tar = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 CSV = "0.10.15"

src/Sleipnir.jl

@@ -21,6 +21,8 @@ using JSON
 using CodecZlib
 using Tar
 import NCDatasets
+using Unitful: m, rad, °
+using CoordRefSystems
 
 # ##############################################
 # ############    PARAMETERS     ###############

src/glaciers/climate/climate2D_utils.jl

@@ -155,8 +155,8 @@ function downscale_2D_climate(climate_step::Dict, glacier::Glacier2D)
                        rain=rain_2D,
                        gradient=Float64(climate_step["gradient"]),
                        avg_gradient=Float64(climate_step["avg_gradient"]),
-                       x=glacier.S_coords["x"],
-                       y=glacier.S_coords["y"],
+                       x=glacier.Coords["lon"],
+                       y=glacier.Coords["lat"],
                        ref_hgt=Float64(climate_step["ref_hgt"]))
 
     # Apply temperature gradients and compute snow/rain fraction for the selected period

src/glaciers/glacier/Glacier1D.jl

@@ -19,7 +19,7 @@ mutable struct Glacier1D{F <: AbstractFloat, I <: Integer} <: AbstractGlacier
     λ::Union{Vector{F}, Nothing}
     slope::Union{Vector{F}, Nothing}
     dist_border::Union{Vector{F}, Nothing}
-    S_coords::Union{Dict{String, Vector{Float64}}, Nothing}
+    Coords::Union{Dict{String, Vector{Float64}}, Nothing}
     Δx::Union{F, Nothing}
     Δy::Union{F, Nothing}
     nx::Union{I, Nothing}
@@ -41,7 +41,7 @@ function Glacier1D(;
     λ::Union{Vector{F}, Nothing} = nothing,
     slope::Union{Vector{F}, Nothing} = nothing,
     dist_border::Union{Vector{F}, Nothing} = nothing,
-    S_coords::Union{Dict{String, Vector{Float64}}, Nothing} = nothing,
+    Coords::Union{Dict{String, Vector{Float64}}, Nothing} = nothing,
     Δx::Union{F, Nothing} = nothing,
     Δy::Union{F, Nothing} = nothing,
     nx::Union{I, Nothing} = nothing,
@@ -64,7 +64,7 @@ function Glacier1D(;
     λ::Union{Vector{F}, Nothing} = nothing,
     slope::Union{Vector{F}, Nothing} = nothing,
     dist_border::Union{Vector{F}, Nothing} = nothing,
-    S_coords::Union{Dict{String, Vector{Float64}}, Nothing} = nothing,
+    Coords::Union{Dict{String, Vector{Float64}}, Nothing} = nothing,
     Δx::Union{F, Nothing} = nothing,
     Δy::Union{F, Nothing} = nothing,
     nx::Union{I, Nothing} = nothing,
@@ -74,7 +74,7 @@ function Glacier1D(;
     # Define default float and integer type for constructor
     ft = Float64
     it = Int64
-    return Glacier1D{ft,it}(rgi_id, gdir, climate, H₀, S, B, V, A, C, n, w₀, λ, slope, dist_border, S_coords, Δx, Δy, nx, ny)
+    return Glacier1D{ft,it}(rgi_id, gdir, climate, H₀, S, B, V, A, C, n, w₀, λ, slope, dist_border, Coords, Δx, Δy, nx, ny)
 end
 
 ###############################################
@@ -85,7 +85,7 @@ Base.:(==)(a::Glacier1D, b::Glacier1D) = a.rgi_id == b.rgi_id && a.gdir == b.gdi
                                       a.H₀ == b.H₀ && a.S == b.S && a.B == b.B && a.V == b.V &&
                                       a.A == b.A && a.C == b.C && a.n == b.n && a.w₀ == b.w₀ && a.λ == b.λ &&
                                       a.slope == b.slope && a.dist_border == b.dist_border && a.rgi_id == b.rgi_id &&
-                                      a.S_coords == b.S_coords && a.Δx == b.Δx && a.Δy == b.Δy && a.Δx == b.Δx && a.nx == b.nx && a.ny == b.ny
+                                      a.Coords == b.Coords && a.Δx == b.Δx && a.Δy == b.Δy && a.Δx == b.Δx && a.nx == b.nx && a.ny == b.ny
 
 include("glacier1D_utils.jl")
 include("../climate/climate1D_utils.jl")

src/glaciers/glacier/Glacier2D.jl

@@ -20,7 +20,7 @@ mutable struct Glacier2D{F <: AbstractFloat, I <: Integer} <: AbstractGlacier
     n::Union{F, Nothing}
     slope::Union{Matrix{F}, Nothing}
     dist_border::Union{Matrix{F}, Nothing}
-    S_coords::Union{Dict{String, Vector{Float64}}, Nothing}
+    Coords::Union{Dict{String, Vector{Float64}}, Nothing}
     Δx::Union{F, Nothing}
     Δy::Union{F, Nothing}
     nx::Union{I, Nothing}
@@ -45,7 +45,7 @@ function Glacier2D(;
     n::Union{F, Nothing} = nothing,
     slope::Union{Matrix{F}, Nothing} = nothing,
     dist_border::Union{Matrix{F}, Nothing} = nothing,
-    S_coords::Union{Dict{String, Vector{Float64}}, Nothing} = nothing,
+    Coords::Union{Dict{String, Vector{Float64}}, Nothing} = nothing,
     Δx::Union{F, Nothing} = nothing,
     Δy::Union{F, Nothing} = nothing,
     nx::Union{I, Nothing} = nothing,
@@ -71,7 +71,7 @@ function Glacier2D(;
     n::Union{F, Nothing} = nothing,
     slope::Union{Matrix{F}, Nothing} = nothing,
     dist_border::Union{Matrix{F}, Nothing} = nothing,
-    S_coords::Union{Dict{String, Vector{Float64}}, Nothing} = nothing,
+    Coords::Union{Dict{String, Vector{Float64}}, Nothing} = nothing,
     Δx::Union{F, Nothing} = nothing,
     Δy::Union{F, Nothing} = nothing,
     nx::Union{I, Nothing} = nothing,
@@ -84,7 +84,7 @@ function Glacier2D(;
     ft = typeof(Δx)
     it = typeof(nx)
 
-    return Glacier2D{ft,it}(rgi_id, climate, H₀, H_glathida, S, B, V, Vx, Vy, A, C, n, slope, dist_border, S_coords, Δx, Δy, nx, ny, cenlon, cenlat)
+    return Glacier2D{ft,it}(rgi_id, climate, H₀, H_glathida, S, B, V, Vx, Vy, A, C, n, slope, dist_border, Coords, Δx, Δy, nx, ny, cenlon, cenlat)
 end
 
 ###############################################
@@ -96,7 +96,7 @@ Base.:(==)(a::Glacier2D, b::Glacier2D) = a.rgi_id == b.rgi_id && a.climate == b.
                                       a.H₀ == b.H₀ && a.H_glathida == b.H_glathida && a.S == b.S && a.B == b.B && a.V == b.V &&
                                       a.A == b.A && a.C == b.C && a.n == b.n &&
                                       a.slope == b.slope && a.dist_border == b.dist_border &&
-                                      a.S_coords == b.S_coords && a.Δx == b.Δx && a.Δy == b.Δy && a.nx == b.nx && a.ny == b.ny &&
+                                      a.Coords == b.Coords && a.Δx == b.Δx && a.Δy == b.Δy && a.nx == b.nx && a.ny == b.ny &&
                                       a.cenlon == b.cenlon && a.cenlat == b.cenlat
 
 
@@ -105,9 +105,25 @@ Base.:(≈)(a::Glacier2D, b::Glacier2D) = a.rgi_id == b.rgi_id && a.climate == b
                                         safe_approx(a.S, b.S) && safe_approx(a.B, b.B) && safe_approx(a.V, b.V) &&
                                         safe_approx(a.A, b.A) && safe_approx(a.C, b.C) && safe_approx(a.n, b.n) &&
                                         isapprox(a.slope, b.slope; rtol=1e-3) && safe_approx(a.dist_border, b.dist_border) &&
-                                        safe_approx(a.S_coords, b.S_coords) && safe_approx(a.Δx, b.Δx) && safe_approx(a.Δy, b.Δy) &&
+                                        safe_approx(a.Coords, b.Coords) && safe_approx(a.Δx, b.Δx) && safe_approx(a.Δy, b.Δy) &&
                                         safe_approx(a.nx, b.nx) && safe_approx(a.ny, b.ny) &&
                                         safe_approx(a.cenlon, b.cenlon) && safe_approx(a.cenlat, b.cenlat)
 
+# Display setup
+function Base.show(io::IO, glacier::Glacier2D)
+    println(io) 
+    println("Glacier:  $(glacier.rgi_id)")
+    println("Available fields:")
+    for key in fieldnames(Glacier2D)
+        value = getproperty(glacier, key)
+        if !isnothing(value)
+            println("  * ", string(key))
+        end
+    end
+end
+# Vectorial form 
+Base.show(io::IO, glaciers::Vector{Glacier2D}) = Base.show.(Ref(io), glaciers)
+
+
 include("glacier2D_utils.jl")
 include("../climate/climate2D_utils.jl")

src/glaciers/glacier/glacier2D_utils.jl

@@ -32,12 +32,15 @@ function initialize_glaciers(rgi_ids::Vector{String}, params::Parameters; test=f
     # Generate raw climate data if necessary
     if params.simulation.test_mode
         map((rgi_id) -> generate_raw_climate_files(rgi_id, params.simulation), rgi_ids) # avoid GitHub CI issue
+        # @infiltrate
+        # glaciers::Vector{Glacier2D} = map((rgi_id) -> initialize_glacier(rgi_id, params; smoothing=false, test=test), rgi_ids)
     else
         pmap((rgi_id) -> generate_raw_climate_files(rgi_id, params.simulation), rgi_ids)
+        # glaciers::Vector{Glacier2D} = pmap((rgi_id) -> initialize_glacier(rgi_id, params; smoothing=false, test=test), rgi_ids)
     end
-
+        
     glaciers::Vector{Glacier2D} = pmap((rgi_id) -> initialize_glacier(rgi_id, params; smoothing=false, test=test), rgi_ids)
-    
+
     if params.simulation.use_glathida_data == true
 
         # Obtain H_glathida values for the valid RGI IDs
@@ -122,9 +125,18 @@ function initialize_glacier_data(rgi_id::String, params::Parameters; smoothing=f
             # H_mask = (dist_border .< 20.0) .&& (S .> maximum(S)*0.7)
             # H₀[H_mask] .= 0.0
 
+        # Mercator Projection 
+        params_projection = parse_proj(glacier_grid["proj"])
+        transform(X,Y) = UTMercantor(X, Y; k=params_projection["k"], cenlon=params_projection["lon_0"], cenlat=params_projection["lat_0"], 
+                        x0=params_projection["x_0"], y0=params_projection["y_0"])
+        easting = dims(glacier_gd, 1).val
+        northing = dims(glacier_gd, 2).val
+        latitudes = map(x -> x.lat.val, transform.(Ref(mean(easting)), northing))
+        longitudes = map(x -> x.lon.val, transform.(easting, Ref(mean(northing))))
+
         B = glacier_gd.topo.data .- H₀ # bedrock
 
-        S_coords = Dict{String,Vector{Float64}}("x"=> dims(glacier_gd, 1).val, "y"=> dims(glacier_gd, 2).val)
+        Coords = Dict{String,Vector{Float64}}("lon"=> longitudes, "lat"=> latitudes)
         S::Matrix{Float64} = glacier_gd.topo.data
         #smooth!(S)
 
@@ -152,7 +164,7 @@ function initialize_glacier_data(rgi_id::String, params::Parameters; smoothing=f
                           H₀ = H₀, S = S, B = B, V = V, Vx = Vx, Vy = Vy,
                           A = 4e-17, C = 0.0, n = 3.0,
                           slope = slope, dist_border = dist_border,
-                          S_coords = S_coords, Δx=Δx, Δy=Δy, nx=nx, ny=ny,
+                          Coords = Coords, Δx=Δx, Δy=Δy, nx=nx, ny=ny,
                           cenlon = glacier_grid["x0y0"][1] , cenlat = glacier_grid["x0y0"][2])
         return glacier
 
@@ -313,6 +325,44 @@ function fillZeros(A, fill=NaN)
     return @. ifelse(iszero(A), fill, A)
 end
 
+"""
+    parse_proj(proj::String)
+
+Parses the string containing the information of the projection to filter for important information
+"+proj=tmerc +lat_0=0 +lon_0=6.985 +k=0.9996 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
+"""
+function parse_proj(proj::String)
+    res = Dict()
+    ℓ = split(proj, (' ', '+', '='))
+    ℓ = ℓ[ℓ .!= ""]
+    for (i, key) in enumerate(ℓ)
+        if key ∈ ["lat_0", "lon_0", "k", "x_0", "y_0"]
+            res[key] = parse(Float64, ℓ[i+1])
+        end
+    end
+    return res
+end
+
+"""
+
+Transverse Mercantor Projection
+"""
+function UTMercantor(x::F, y::F; k=0.9996, cenlon=0.0, cenlat=0.0, x0=0.0, y0=0.0) where {F <: AbstractFloat}
+
+    # Convert to right units 
+    lonₒ = cenlon * 1.0°
+    latₒ = cenlat * 1.0°
+    xₒ = x0 * 1.0m
+    yₒ = y0 * 1.0m
+
+    # Define shift in new coordinate system
+    S = CoordRefSystems.Shift(; lonₒ, xₒ, yₒ)
+
+    datum = WGS84Latest
+    location = TransverseMercator{k, latₒ, datum, S}(x, y)
+
+    return convert(LatLon, location)
+end
 """
     smooth!(A)