Add support for LMR data to shop (#1257)

docs/whats-new.rst

@@ -16,6 +16,8 @@ Enhancements
 - Added ``cook_rgidf()`` function in ``oggm.utils`` to simplify the use
   of a non-RGI glacier inventory in OGGM (:pull:`1251`).
   By `Li Fei <https://github.com/Keeptg>`_
+- Added support for last millenium reanalysis data to the shop (:pull:`1257`).
+  By `Fabien Maussion <https://github.com/fmaussion>`_
 
 Bug fixes
 ~~~~~~~~~

oggm/shop/gcm_climate.py

@@ -5,6 +5,7 @@
 import warnings
 
 # External libs
+import cftime
 import numpy as np
 import netCDF4
 import xarray as xr
@@ -86,79 +87,76 @@ def process_gcm_data(gdir, filesuffix='', prcp=None, temp=None,
     assert len(prcp) // 12 == len(prcp) / 12, 'Somehow we didn\'t get full years'
     assert len(temp) // 12 == len(temp) / 12, 'Somehow we didn\'t get full years'
 
-    # Get CRU to apply the anomaly to
+    # Get the reference data to apply the anomaly to
     fpath = gdir.get_filepath('climate_historical')
-    ds_cru = xr.open_dataset(fpath)
+    with xr.open_dataset(fpath) as ds_ref:
+
+        ds_ref = ds_ref.sel(time=slice(*year_range))
+
+        # compute monthly anomalies
+        # of temp
+        if scale_stddev:
+            # This is a bit more arithmetic
+            ts_tmp_sel = temp.sel(time=slice(*year_range))
+            ts_tmp_std = ts_tmp_sel.groupby('time.month').std(dim='time')
+            std_fac = ds_ref.temp.groupby('time.month').std(dim='time') / ts_tmp_std
+            std_fac = std_fac.roll(month=13-sm, roll_coords=True)
+            std_fac = np.tile(std_fac.data, len(temp) // 12)
+            # We need an even number of years for this to work
+            if ((len(ts_tmp_sel) // 12) % 2) == 1:
+                raise InvalidParamsError('We need an even number of years '
+                                         'for this to work')
+            win_size = len(ts_tmp_sel) + 1
+
+            def roll_func(x, axis=None):
+                x = x[:, ::12]
+                n = len(x[0, :]) // 2
+                xm = np.nanmean(x, axis=axis)
+                return xm + (x[:, n] - xm) * std_fac
+
+            temp = temp.rolling(time=win_size, center=True,
+                                min_periods=1).reduce(roll_func)
 
-    # Add CRU clim
-    dscru = ds_cru.sel(time=slice(*year_range))
-
-    # compute monthly anomalies
-    # of temp
-    if scale_stddev:
-        # This is a bit more arithmetic
         ts_tmp_sel = temp.sel(time=slice(*year_range))
-        ts_tmp_std = ts_tmp_sel.groupby('time.month').std(dim='time')
-        std_fac = dscru.temp.groupby('time.month').std(dim='time') / ts_tmp_std
-        std_fac = std_fac.roll(month=13-sm, roll_coords=True)
-        std_fac = np.tile(std_fac.data, len(temp) // 12)
-        # We need an even number of years for this to work
-        if ((len(ts_tmp_sel) // 12) % 2) == 1:
-            raise InvalidParamsError('We need an even number of years '
-                                     'for this to work')
-        win_size = len(ts_tmp_sel) + 1
-
-        def roll_func(x, axis=None):
-            x = x[:, ::12]
-            n = len(x[0, :]) // 2
-            xm = np.nanmean(x, axis=axis)
-            return xm + (x[:, n] - xm) * std_fac
-
-        temp = temp.rolling(time=win_size, center=True,
-                            min_periods=1).reduce(roll_func)
-
-    ts_tmp_sel = temp.sel(time=slice(*year_range))
-    ts_tmp_avg = ts_tmp_sel.groupby('time.month').mean(dim='time')
-    ts_tmp = temp.groupby('time.month') - ts_tmp_avg
-    # of precip -- scaled anomalies
-    ts_pre_avg = prcp.sel(time=slice(*year_range))
-    ts_pre_avg = ts_pre_avg.groupby('time.month').mean(dim='time')
-    ts_pre_ano = prcp.groupby('time.month') - ts_pre_avg
-    # scaled anomalies is the default. Standard anomalies above
-    # are used later for where ts_pre_avg == 0
-    ts_pre = prcp.groupby('time.month') / ts_pre_avg
-
-    # for temp
-    loc_tmp = dscru.temp.groupby('time.month').mean()
-    ts_tmp = ts_tmp.groupby('time.month') + loc_tmp
-
-    # for prcp
-    loc_pre = dscru.prcp.groupby('time.month').mean()
-    # scaled anomalies
-    ts_pre = ts_pre.groupby('time.month') * loc_pre
-    # standard anomalies
-    ts_pre_ano = ts_pre_ano.groupby('time.month') + loc_pre
-    # Correct infinite values with standard anomalies
-    ts_pre.values = np.where(np.isfinite(ts_pre.values),
-                             ts_pre.values,
-                             ts_pre_ano.values)
-    # The previous step might create negative values (unlikely). Clip them
-    ts_pre.values = utils.clip_min(ts_pre.values, 0)
-
-    assert np.all(np.isfinite(ts_pre.values))
-    assert np.all(np.isfinite(ts_tmp.values))
-
-    gdir.write_monthly_climate_file(temp.time.values,
-                                    ts_pre.values, ts_tmp.values,
-                                    float(dscru.ref_hgt),
-                                    prcp.lon.values, prcp.lat.values,
-                                    time_unit=time_unit,
-                                    calendar=calendar,
-                                    file_name='gcm_data',
-                                    source=source,
-                                    filesuffix=filesuffix)
-
-    ds_cru.close()
+        ts_tmp_avg = ts_tmp_sel.groupby('time.month').mean(dim='time')
+        ts_tmp = temp.groupby('time.month') - ts_tmp_avg
+        # of precip -- scaled anomalies
+        ts_pre_avg = prcp.sel(time=slice(*year_range))
+        ts_pre_avg = ts_pre_avg.groupby('time.month').mean(dim='time')
+        ts_pre_ano = prcp.groupby('time.month') - ts_pre_avg
+        # scaled anomalies is the default. Standard anomalies above
+        # are used later for where ts_pre_avg == 0
+        ts_pre = prcp.groupby('time.month') / ts_pre_avg
+
+        # for temp
+        loc_tmp = ds_ref.temp.groupby('time.month').mean()
+        ts_tmp = ts_tmp.groupby('time.month') + loc_tmp
+
+        # for prcp
+        loc_pre = ds_ref.prcp.groupby('time.month').mean()
+        # scaled anomalies
+        ts_pre = ts_pre.groupby('time.month') * loc_pre
+        # standard anomalies
+        ts_pre_ano = ts_pre_ano.groupby('time.month') + loc_pre
+        # Correct infinite values with standard anomalies
+        ts_pre.values = np.where(np.isfinite(ts_pre.values),
+                                 ts_pre.values,
+                                 ts_pre_ano.values)
+        # The previous step might create negative values (unlikely). Clip them
+        ts_pre.values = utils.clip_min(ts_pre.values, 0)
+
+        assert np.all(np.isfinite(ts_pre.values))
+        assert np.all(np.isfinite(ts_tmp.values))
+
+        gdir.write_monthly_climate_file(temp.time.values,
+                                        ts_pre.values, ts_tmp.values,
+                                        float(ds_ref.ref_hgt),
+                                        prcp.lon.values, prcp.lat.values,
+                                        time_unit=time_unit,
+                                        calendar=calendar,
+                                        file_name='gcm_data',
+                                        source=source,
+                                        filesuffix=filesuffix)
 
 
 @entity_task(log, writes=['gcm_data'])
@@ -332,3 +330,97 @@ def process_cmip_data(gdir, filesuffix='', fpath_temp=None,
 
     process_gcm_data(gdir, filesuffix=filesuffix, prcp=precip, temp=temp,
                      source=filesuffix, **kwargs)
+
+
+@entity_task(log, writes=['gcm_data'])
+def process_lmr_data(gdir, fpath_temp=None, fpath_precip=None,
+                     year_range=('1951', '1980'), filesuffix='', **kwargs):
+    """Read, process and store the Last Millennium Reanalysis (LMR) data for this glacier.
+
+    LMR data: https://atmos.washington.edu/~hakim/lmr/LMRv2/
+
+    LMR data is annualised in anomaly format relative to 1951-1980. We
+    create synthetic timeseries from the reference data.
+
+    It stores the data in a format that can be used by the OGGM mass balance
+    model and in the glacier directory.
+
+    Parameters
+    ----------
+    fpath_temp : str
+        path to the temp file (default: LMR v2.1 from server above)
+    fpath_precip : str
+        path to the precip file (default: LMR v2.1 from server above)
+    year_range : tuple of str
+        the year range for which you want to compute the anomalies. Default
+        for LMR is `('1951', '1980')`
+    filesuffix : str
+        append a suffix to the filename (useful for ensemble experiments).
+
+    **kwargs: any kwarg to be passed to ref:`process_gcm_data`
+    """
+
+    # Get the path of GCM temperature & precipitation data
+    base_url = 'https://atmos.washington.edu/%7Ehakim/lmr/LMRv2/'
+    if fpath_temp is None:
+        fpath_temp = utils.file_downloader(base_url + 'air_MCruns_ensemble_mean_LMRv2.1.nc')
+    if fpath_precip is None:
+        fpath_precip = utils.file_downloader(base_url + 'prate_MCruns_ensemble_mean_LMRv2.1.nc')
+
+    # Glacier location
+    glon = gdir.cenlon
+    glat = gdir.cenlat
+
+    # Read the GCM files
+    with xr.open_dataset(fpath_temp, use_cftime=True) as tempds, \
+            xr.open_dataset(fpath_precip, use_cftime=True) as precipds:
+
+        # Check longitude conventions
+        if tempds.lon.min() >= 0 and glon <= 0:
+            glon += 360
+
+        # Take the closest to the glacier
+        # Should we consider GCM interpolation?
+        temp = tempds.air.sel(lat=glat, lon=glon, method='nearest')
+        precip = precipds.prate.sel(lat=glat, lon=glon, method='nearest')
+
+        # Currently we just take the mean of the ensemble, although
+        # this is probably not advised. The GCM climate will correct
+        # anyways
+        temp = temp.mean(dim='MCrun')
+        precip = precip.mean(dim='MCrun')
+
+        # Precip unit is kg/m^2/s we convert to mm month since we apply the anomaly after
+        precip = precip * 30.5 * (60 * 60 * 24)
+
+        # Back to [-180, 180] for OGGM
+        temp.lon.values = temp.lon if temp.lon <= 180 else temp.lon - 360
+        precip.lon.values = precip.lon if precip.lon <= 180 else precip.lon - 360
+
+    # OK now we have to turn these annual timeseries in monthly data
+    # We take the ref climate
+    fpath = gdir.get_filepath('climate_historical')
+    with xr.open_dataset(fpath) as ds_ref:
+        ds_ref = ds_ref.sel(time=slice(*year_range))
+
+        loc_tmp = ds_ref.temp.groupby('time.month').mean()
+        loc_pre = ds_ref.prcp.groupby('time.month').mean()
+
+        # Make time coord
+        t = np.cumsum([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] * len(temp))
+        t = cftime.num2date(np.append([0], t[:-1]), 'days since 0000-01-01 00:00:00',
+                            calendar='noleap')
+
+        temp = xr.DataArray((loc_tmp.data + temp.data[:, np.newaxis]).flatten(),
+                            coords={'time': t, 'lon': temp.lon, 'lat': temp.lat},
+                            dims=('time',))
+
+        # For precip the std dev is very small - lets keep it as is for now but
+        # this is a bit ridiculous. We clip to zero here to be sure
+        precip = utils.clip_min((loc_pre.data + precip.data[:, np.newaxis]).flatten(), 0)
+        precip = xr.DataArray(precip, dims=('time',),
+                              coords={'time': t, 'lon': temp.lon, 'lat': temp.lat})
+
+    process_gcm_data(gdir, filesuffix=filesuffix, prcp=precip, temp=temp,
+                     year_range=year_range, calendar='noleap',
+                     source='lmr', **kwargs)

oggm/tests/test_prepro.py

@@ -3427,7 +3427,8 @@ def test_process_cmip5(self):
                                        scesm.prcp.mean(),
                                        rtol=1e-3)
 
-            # Here no std dev!
+            # Here also std dev! But its not perfect because std_dev
+            # is preserved over 31 years
             _scru = scru.groupby('time.month').std(dim='time')
             _scesm = scesm.groupby('time.month').std(dim='time')
             assert np.allclose(_scru.temp, _scesm.temp, rtol=1e-2)
@@ -3450,6 +3451,59 @@ def test_process_cmip5(self):
             # N more than 30%? (silly test)
             np.testing.assert_allclose(scmip1.prcp, scmip2.prcp, rtol=0.3)
 
+    def test_process_lmr(self):
+
+        hef_file = get_demo_file('Hintereisferner_RGI5.shp')
+        entity = gpd.read_file(hef_file).iloc[0]
+
+        gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir)
+        gis.define_glacier_region(gdir)
+        tasks.process_cru_data(gdir)
+
+        ci = gdir.get_climate_info()
+        self.assertEqual(ci['baseline_hydro_yr_0'], 1902)
+        self.assertEqual(ci['baseline_hydro_yr_1'], 2014)
+
+        fpath_temp = get_demo_file('air_MCruns_ensemble_mean_LMRv2.1.nc')
+        fpath_precip = get_demo_file('prate_MCruns_ensemble_mean_LMRv2.1.nc')
+        gcm_climate.process_lmr_data(gdir, fpath_temp=fpath_temp,
+                                     fpath_precip=fpath_precip)
+
+        fh = gdir.get_filepath('climate_historical')
+        fcmip = gdir.get_filepath('gcm_data')
+        with xr.open_dataset(fh) as cru, xr.open_dataset(fcmip) as cmip:
+
+            # Let's do some basic checks
+            scru = cru.sel(time=slice('1951', '1980'))
+            scesm = cmip.sel(time=slice('1951', '1980'))
+            # 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)
+
+            # Here also std dev! But its not perfect because std_dev
+            # is preserved over 31 years
+            _scru = scru.groupby('time.month').std(dim='time')
+            _scesm = scesm.groupby('time.month').std(dim='time')
+            np.testing.assert_allclose(_scru.temp, _scesm.temp, rtol=0.15)
+
+            # And also the annual cycle
+            scru = scru.groupby('time.month').mean(dim='time')
+            scesm = scesm.groupby('time.month').mean(dim='time')
+            np.testing.assert_allclose(scru.temp, scesm.temp, rtol=1e-3)
+            np.testing.assert_allclose(scru.prcp, scesm.prcp, rtol=1e-3)
+
+            # How did the annual cycle change with time?
+            scmip1 = cmip.sel(time=slice('1970', '1999'))
+            scmip2 = cmip.sel(time=slice('1800', '1829'))
+            scmip1 = scmip1.groupby('time.month').mean(dim='time')
+            scmip2 = scmip2.groupby('time.month').mean(dim='time')
+            # It has warmed
+            assert scmip2.temp.mean() < scmip1.temp.mean()
+
     def test_process_cmip5_scale(self):
 
         hef_file = get_demo_file('Hintereisferner_RGI5.shp')

oggm/utils/_downloads.py

@@ -69,7 +69,7 @@
 # The given commit will be downloaded from github and used as source for
 # all sample data
 SAMPLE_DATA_GH_REPO = 'OGGM/oggm-sample-data'
-SAMPLE_DATA_COMMIT = 'c631560e5a84008e95b5bca3809992e0d60472cc'
+SAMPLE_DATA_COMMIT = '7c4d9ddb9cd855646f63d66d39324865acbbba0b'
 
 GDIR_L1L2_URL = ('https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.4/'
                  'L1-L2_files/centerlines/')