diff --git a/oggm/shop/its_live.py b/oggm/shop/its_live.py index 5cb473142..df7d8e147 100644 --- a/oggm/shop/its_live.py +++ b/oggm/shop/its_live.py @@ -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', @@ -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') @@ -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) @@ -136,12 +138,13 @@ 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') @@ -149,6 +152,8 @@ def _reproject_and_scale(gdir, do_error=False): 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: @@ -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: @@ -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: @@ -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 diff --git a/oggm/tests/test_shop.py b/oggm/tests/test_shop.py index 83a84e2dd..061aad7d4 100644 --- a/oggm/tests/test_shop.py +++ b/oggm/tests/test_shop.py @@ -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) @@ -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: @@ -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) @@ -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 @@ -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 @@ -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