Matching coordinates of glaciers with OGGM data

[facusapienza21]

I am trying to retrieve the exact latitude and longitude of the different glaciers but I am having some degrees of latitude/longitude discrepancies and I am not sure why...

Right now, the glacier_grid (defined here) includes some information of how to compute the latitude/longitude based on UTC projections.

As it was explained by @fmaussion in the OGGM Slack (message in support the 15th Jan 2023), "the glacier directories are in a Transverse Mercator map projection centered on the glacier. It's exactly the same as UTM, but with slightly different project parameters for each glacier."

I am trying to recover the exact coordinates based on this logic using Geodesy.jl, but I am not having much success. Here an example with Bossons glacier (RGI60-11.00773). Using the central latitude and longitude obtained from glacier_grid,

using Geodesy

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)
# Returns latitude=46.87433853329805, longitude=10.124544115391744

(the actual coordines for Bossons should be latitude~=45.86, longitude~=6.86)

I would like to include this data in Glacier2D or create a simple function that allows to recover latitude and longitude from it.

@albangossard @JordiBolibar by any chance did you encounter a similar problem when doing this? I don't think @JordiBolibar and I encounter this problem before, since we never really needed the exact latitudes and longitudes at each pixel.

[fmaussion]

Mh, a few things come to mind:

• you don't seem to be using the lon_0 parameter for the projection. In python the projection parameters (lon_0, k) need to be passed when you instanciate the proj
• x0y0 refers to the upper left corner of the map - not the center of the glacier. from this and dxdy you can reconstruct the entire grid coordinates

[fmaussion]

Ups sorry now I see you pass it in the second line - ok yeah then I don't know.

If you want we have UTM grids as well? But I'm not sure if it helps

[facusapienza21]

Following @albangossard suggestion, I swapped to CoordRefSystems.jl, which I find more intuitive. The problem does not seems to be coming from the software side, but rather from the glacier data.

Here an example of a wrong answer obtained for Bossons glacier (RGI RGI60-11.00773"):

using Unitful: m, rad, °, ppm
using CoordRefSystems

k₀ = 0.9996
xₒ = 0.0m
yₒ = 0.0m

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°

The answer here is incorrect. All these parameters are extracted from the gdirs, specially the shift in longitude and k.

Now, if we do the same for 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)

and 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

then the answers are correct.

Is it possible @fmaussion that the problem is coming from the data? I checked a couple of times the RGI ID and I think it is correct... maybe I am missing something very naive...

[JordiBolibar]

I just want to mention that I'm really happy to see that you are using Unitful.jl. It's been on the wishlist for a while. It would be awesome to extend it to the whole ecosystem. I even created an issue for that a while ago :)

[fmaussion]

The RGI ID of Bossons is RGI60-11.03646 - "RGI60-11.00773" is somewhere at the swiss / Austrian border!

[JordiBolibar]

I'm afraid this could be my bad, I used GitHub copilot to generate that list as a quick example and didn't check it in detail...

[facusapienza21]

Jajaja and you still want me to use ChatGPT and its relative cousins, @JordiBolibar ? Just a joke.

I don't think this was actually your fault. I double checked multiple times the RGI and for some reason I am using the GLIMS viewer interface wrong:

I searched for Bossons and this is the answer, reason why I used this RGI ID. Maybe @JordiBolibar you made the same mistake?

[fmaussion]

This is RGI7 - for RGI6 one needs to select the layer

[facusapienza21]

mmm not sure what that means... I guess I should read more carefully the docs XD

[fmaussion]

https://tutorials.oggm.org/stable/notebooks/tutorials/working_with_rgi.html#find-out-the-id-of-a-glacier