Skip to content

Commit

Permalink
Fix tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@fmaussion
fmaussion committed Nov 20, 2024
1 parent 86e67d3 commit 79c95d3
Show file tree
Hide file tree
Showing 2 changed files with 41 additions and 34 deletions.
51 changes: 28 additions & 23 deletions oggm/shop/its_live.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,15 +26,6 @@

regions = [f'RGI{reg:02d}A' for reg in range(1, 20)]

region_files = {}
for reg in regions:
d = {}
for var in ['v', 'vx', 'vy', 'v_error', 'vx_error', 'vy_error']:
d[var] = base_url + f'{reg}_0000_v02_{var}.tif'
region_files[reg] = d

region_grids = {}

rgi_region_links = {'01': 'RGI01A',
'02': ['RGI01A', 'RGI02A'],
'03': 'RGI03A',
Expand All @@ -54,14 +45,24 @@
}


def find_region(gdir):
def _find_region(gdir):
return rgi_region_links.get(gdir.rgi_region, None)


def _region_files():
region_files = {}
for reg in regions:
d = {}
for var in ['v', 'vx', 'vy', 'v_error', 'vx_error', 'vy_error']:
d[var] = base_url + f'{reg}_0000_v02_{var}.tif'
region_files[reg] = d
return region_files


def _reproject_and_scale(gdir, do_error=False):
"""Reproject and scale itslive data, avoid code duplication for error"""

reg = find_region(gdir)
reg = _find_region(gdir)
if reg is None:
raise InvalidWorkflowError('There does not seem to be its_live data '
'available for this glacier')
Expand All @@ -75,9 +76,10 @@ def _reproject_and_scale(gdir, do_error=False):
vny += '_error'

with utils.get_lock():
fv = utils.file_downloader(region_files[reg][vn])
fx = utils.file_downloader(region_files[reg][vnx])
fy = utils.file_downloader(region_files[reg][vny])
reg_files = _region_files()
fv = utils.file_downloader(reg_files[reg][vn])
fx = utils.file_downloader(reg_files[reg][vnx])
fy = utils.file_downloader(reg_files[reg][vny])

# Open the files
dsv = salem.GeoTiff(fv)
Expand Down Expand Up @@ -136,19 +138,22 @@ def _reproject_and_scale(gdir, do_error=False):
with warnings.catch_warnings():
warnings.filterwarnings("ignore", category=RuntimeWarning)
new_vel = np.sqrt(vx**2 + vy**2)
p_ok = new_vel > 0 # avoid div by zero
vx[p_ok] = vx[p_ok] * vel[p_ok] / new_vel[p_ok]
vy[p_ok] = vy[p_ok] * vel[p_ok] / new_vel[p_ok]
p_ok = np.isfinite(new_vel) & (new_vel > 0) # avoid div by zero
scaler = vel[p_ok] / new_vel[p_ok]
vx[p_ok] = vx[p_ok] * scaler
vy[p_ok] = vy[p_ok] * scaler

# And transform to local map
v = grid_gla.map_gridded_data(vel, grid=grid_vel, interp='linear')
vv = grid_gla.map_gridded_data(vel, grid=grid_vel, interp='linear')
vx = grid_gla.map_gridded_data(vx, grid=grid_vel, interp='linear')
vy = grid_gla.map_gridded_data(vy, grid=grid_vel, interp='linear')

# Write
with utils.ncDataset(gdir.get_filepath('gridded_data'), 'a') as nc:

vn = 'itslive_v'
if do_error:
vn += '_error'
if vn in nc.variables:
v = nc.variables[vn]
else:
Expand All @@ -158,11 +163,11 @@ def _reproject_and_scale(gdir, do_error=False):
if do_error:
ln = 'Uncertainty of ' + ln
v.long_name = ln
v[:] = v.filled(np.nan)
v[:] = vv.filled(np.nan)

vn = 'itslive_vx'
if do_error:
vn += '_err'
vn += '_error'
if vn in nc.variables:
v = nc.variables[vn]
else:
Expand All @@ -176,7 +181,7 @@ def _reproject_and_scale(gdir, do_error=False):

