Skip to content

Commit ad42322

Browse files
chloe-hancockchloe-hancock
authored
OGGM Calibration to RMSD functionality (#1827)
1 parent ee25633 commit ad42322

File tree

3 files changed

+550
-1
lines changed

3 files changed

+550
-1
lines changed

‎docs/whats-new.rst‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -61,8 +61,14 @@ Enhancements
6161
``climatic_mb_myr`` -> ``climatic_mb`` and ``flux_divergence_myr`` ->
6262
``flux_divergence`` (:pull:`1792`).
6363
By `Patrick Schmitt <https://github.com/pat-schmitt>`_
64+
<<<<<<< rmsd_calibration
65+
- Added a new way to calibrate the mass balance model with MB
66+
timeseries (:pull:`1827`).
67+
By `Chloe Hancock <https://github.com/chloe-hancock>`_
68+
=======
6469
- Fixed ``TimeoutError`` when running tests when urlopen hangs on DNS lookup (:pull:`1813`).
6570
By `Nicolas Gampierakis <https://github.com/gampnico>`_.
71+
>>>>>>> master
6672

6773
Bug fixes
6874
~~~~~~~~~

‎oggm/core/massbalance.py‎

Lines changed: 250 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,7 @@
1717
floatyear_to_date, date_to_floatyear, get_demo_file,
1818
monthly_timeseries, ncDataset, get_temp_bias_dataframe,
1919
clip_min, clip_max, clip_array, clip_scalar,
20-
weighted_average_1d, lazy_property, set_array_type)
20+
weighted_average_1d, lazy_property, set_array_type, rmsd)
2121
from oggm.exceptions import (InvalidWorkflowError, InvalidParamsError,
2222
MassBalanceCalibrationError)
2323
from oggm import entity_task
@@ -1517,6 +1517,255 @@ def mb_calibration_from_wgms_mb(gdir, **kwargs):
15171517
**kwargs)
15181518

15191519

