Moving some exploratory code here

Julia/CondaPkg.toml

@@ -0,0 +1,7 @@
+[deps]
+cartopy = ""
+pandas = ""
+matplotlib = ""
+xarray = ""
+seaborn = ""
+numpy = ""

Julia/examples/explore.jl

@@ -0,0 +1,91 @@
+using Pkg; Pkg.activate(".")
+using ExploreODINN
+
+# Exploration of Glacier data is done with OGGM: https://tutorials.oggm.org/stable/notebooks/tutorials/working_with_rgi.html
+# As there explained, we can find an individual glacier using the GLIMS glacier viewer: https://www.glims.org/maps/glims
+
+
+# Data corresponding to Mont Blanc Massif
+# Raw data
+nc_file_multi = "../../data/DATASET_FACU/c_x01225_y03675_all_filt-multi.nc"
+# Interpolated dataset 
+nc_file_interp = "../../data/DATASET_FACU/c_x01225_y03675_interp.nc"
+
+# Read the dataset
+data = preprocessing(nc_file_multi)
+
+# Let's make a plot with the timeseries for a grid point with good measurement
+
+# Bossons Glacier
+# idx, idy = 165, 95
+# Brenva Glacier 
+idx, idy = 225, 235
+
+plot_timeseries(data, idx, idy; ignore_zeros=true, saveas="figures/timeseries_vabs.pdf")
+plot_timeseries(data, idx-1, idy+1; ignore_zeros=true, saveas="figures/timeseries_vabs_2.pdf")
+plot_timeseries(data, idx+1, idy-1; ignore_zeros=true, saveas="figures/timeseries_vabs_3.pdf")
+
+
+# Plot the total number of observations
+plot_count(data, saveas="figures/v_count.pdf");
+
+
+# facubsapienza
+# FacuEarthData_2025
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+# id_date = 70
+# @show date1[id_date]
+
+
+# zs1 = vx[:, :, id_date]
+# zs2 = vy[:, :, id_date]
+# zs3 = v_abs[:, :, id_date]
+
+# vmax = maximum(x -> isnan(x) ? -Inf : x, zs3)
+
+# # zs3 = replace!(x -> x==0.0 ? NaN : x, Array{Union{Float64, Nothing}}(zs3))
+# zs1 = replace!(x -> x==0.0 ? NaN : x, zs1)
+# zs2 = replace!(x -> x==0.0 ? NaN : x, zs2)
+# zs3 = replace!(x -> x==0.0 ? NaN : x, zs3)
+
+# joint_limits = (-vmax, vmax)  # here we pick the limits manually for simplicity instead of computing them
+# abs_limits = (0, vmax)
+
+# fig = Figure(size=(2000, 600), axis=(;title="dsds"));
+
+# ax1, hm1 = heatmap(fig[1,1], xs, ys, zs1,  colorrange = abs_limits, colormap = Reverse(:acton25), nan_color = :transparent, axis=(; title="V Abs (Data: $(date1[id_date]))"));
+# ax2, hm2 = heatmap(fig[1, end+1], xs, ys, zs2, colorrange = joint_limits, colormap = :broc25, axis=(; title="Vx"));
+# ax3, hm3 = heatmap(fig[1, end+1], xs, ys, zs3, colorrange = joint_limits, colormap = :broc25, axis=(; title="Vy"));
+
+# # CairoMakie.scatter!(ax1, Point(xs[idx], ys[idy]); marker='🦜', markersize=40, color = :purple)
+# CairoMakie.scatter!(ax1, Point(xs[idx], ys[idy]); marker='x', markersize=40, color = :purple)
+
+# # Colorbar(fig[:, end+1], hm1)                     # These three
+# # Colorbar(fig[:, end+1], hm2)                     # colorbars are
+# Colorbar(fig[:, end+1], hm1);                      # colorbars are
+# Colorbar(fig[:, end+1], hm2);  # equivalent
+
+# save("figures/v_heatmap.pdf", fig)
+
+

Julia/examples/geodesy.jl