vn = 'itslive_vy'
if do_error:
vn += '_err'
vn += '_error'
if vn in nc.variables:
v = nc.variables[vn]
else:
Expand Down Expand Up @@ -252,8 +257,8 @@ def itslive_statistics(gdir):
warnings.filterwarnings('ignore', category=RuntimeWarning)
d['itslive_avg_vel'] = np.nanmean(v)
d['itslive_max_vel'] = np.nanmax(v)
d['itslive_perc_cov'] = (float((~v.isnull()).sum() * gdir.grid.dx ** 2 * 1e-6) /
gridded_area)
d['itslive_perc_cov'] = float(((~v.isnull()).sum() * gdir.grid.dx ** 2 * 1e-6) /
gridded_area)
except (FileNotFoundError, AttributeError, KeyError):
pass

Expand Down
24 changes: 13 additions & 11 deletions oggm/tests/test_shop.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,11 +40,9 @@ def test_repro_to_glacier(self, class_case_dir, monkeypatch):
tasks.glacier_masks(gdir)

# use our files
region_files = {'ALA': {
'vx': get_demo_file('crop_ALA_G0120_0000_vx.tif'),
'vy': get_demo_file('crop_ALA_G0120_0000_vy.tif')}}
monkeypatch.setattr(its_live, 'region_files', region_files)
monkeypatch.setattr(utils, 'file_downloader', lambda x: x)
base_url = ('https://cluster.klima.uni-bremen.de/~oggm/test_files/'
'itslive_v2/ITS_LIVE_velocity_120m_')
monkeypatch.setattr(its_live, 'base_url', base_url)

with warnings.catch_warnings():
warnings.filterwarnings("ignore", category=RuntimeWarning)
Expand All @@ -60,11 +58,14 @@ def test_repro_to_glacier(self, class_case_dir, monkeypatch):
assert np.nanmin(vel) < 2

# We reproject with rasterio and check no big diff
reg_files = its_live._region_files()
cfg.BASENAMES['itslive_vx'] = ('itslive_vx.tif', '')
cfg.BASENAMES['itslive_vy'] = ('itslive_vy.tif', '')
gis.rasterio_to_gdir(gdir, region_files['ALA']['vx'], 'itslive_vx',
gis.rasterio_to_gdir(gdir, reg_files['RGI01A']['vx'],
'itslive_vx',
resampling='bilinear')
gis.rasterio_to_gdir(gdir, region_files['ALA']['vy'], 'itslive_vy',
gis.rasterio_to_gdir(gdir, reg_files['RGI01A']['vy'],
'itslive_vy',
resampling='bilinear')

with rioxr.open_rasterio(gdir.get_filepath('itslive_vx')) as da:
Expand All @@ -73,8 +74,9 @@ def test_repro_to_glacier(self, class_case_dir, monkeypatch):
_vy = da.where(mask).data.squeeze()

_vel = np.sqrt(_vx**2 + _vy**2)
mask = mask & np.isfinite(_vel) & np.isfinite(vel)
np.testing.assert_allclose(utils.rmsd(vel[mask], _vel[mask]), 0,
atol=40)
atol=100)
np.testing.assert_allclose(utils.md(vel[mask], _vel[mask]), 0,
atol=8)

Expand Down Expand Up @@ -138,7 +140,7 @@ def test_thickvel_to_glacier(self, class_case_dir, monkeypatch):
thick = ds.millan_ice_thickness.where(mask).data

# Simply some coverage and sanity checks
assert np.isfinite(thick).sum() / mask.sum() > 0.98
assert np.isfinite(thick).sum() / mask.sum() > 0.92
assert np.nanmax(thick) > 800
assert np.nansum(thick) * gdir.grid.dx**2 * 1e-9 > 174

Expand All @@ -152,7 +154,7 @@ def test_thickvel_to_glacier(self, class_case_dir, monkeypatch):
vy = ds.millan_vy.where(mask).data

# Simply some coverage and sanity checks
assert np.isfinite(v).sum() / mask.sum() > 0.98
assert np.isfinite(v).sum() / mask.sum() > 0.92
assert np.nanmax(v) > 2000
assert np.nanmax(vx) > 500
assert np.nanmax(vy) > 400
Expand All @@ -162,7 +164,7 @@ def test_thickvel_to_glacier(self, class_case_dir, monkeypatch):
assert df['millan_avg_vel'] > 180
assert df['millan_max_vel'] > 2000
assert df['millan_perc_cov'] > 0.95
assert df['millan_vel_perc_cov'] > 0.97
assert df['millan_vel_perc_cov'] > 0.92
assert df['millan_vol_km3'] > 174


Expand Down

0 comments on commit 79c95d3

Please sign in to comment.