Skip to content
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

Jump to bottom

Update velocity and bedmachine in shop #1753

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Update velocity and bedmachine in shop #1753

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
86e67d3
Update velocity and bedmachine in shop
fmaussion Nov 20, 2024
79c95d3
Fix tests
fmaussion Nov 20, 2024
d2b53d7
add bedmachine to prepro
fmaussion Nov 21, 2024
28139ae
fix integer
fmaussion Nov 21, 2024
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
Prev Previous commit
Next Next commit
Fix tests
  • Loading branch information
@fmaussion
fmaussion committed Nov 20, 2024
commit 79c95d3522080d06fe20e49617280a5c0ba0887d
51 changes: 28 additions & 23 deletions oggm/shop/its_live.py
View file
Open in desktop
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
Open in desktop
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