@@ -0,0 +1,83 @@
+
+
+# Playing with coodinate projection 
+
+using Geodesy
+
+projection = LLAfromUTMZ(wgs84)
+
+x_bossom = UTMZ(-682.471, 5.1912e6, 0.0, 31,true)
+
+tm = Geodesy.TransverseMercator(WGS84)
+
+x₀ = 10.1335
+x = -682.471
+y = 5.1912e6
+
+Geodesy.transverse_mercator_reverse(x₀, x, y , 0.9996, tm)
+
+
+###
+using Unitful: m, rad, °
+using CoordRefSystems
+
+k₀ = 0.9996
+
+# xₒ = 500000.0m
+xₒ = 0.0m
+yₒ = 0.0m
+# yₒ = hemisphere == :north ? 0.0m : 10000000.0m
+
+lonₒ = 8.01919°
+# lonₒ = 10.1335°
+latₒ = 0.0°
+
+
+S = CoordRefSystems.Shift(; lonₒ, xₒ, yₒ)
+
+x = -8040.94
+y = 5.15794e6
+x_glacier = TransverseMercator{k₀,latₒ,datum,S}(x, y)
+
+convert(LatLon, x_glacier)
+
+# Clean for each glacier
+
+k₀ = 0.9996
+xₒ = 0.0m
+yₒ = 0.0m
+
+# Aletschgletscher
+
+lonₒ = 8.01919°
+latₒ = 0.0°
+x = -8040.94
+y = 5.15794e6
+S = CoordRefSystems.Shift(; lonₒ, xₒ, yₒ)
+x_glacier = TransverseMercator{k₀,latₒ,WGS84Latest,S}(x, y)
+
+convert(LatLon, x_glacier)
+# Answer: Lat = 46.574977403312715°, Lon = 7.914252983649885° (correct for upper left corner)
+
+# Bossons 
+lonₒ = 10.1335°
+latₒ = 0.0°
+x = -682.471
+y = 5.1912e6
+S = CoordRefSystems.Shift(; lonₒ, xₒ, yₒ)
+x_glacier = TransverseMercator{k₀,latₒ,WGS84Latest,S}(x, y)
+
+convert(LatLon, x_glacier)
+# Answer: Lat = 46.874338533298065°, Lon = 10.124544115391744°
+
+
+# Glacier d’Argentière
+lonₒ = 6.985°
+latₒ = 0.0°
+x = -3574.0
+y = 5.09221e6
+S = CoordRefSystems.Shift(; lonₒ, xₒ, yₒ)
+x_glacier = TransverseMercator{k₀,latₒ,WGS84Latest,S}(x, y)
+
+convert(LatLon, x_glacier)
+# Answer: Lat = 45.98345259928449°, Lon = 6.938857312258543

Julia/examples/glacier_oggm.jl