1520+
@entity_task(log, writes=['mb_calib'])
1521+
def mb_calibration_to_rmsd(gdir, *,
1522+
ref_df=None,
1523+
write_to_gdir=True,
1524+
overwrite_gdir=False,
1525+
use_2d_mb=False,
1526+
calibrate_params=('melt_f',),
1527+
melt_f=None,
1528+
melt_f_min=None,
1529+
melt_f_max=None,
1530+
prcp_fac=None,
1531+
prcp_fac_min=None,
1532+
prcp_fac_max=None,
1533+
temp_bias=None,
1534+
temp_bias_min=None,
1535+
temp_bias_max=None,
1536+
mb_model_class=MonthlyTIModel,
1537+
filesuffix='',):
1538+
"""Determine the MB parameters by minimising RMSD to a reference timeseries
1539+
1540+
This calibrates the mass balance parameters using interannual
1541+
MB data from the WGMS data over a given period. This calibration uses
1542+
differential evolution to calibrate all given parameters to minimize
1543+
the RMSD as much as possible.
1544+
1545+
This function is useful to calibrate all three parameters at once,
1546+
on glaciers where WGMS or other in-situ observations are available.
1547+
This is achieved by minimising the RMSD between the reference MB
1548+
timeseries and the modelled MB timeseries over the period of available
1549+
observations. The minimisiation technique chosen here is differential
1550+
evolution, which is a global optimization technique that does not
1551+
require the function to be differentiable. This makes it
1552+
suitable for our problem, where the relationship between the parameters
1553+
and the MB timeseries can be complex and non-linear, and we are able
1554+
to calibrate all three parameters at once.
1555+
1556+
Note that this does not compute the apparent mass balance at
1557+
the same time - users need to run `apparent_mb_from_any_mb after`
1558+
calibration.
1559+
1560+
Parameters
1561+
----------
1562+
gdir : :py:class:`oggm.GlacierDirectory`
1563+
the glacier directory to calibrate
1564+
ref_df : pandas dataframe, required
1565+
the dataframe of annual mass balance values from the wgms data
1566+
(units: kg m-2 yr-1).
1567+
It is required here - if you want to use available observations,
1568+
write_to_gdir : bool
1569+
whether to write the results of the calibration to the glacier
1570+
directory. If True (the default), this will be saved as `mb_calib.json`
1571+
and be used by the MassBalanceModel class as parameters in subsequent
1572+
tasks.
1573+
overwrite_gdir : bool
1574+
if a `mb_calib.json` exists, this task won't overwrite it per default.
1575+
Set this to True to enforce overwriting (i.e. with consequences for the
1576+
future workflow).
1577+
use_2d_mb : bool
1578+
Set to True if the mass balance calibration has to be done of the 2D mask
1579+
of the glacier (for fully distributed runs only).
1580+
mb_model_class : MassBalanceModel class
1581+
the MassBalanceModel to use for the calibration. Needs to use the
1582+
same parameters as MonthlyTIModel (the default): melt_f,
1583+
temp_bias, prcp_fac.
1584+
calibrate_params : tuple
1585+
the parameter(s) that will be used in the calibration, it must be at least one of:
1586+
'melt_f', 'temp_bias', 'prcp_fac'. Defaults to ('melt_f',)
1587+
melt_f: float
1588+
the default value to use as melt factor (or the starting value when
1589+
optimizing MB). Defaults to cfg.PARAMS['melt_f'].
1590+
melt_f_min: float
1591+
the minimum accepted value for the melt factor during optimisation.
1592+
Defaults to cfg.PARAMS['melt_f_min'].
1593+
melt_f_max: float
1594+
the maximum accepted value for the melt factor during optimisation.
1595+
Defaults to cfg.PARAMS['melt_f_max'].
1596+
prcp_fac: float
1597+
the default value to use as precipitation scaling factor
1598+
(or the starting value when optimizing MB). Defaults to the method
1599+
chosen in `params.cfg` (winter prcp or global factor).
1600+
prcp_fac_min: float
1601+
the minimum accepted value for the precipitation scaling factor during
1602+
optimisation. Defaults to cfg.PARAMS['prcp_fac_min'].
1603+
prcp_fac_max: float
1604+
the maximum accepted value for the precipitation scaling factor during
1605+
optimisation. Defaults to cfg.PARAMS['prcp_fac_max'].
1606+
temp_bias: float
1607+
the default value to use as temperature bias (or the starting value when
1608+
optimizing MB). Defaults to 0.
1609+
temp_bias_min: float
1610+
the minimum accepted value for the temperature bias during optimisation.
1611+
Defaults to cfg.PARAMS['temp_bias_min'].
1612+
temp_bias_max: float
1613+
the maximum accepted value for the temperature bias during optimisation.
1614+
Defaults to cfg.PARAMS['temp_bias_max'].
1615+
filesuffix: str
1616+
add a filesuffix to mb_calib.json. This could be useful for sensitivity
1617+
analyses with MB models, if they need to fetch other sets of params for
1618+
example.
1619+
"""
1620+
1621+
# Param constraints
1622+
if melt_f_min is None:
1623+
melt_f_min = cfg.PARAMS['melt_f_min']
1624+
if melt_f_max is None:
1625+
melt_f_max = cfg.PARAMS['melt_f_max']
1626+
if prcp_fac_min is None:
1627+
prcp_fac_min = cfg.PARAMS['prcp_fac_min']
1628+
if prcp_fac_max is None:
1629+
prcp_fac_max = cfg.PARAMS['prcp_fac_max']
1630+
if temp_bias_min is None:
1631+
temp_bias_min = cfg.PARAMS['temp_bias_min']
1632+
if temp_bias_max is None:
1633+
temp_bias_max = cfg.PARAMS['temp_bias_max']
1634+
1635+
if not use_2d_mb:
1636+
fls = gdir.read_pickle('inversion_flowlines')
1637+
else:
1638+
# if the 2D data is used, the flowline is not needed.
1639+
fls = None
1640+
# get the 2D data
1641+
fp = gdir.get_filepath('gridded_data')
1642+
with xr.open_dataset(fp) as ds:
1643+
# 'topo' instead of 'topo_smoothed'?
1644+
heights = ds.topo_smoothed.data[ds.glacier_mask.data == 1]
1645+
widths = np.ones(len(heights))
1646+
1647+
# Climate period
1648+
ref_mb_years = ref_df.index.values
1649+
years = ref_mb_years
1650+
1651+
# Do we have a calving glacier?
1652+
cmb = calving_mb(gdir)
1653+
if cmb != 0:
1654+
raise NotImplementedError('Calving with geodetic MB is not implemented '
1655+
'yet, but it should actually work. Well keep '
1656+
'you posted!')
1657+
1658+
# Ok, regardless on how we want to calibrate, we start with defaults
1659+
if melt_f is None:
1660+
melt_f = cfg.PARAMS['melt_f']
1661+
1662+
if prcp_fac is None:
1663+
if cfg.PARAMS['use_winter_prcp_fac']:
1664+
# Some sanity check
1665+
if cfg.PARAMS['prcp_fac'] is not None:
1666+
raise InvalidWorkflowError("Set PARAMS['prcp_fac'] to None "
1667+
"if using PARAMS['winter_prcp_factor'].")
1668+
prcp_fac = decide_winter_precip_factor(gdir)
1669+
else:
1670+
prcp_fac = cfg.PARAMS['prcp_fac']
1671+
if prcp_fac is None:
1672+
raise InvalidWorkflowError("Set either PARAMS['use_winter_prcp_fac'] "
1673+
"or PARAMS['winter_prcp_factor'].")
1674+
1675+
if temp_bias is None:
1676+
temp_bias = 0
1677+
1678+
# Create the MB model we will calibrate
1679+
mb_mod = mb_model_class(gdir,
1680+
melt_f=melt_f,
1681+
temp_bias=temp_bias,
1682+
prcp_fac=prcp_fac,
1683+
check_calib_params=False)
1684+
1685+
# Check that the years are available
1686+
for y in years:
1687+
if not mb_mod.is_year_valid(y):
1688+
raise ValueError(f'year {y} out of the valid time bounds: '
1689+
f'[{mb_mod.ys}, {mb_mod.ye}]')
1690+
1691+
# Check that the calibrate params are valid
1692+
for param in calibrate_params:
1693+
if param not in ('melt_f', 'prcp_fac', 'temp_bias'):
1694+
raise InvalidParamsError("calibrate_params must be a tuple with any of "
1695+
"'melt_f', 'prcp_fac', 'temp_bias'")
1696+
1697+
# Set the bounds for the optimization
1698+
bounds = []
1699+
for param in calibrate_params:
1700+
if param == 'prcp_fac':
1701+
bounds.append((prcp_fac_min, prcp_fac_max))
1702+
elif param == 'melt_f':
1703+
bounds.append((melt_f_min, melt_f_max))
1704+
elif param == 'temp_bias':
1705+
bounds.append((temp_bias_min, temp_bias_max))
1706+
1707+
# Optimises all three mass balance parameters at the same time to minimize the RMSD between the simulated and reference MB timeseries
1708+
def rmsd_cost_function(x, *model_attrs: tuple):
1709+
for i, model_attr in enumerate(model_attrs):
1710+
setattr(mb_mod, model_attr, x[i])
1711+
1712+
if use_2d_mb:
1713+
sim_out = mb_mod.get_specific_mb(heights=heights, widths=widths, year=years)
1714+
else:
1715+
sim_out = mb_mod.get_specific_mb(fls=fls, year=years)
1716+
1717+
return rmsd(ref_df, sim_out)
1718+
1719+
try:
1720+
res = optimize.differential_evolution(rmsd_cost_function,
1721+
bounds=bounds,
1722+
tol=1e-8,
1723+
maxiter=5000,
1724+
args=(calibrate_params),
1725+
)
1726+
1727+
calib_params = res.x
1728+
1729+
# Assign parameters
1730+
for i, param in enumerate(calibrate_params):
1731+
if param == 'prcp_fac':
1732+
prcp_fac = calib_params[i]
1733+
elif param == 'melt_f':
1734+
melt_f = calib_params[i]
1735+
elif param == 'temp_bias':
1736+
temp_bias = calib_params[i]
1737+
1738+
except ValueError:
1739+
raise RuntimeError(f'{gdir.rgi_id}: could not minimise the rmsd. '
1740+
f'Try another technique.')
1741+
1742+
# Store parameters
1743+
df = gdir.read_json('mb_calib', allow_empty=True)
1744+
df['rgi_id'] = gdir.rgi_id
1745+
df['bias'] = 0
1746+
df['melt_f'] = melt_f
1747+
df['prcp_fac'] = prcp_fac
1748+
df['temp_bias'] = temp_bias
1749+
# What did we try to match?
1750+
df['reference_mb'] = ref_df.values.mean()
1751+
df['reference_period'] = str(ref_mb_years)
1752+
df['rmsd'] = res.fun
1753+
1754+
# Add the climate related params to the GlacierDir to make sure
1755+
# other tools cannot fool around without re-calibration
1756+
df['mb_global_params'] = {k: cfg.PARAMS[k] for k in MB_GLOBAL_PARAMS}
1757+
df['baseline_climate_source'] = gdir.get_climate_info()['baseline_climate_source']
1758+
# Write
1759+
if write_to_gdir:
1760+
if gdir.has_file('mb_calib', filesuffix=filesuffix) and not overwrite_gdir:
1761+
raise InvalidWorkflowError('`mb_calib.json` already exists for '
1762+
'this repository. Set `overwrite_gdir` '
1763+
'to True if you want to overwrite '
1764+
'a previous calibration.')
1765+
gdir.write_json(df, 'mb_calib', filesuffix=filesuffix)
1766+
return df
1767+
1768+
15201769
@entity_task(log, writes=['mb_calib'])
15211770
def mb_calibration_from_geodetic_mb(gdir, *,
15221771
ref_period=None,

0 commit comments

Comments
 (0)