Skip to content

small bug in climate selection for glaciers near 0 longitude #1825

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

Open
lilianschuster wants to merge 2 commits into
base: master
Choose a base branch
from
+57 −15
Open

small bug in climate selection for glaciers near 0 longitude #1825

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
f63705e
found small bug in climate selection for glaciers near 0 longitude
lilianschuster Dec 8, 2025
4dc424b
updated tests to work w. gridpoint near zero
lilianschuster Dec 8, 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
44 changes: 37 additions & 7 deletions oggm/shop/ecmwf.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -134,7 +134,6 @@ def process_ecmwf_data(gdir, dataset=None, ensemble_member=0,
if dataset is None:
dataset = cfg.PARAMS['baseline_climate']

# Use xarray to read the data
lon = gdir.cenlon + 360 if gdir.cenlon < 0 else gdir.cenlon
lat = gdir.cenlat
with xr.open_dataset(get_ecmwf_file(dataset, 'tmp')) as ds:
Expand All @@ -147,13 +146,20 @@ def process_ecmwf_data(gdir, dataset=None, ensemble_member=0,
assert ds['time.month'][-1] == 5
y1 -= 1
ds = ds.sel(time=slice(f'{y0}-01-01', f'{y1}-12-01'))
lon_sp = lon
if np.any(ds.longitude==0):
res = float(360-ds.longitude[-1]) # ("flattened") grid resolution
wrap_limit = 360 - res / 2
if lon > wrap_limit:
lon_sp = lon - 360
if dataset == 'CERA':
ds = ds.sel(number=ensemble_member)
try:
ds = ds.sel(longitude=lon, latitude=lat, method='nearest')
ds = ds.sel(longitude=lon_sp, latitude=lat, method='nearest')
except (ValueError, KeyError):
# Flattened ERA5
c = (ds.longitude - lon)**2 + (ds.latitude - lat)**2
# there is an ERA5 gridpoint at longitude 0 (i.e., equal to 360)
c = (ds.longitude - lon_sp) ** 2 + (ds.latitude - lat) ** 2
ds = ds.isel(points=np.argmin(c.data))
temp = ds['t2m'].data - 273.15
time = ds.time.data
Expand All @@ -163,28 +169,40 @@ def process_ecmwf_data(gdir, dataset=None, ensemble_member=0,
with xr.open_dataset(get_ecmwf_file(dataset, 'pre')) as ds:
_check_ds_validity(ds)
ds = ds.sel(time=slice(f'{y0}-01-01', f'{y1}-12-01'))
lon_sp = lon
if np.any(ds.longitude==0):
res = float(360-ds.longitude[-1]) # ("flattened") grid resolution
wrap_limit = 360 - res / 2
if lon > wrap_limit:
lon_sp = lon - 360
if dataset == 'CERA':
ds = ds.sel(number=ensemble_member)
try:
ds = ds.sel(longitude=lon, latitude=lat, method='nearest')
ds = ds.sel(longitude=lon_sp, latitude=lat, method='nearest')
except (ValueError, KeyError):
# Flattened ERA5
c = (ds.longitude - lon)**2 + (ds.latitude - lat)**2
c = (ds.longitude - lon_sp) ** 2 + (ds.latitude - lat) ** 2
ds = ds.isel(points=np.argmin(c.data))
prcp = ds['tp'].data * 1000 * ds['time.daysinmonth']
with xr.open_dataset(get_ecmwf_file(dataset, 'inv')) as ds:
_check_ds_validity(ds)
lon_sp = lon
if np.any(ds.longitude==0):
res = float(360-ds.longitude[-1]) # ("flattened") grid resolution
wrap_limit = 360 - res / 2
if lon > wrap_limit:
lon_sp = lon - 360
try:
ds = ds.isel(time=0)
except (ValueError):
# new inv flattening files do not have any
# time dependencies anymore
pass
try:
ds = ds.sel(longitude=lon, latitude=lat, method='nearest')
ds = ds.sel(longitude=lon_sp, latitude=lat, method='nearest')
except (ValueError, KeyError):
# Flattened ERA5
c = (ds.longitude - lon)**2 + (ds.latitude - lat)**2
c = (ds.longitude - lon_sp) ** 2 + (ds.latitude - lat) ** 2
ds = ds.isel(points=np.argmin(c.data))
hgt = ds['z'].data / cfg.G

Expand All @@ -194,10 +212,22 @@ def process_ecmwf_data(gdir, dataset=None, ensemble_member=0,
with xr.open_dataset(get_ecmwf_file(dataset, 'lapserates')) as ds:
_check_ds_validity(ds)
ds = ds.sel(time=slice(f'{y0}-01-01', f'{y1}-12-01'))
lon_sp = lon
if np.any(ds.longitude == 0):
res = float(360 - ds.longitude[-1]) # ("flattened") grid resolution
wrap_limit = 360 - res / 2
if lon > wrap_limit:
lon_sp = lon - 360
ds = ds.sel(longitude=lon, latitude=lat, method='nearest')

with xr.open_dataset(get_ecmwf_file(dataset, 'tempstd')) as ds:
_check_ds_validity(ds)
lon_sp = lon
if np.any(ds.longitude == 0):
res = float(360 - ds.longitude[-1]) # ("flattened") grid resolution
wrap_limit = 360 - res / 2
if lon > wrap_limit:
lon_sp = lon - 360
ds = ds.sel(time=slice(f'{y0}-01-01', f'{y1}-12-01'))
ds = ds.sel(longitude=lon, latitude=lat, method='nearest')
temp_std = ds['t2m_std'].data
Expand Down
28 changes: 20 additions & 8 deletions oggm/tests/test_shop.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -313,7 +313,7 @@ def test_glacier_gridpoint_selection(self):
# https://nbviewer.org/urls/cluster.klima.uni-bremen.de/
# ~lschuster/example_ipynb/flatten_glacier_gridpoint_tests.ipynb
if d == 'GSWP3_W5E5':
res = 0.5
res = 0.5 # 0.5° resolution
with xr.open_dataset(w5e5.get_gswp3_w5e5_file(d, 'inv')) as dinv:
dinv = dinv.load()
else:
Expand All @@ -328,11 +328,24 @@ def test_glacier_gridpoint_selection(self):
(51.495, 30.9010), # RGI60-12.01691
(0, 0), # random gridpoint not near to a glacier
( -141.670274+360, 69.166921), # in RGI7C, not in RGI6
(-66.855668+360, -67.535551) # only in RGI7G, not in RGI6 or in RGI 7C
(-66.855668+360, -67.535551), # only in RGI7G, not in RGI6 or in RGI 7C
(-0.039683+360, 42.695419 ), # RGI60-11.03228 also check glacier near longitude 0
(-179.915527+360, 66.276108) # RGI60-10.05049 (near -180 longitude)
]:
lon, lat = coord
# get the distances to the glacier coordinate
c = (dinv.longitude - lon) ** 2 + (dinv.latitude - lat) ** 2

#this does not work if it is flattened...
# res = float(dinv.longitude[1] - dinv.longitude[0]) # grid resolution
wrap_limit = 360 - res / 2
lon_sp = lon
if (lon > wrap_limit) and np.any(dinv.longitude==0):
# there is an ERA5 gridpoint at longitude 0 (i.e., equal to 360)
# we want to actually select longitude 0 if we are near to 360, so we reconvert it
lon_sp = lon - 360
# W5E5 gridpoints are only at 0.25°, 0.75°..., not exactly at 0 and do not encounter that issue
assert d != 'GSWP3-W5E5'
c = (dinv.longitude - lon_sp) ** 2 + (dinv.latitude - lat) ** 2
c = c.to_dataframe('distance').sort_values('distance')
# select the nearest climate point from the flattened glacier gridpoint
lat_near, lon_near, dist = c.iloc[0]
Expand All @@ -343,14 +356,13 @@ def test_glacier_gridpoint_selection(self):
assert np.abs(lon_near - lon) <= res/2
assert dist <= ((res/2) ** 2 + (res/2)** 2) ** 0.5
# for glaciers the next gridpoint should be the nearest
# (GSWP3-W5E5 resolution is 0.5°)
else:
assert np.abs(lat_near - lat) <= res/2
assert np.abs(lon_near - lon) <= res/2
assert np.abs(lon_near - lon_sp) <= res/2
assert dist <= ((res/2) ** 2 + (res/2) ** 2) ** 0.5

# this only contains data for five glaciers, let's still check some basics
# both glaciers are not at latitude or longitude 0
# this only contains data for a few glaciers, let's still check some basics
# both glaciers are not at latitude or longitude 0, but for one glacier we need to pick longitude 0
if d == 'GSWP3_W5E5':
with xr.open_dataset(w5e5.get_gswp3_w5e5_file(d, 'temp_std')) as dtemp_std:
assert np.all(dtemp_std.latitude != 0)
Expand All @@ -361,7 +373,7 @@ def test_glacier_gridpoint_selection(self):
else:
with xr.open_dataset(ecmwf.get_ecmwf_file('ERA5', 'tmp')) as dtemp:
assert np.all(dtemp.latitude != 0)
assert np.all(dtemp.longitude != 0)
assert np.all(dtemp.longitude >= 0)
assert dtemp.isel(time=0).t2m.std() > 0
assert dtemp.longitude.std() > 0
assert dtemp.latitude.std() > 0
Expand Down