@@ -0,0 +1,125 @@
+# Let's retrieve some more data from OGGM for one of the glaciers!
+using Pkg; Pkg.activate(".")
+
+
+using ExploreODINN
+# include("../src/ExploreODINN.jl")
+using Sleipnir, Huginn
+using Plots
+using Plots.PlotMeasures
+using MakieCore: heatmap
+using CairoMakie
+using Statistics
+using Rasters
+using Dates
+using PythonCall
+
+working_dir = joinpath(homedir(), ".ODINN/ODINN_prepro")
+
+# Bossons glacier
+rgi_ids = ["RGI60-11.03646"] 
+
+# Aletschgletscher
+# rgi_ids = ["RGI60-11.01450"] 
+
+# Glacier d’Argentière
+# rgi_ids = ["RGI60-11.03638"] 
+
+# Glacier Mar de Glace
+# rgi_ids = ["RGI60-11.03643"] 
+
+# rgi_ids = ["RGI60-11.00872"]
+
+# Perito Moreno 
+# rgi_ids = ["RGI60-17.00312"] 
+
+# Fix this
+rgi_paths = get_rgi_paths()
+# rgi_paths = Dict("RGI60-11.03643" => "per_glacier/RGI60-11/RGI60-11.03/RGI60-11.03643")
+rgi_paths = Dict(k => rgi_paths[k] for k in rgi_ids)
+
+params = Huginn.Parameters(simulation = SimulationParameters(use_MB=false,
+                                                          velocities=false,
+                                                          tspan=(2010.0, 2015.0),
+                                                          working_dir = working_dir,
+                                                          multiprocessing = false,
+                                                        #   working_dir = Huginn.root_dir,
+                                                          test_mode = true,
+                                                          rgi_paths = rgi_paths),
+                        solver = SolverParameters(reltol=1e-12)
+                        )
+
+# Define new model
+model = Huginn.Model(iceflow = SIA2Dmodel(params), mass_balance = nothing)
+
+# We retrieve some glaciers for the simulation
+glaciers = initialize_glaciers(rgi_ids, params)
+
+# Let's grab a single glacier
+glacier = glaciers[1]
+
+lats = glacier.Coords["lat"]
+lons = glacier.Coords["lon"]
+
+# Ice Surface Velocity Data
+
+nc_file_multi = "../../data/DATASET_FACU/c_x01225_y03675_all_filt-multi.nc"
+ice_surface_data = preprocessing(nc_file_multi)
+
+v_count = mean((ice_surface_data.vabs .!== 0.0) .* (.!isnan.(ice_surface_data.vabs)), dims=3)[:,:,1]
+
+# Python Plot
+
+fig, axes = plt[].subplots()
+fig.set_size_inches(10, 10)
+caxes = axes.inset_axes([1.04, 0.06, 0.03, 0.4])
+
+# axes.set_aspect("equal")
+axes.set_xlabel("Longitude")
+axes.set_ylabel("Latitude")
+axes.spines[pylist(["right", "top"])].set_visible(false)
+
+xs = lons #range(1, 48, length=48)
+ys = lats # range(1, 79, length=79)
+xs_mesh = xs' .* ones(length(ys))
+ys_mesh = ones(length(xs))' .* ys
+
+Δdist = 0.4 * (glacier.Δx^2 + glacier.Δy^2)^0.5
+
+H_masked = copy(glacier.H₀)
+H_masked[H_masked .< 0.01] .= NaN
+
+# V_abs = (glacier.Vx.^2 + glacier.Vy.^2).^0.5
+# V_abs[glacier.H₀ .< 0.01] .= NaN
+
+# H_masked = np[].ma.array(glacier.H₀', mask=0.0)
+
+# ColorGrid = axes.pcolormesh(xs_mesh, ys_mesh, H_masked', alpha=0.20)
+# ColorGrid = axes.pcolormesh(xs_mesh, ys_mesh, H_masked', alpha=0.20)
+ColorGrid = axes.pcolormesh(ice_surface_data.lon, ice_surface_data.lat, v_count, alpha=0.70)
+# axes.scatter(ice_surface_data.lon', ice_surface_data.lat', s=1.0, c="k")
+# ColorGrid = axes.pcolormesh(xs_mesh, ys_mesh, V_abs', alpha=0.95)
+
+ContourLines = axes.contour(xs_mesh, ys_mesh, H_masked', 20, colors='k', levels=[Δdist, 50, 100, 150, 200, 250])  
+
+cbar = plt[].colorbar(ColorGrid, cax=caxes, orientation="vertical")
+
+plt[].savefig("./figures/thickness.pdf", dpi=300, format="pdf", bbox_inches="tight")
+
+
+# Temperature plot
+plot_temp = Plots.plot(glacier.climate.raw_climate[At(DateTime(2012):Day(1):DateTime(2019))][:temp])
+Plots.savefig(plot_temp, "./figures/temperature.pdf")
+# plt[].savefig("./figures/temperature.pdf", dpi=300, format="pdf", bbox_inches="tight")
+
+# Precipitation plot
+plot_prcp = Plots.plot(glacier.climate.raw_climate[At(DateTime(2012):Day(1):DateTime(2019))][:prcp])
+# plt[].savefig("./figures/precipitation.pdf", dpi=300, format="pdf", bbox_inches="tight")
+Plots.savefig(plot_prcp, "./figures/precipitation.pdf")
+
+
+
+
+
+# Now, let's plot the contours of the glacier on top of an image 
+