Skip to content

Commit

Permalink
fix cesm warnings in tests (#314)
Browse files Browse the repository at this point in the history 
* fix cesm warnings in tests

* Add some ELA tests

* One more test
  • Loading branch information
@fmaussion
fmaussion committed Oct 6, 2017
1 parent 60f5f5c commit 38ecf1c
Show file tree
Hide file tree
Showing 2 changed files with 109 additions and 78 deletions.
93 changes: 56 additions & 37 deletions oggm/tests/test_models.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1644,39 +1644,50 @@ def test_cesm(self):
gdir = self.gdir

# init
cfg.PATHS['gcm_temp_file'] = get_demo_file('cesm.TREFHT.160001-200512.selection.nc')
cfg.PATHS['gcm_precc_file'] = get_demo_file('cesm.PRECC.160001-200512.selection.nc')
cfg.PATHS['gcm_precl_file'] = get_demo_file('cesm.PRECL.160001-200512.selection.nc')
f = get_demo_file('cesm.TREFHT.160001-200512.selection.nc')
cfg.PATHS['gcm_temp_file'] = f
f = get_demo_file('cesm.PRECC.160001-200512.selection.nc')
cfg.PATHS['gcm_precc_file'] = f
f = get_demo_file('cesm.PRECL.160001-200512.selection.nc')
cfg.PATHS['gcm_precl_file'] = f
climate.process_cesm_data(self.gdir)

# CLimate data
with xr.open_dataset(gdir.get_filepath('climate_monthly')) as hist, \
xr.open_dataset(gdir.get_filepath('cesm_data')) as cesm:

time = pd.period_range(cesm.time.values[0].strftime('%Y-%m-%d'),
cesm.time.values[-1].strftime('%Y-%m-%d'),
freq='M')
cesm['time'] = time
cesm.coords['year'] = ('time', time.year)
cesm.coords['month'] = ('time', time.month)

# Let's do some basic checks
shist = hist.sel(time=slice('1961', '1990'))
scesm = cesm.isel(time=(cesm.year >= 1961) & (cesm.year <= 1990))
# Climate during the chosen period should be the same
np.testing.assert_allclose(shist.temp.mean(),
scesm.temp.mean(),
rtol=1e-3)
np.testing.assert_allclose(shist.prcp.mean(),
scesm.prcp.mean(),
rtol=1e-3)
np.testing.assert_allclose(shist.grad.mean(), scesm.grad.mean())
# And also the anual cycle
scru = shist.groupby('time.month').mean()
scesm = scesm.groupby(scesm.month).mean()
np.testing.assert_allclose(scru.temp, scesm.temp, rtol=5e-3)
np.testing.assert_allclose(scru.prcp, scesm.prcp, rtol=1e-3)
np.testing.assert_allclose(scru.grad, scesm.grad)
# Climate data
with warnings.catch_warnings():
# Long time series are currently a pain pandas
warnings.filterwarnings("ignore",
message='Unable to decode time axis')
fh = gdir.get_filepath('climate_monthly')
fcesm = gdir.get_filepath('cesm_data')
with xr.open_dataset(fh) as hist, xr.open_dataset(fcesm) as cesm:

tv = cesm.time.values
time = pd.period_range(tv[0].strftime('%Y-%m-%d'),
tv[-1].strftime('%Y-%m-%d'),
freq='M')
cesm['time'] = time
cesm.coords['year'] = ('time', time.year)
cesm.coords['month'] = ('time', time.month)

# Let's do some basic checks
shist = hist.sel(time=slice('1961', '1990'))
scesm = cesm.isel(time=((cesm.year >= 1961) &
(cesm.year <= 1990)))
# Climate during the chosen period should be the same
np.testing.assert_allclose(shist.temp.mean(),
scesm.temp.mean(),
rtol=1e-3)
np.testing.assert_allclose(shist.prcp.mean(),
scesm.prcp.mean(),
rtol=1e-3)
np.testing.assert_allclose(shist.grad.mean(),
scesm.grad.mean())
# And also the anual cycle
scru = shist.groupby('time.month').mean()
scesm = scesm.groupby(scesm.month).mean()
np.testing.assert_allclose(scru.temp, scesm.temp, rtol=5e-3)
np.testing.assert_allclose(scru.prcp, scesm.prcp, rtol=1e-3)
np.testing.assert_allclose(scru.grad, scesm.grad)

# Mass balance models
mb_cru = massbalance.PastMassBalanceModel(self.gdir)
Expand All @@ -1698,23 +1709,29 @@ def test_cesm(self):
ts2 = mb_cesm.get_specific_mb(h, w, year=yrs)
if do_plot:
df = pd.DataFrame(index=yrs)
k1 = 'Histalp (mean={:.1f}, stddev={:.1f})'.format(np.mean(ts1), np.std(ts1))
k2 = 'CESM (mean={:.1f}, stddev={:.1f})'.format(np.mean(ts2), np.std(ts2))
k1 = 'Histalp (mean={:.1f}, stddev={:.1f})'.format(np.mean(ts1),
np.std(ts1))
k2 = 'CESM (mean={:.1f}, stddev={:.1f})'.format(np.mean(ts2),
np.std(ts2))
df[k1] = ts1
df[k2] = ts2

df.plot()
plt.plot(yrs, df[k1].rolling(31, center=True, min_periods=15).mean(),
plt.plot(yrs,
df[k1].rolling(31, center=True, min_periods=15).mean(),
color='C0', linewidth=3)
plt.plot(yrs, df[k2].rolling(31, center=True, min_periods=15).mean(),
plt.plot(yrs,
df[k2].rolling(31, center=True, min_periods=15).mean(),
color='C1', linewidth=3)
plt.title('SMB Hintereisferner Histalp VS CESM')
plt.show()

# See what that means for a run
flowline.init_present_time_glacier(gdir)
flowline.run_from_climate_data(gdir, filesuffix='_hist')
flowline.run_from_climate_data(gdir, filename='cesm_data',
flowline.run_from_climate_data(gdir, ys=1961, ye=1990,
filesuffix='_hist')
flowline.run_from_climate_data(gdir, ys=1961, ye=1990,
filename='cesm_data',
filesuffix='_cesm')

ds1 = utils.compile_run_output([gdir], path=False, filesuffix='_hist')
Expand All @@ -1723,6 +1740,8 @@ def test_cesm(self):
assert_allclose(ds1.volume.isel(rgi_id=0, time=-1),
ds2.volume.isel(rgi_id=0, time=-1),
rtol=0.1)
# ELA should be close
assert_allclose(ds1.ela.mean(), ds2.ela.mean(), atol=50)

@is_slow
def test_elevation_feedback(self):
Expand Down
94 changes: 53 additions & 41 deletions oggm/tests/test_prepro.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1970,48 +1970,60 @@ def test_process_cesm(self):
self.assertEqual(ci['hydro_yr_0'], 1902)
self.assertEqual(ci['hydro_yr_1'], 2014)

cfg.PATHS['gcm_temp_file'] = get_demo_file('cesm.TREFHT.160001-200512.selection.nc')
cfg.PATHS['gcm_precc_file'] = get_demo_file('cesm.PRECC.160001-200512.selection.nc')
cfg.PATHS['gcm_precl_file'] = get_demo_file('cesm.PRECL.160001-200512.selection.nc')
f = get_demo_file('cesm.TREFHT.160001-200512.selection.nc')
cfg.PATHS['gcm_temp_file'] = f
f = get_demo_file('cesm.PRECC.160001-200512.selection.nc')
cfg.PATHS['gcm_precc_file'] = f
f = get_demo_file('cesm.PRECL.160001-200512.selection.nc')
cfg.PATHS['gcm_precl_file'] = f
climate.process_cesm_data(gdir)

with xr.open_dataset(gdir.get_filepath('climate_monthly')) as cru, \
xr.open_dataset(gdir.get_filepath('cesm_data')) as cesm:

time = pd.period_range(cesm.time.values[0].strftime('%Y-%m-%d'),
cesm.time.values[-1].strftime('%Y-%m-%d'),
freq='M')
cesm['time'] = time
cesm.coords['year'] = ('time', time.year)
cesm.coords['month'] = ('time', time.month)

# Let's do some basic checks
scru = cru.sel(time=slice('1961', '1990'))
scesm = cesm.isel(time=(cesm.year >= 1961) & (cesm.year <= 1990))
# Climate during the chosen period should be the same
np.testing.assert_allclose(scru.temp.mean(),
scesm.temp.mean(),
rtol=1e-3)
np.testing.assert_allclose(scru.prcp.mean(),
scesm.prcp.mean(),
rtol=1e-3)
np.testing.assert_allclose(scru.grad.mean(), scesm.grad.mean())
# And also the anual cycle
scru = scru.groupby('time.month').mean()
scesm = scesm.groupby(scesm.month).mean()
np.testing.assert_allclose(scru.temp, scesm.temp, rtol=1e-3)
np.testing.assert_allclose(scru.prcp, scesm.prcp, rtol=1e-3)
np.testing.assert_allclose(scru.grad, scesm.grad)

# How did the annua cycle change with time?
scesm1 = cesm.isel(time=(cesm.year >= 1961) & (cesm.year <= 1990))
scesm2 = cesm.isel(time=(cesm.year >= 1661) & (cesm.year <= 1690))
scesm1 = scesm1.groupby(scesm1.month).mean()
scesm2 = scesm2.groupby(scesm2.month).mean()
# No more than one degree? (silly test)
np.testing.assert_allclose(scesm1.temp, scesm2.temp, atol=1)
# N more than 20%? (silly test)
np.testing.assert_allclose(scesm1.prcp, scesm2.prcp, rtol=0.2)
with warnings.catch_warnings():
# Long time series are currently a pain pandas
warnings.filterwarnings("ignore",
message='Unable to decode time axis')
fh = gdir.get_filepath('climate_monthly')
fcesm = gdir.get_filepath('cesm_data')
with xr.open_dataset(fh) as cru, xr.open_dataset(fcesm) as cesm:

tv = cesm.time.values
time = pd.period_range(tv[0].strftime('%Y-%m-%d'),
tv[-1].strftime('%Y-%m-%d'),
freq='M')
cesm['time'] = time
cesm.coords['year'] = ('time', time.year)
cesm.coords['month'] = ('time', time.month)

# Let's do some basic checks
scru = cru.sel(time=slice('1961', '1990'))
scesm = cesm.isel(time=((cesm.year >= 1961) &
(cesm.year <= 1990)))
# Climate during the chosen period should be the same
np.testing.assert_allclose(scru.temp.mean(),
scesm.temp.mean(),
rtol=1e-3)
np.testing.assert_allclose(scru.prcp.mean(),
scesm.prcp.mean(),
rtol=1e-3)
np.testing.assert_allclose(scru.grad.mean(),
scesm.grad.mean())
# And also the anual cycle
scru = scru.groupby('time.month').mean()
scesm = scesm.groupby(scesm.month).mean()
np.testing.assert_allclose(scru.temp, scesm.temp, rtol=1e-3)
np.testing.assert_allclose(scru.prcp, scesm.prcp, rtol=1e-3)
np.testing.assert_allclose(scru.grad, scesm.grad)

# How did the annua cycle change with time?
scesm1 = cesm.isel(time=((cesm.year >= 1961) &
(cesm.year <= 1990)))
scesm2 = cesm.isel(time=((cesm.year >= 1661) &
(cesm.year <= 1690)))
scesm1 = scesm1.groupby(scesm1.month).mean()
scesm2 = scesm2.groupby(scesm2.month).mean()
# No more than one degree? (silly test)
np.testing.assert_allclose(scesm1.temp, scesm2.temp, atol=1)
# N more than 20%? (silly test)
np.testing.assert_allclose(scesm1.prcp, scesm2.prcp, rtol=0.2)


class TestCatching(unittest.TestCase):
Expand Down

0 comments on commit 38ecf1c

Please sign in to comment.