Assimilating records whose PSMs cannot be calibrated

Contents

Assimilating records whose PSMs cannot be calibrated#

In some cases, we may need to assimilate records whose PSMs cannot be calibrated against instrumental observations due to non-overlapping timespans. This tutorial shows how to specify the required parameters and assimilate those records without calling the calibration step of their PSMs.

[1]:

%load_ext autoreload
%autoreload 2

import cfr
print(cfr.__version__)

Create a reconstruction job object#

[2]:

job = cfr.ReconJob(verbose=True)

>>> job.configs:
{}

Load a proxy database#

[3]:

job.load_proxydb('PAGES2kv2')

Filter the proxy database#

[4]:

job.filter_proxydb(by='ptype', keys=['coral'])
# the .plot() method will generate a static map utilizing the Cartopy library
fig, ax = job.proxydb.plot(plot_count=True)

../_images/notebooks_lmr-assim-no-calib_7_0.png

Annualize/seasonalize the proxy database#

[5]:

job.annualize_proxydb(months=[12, 1, 2], ptypes=['coral'], verbose=True)
fig, ax = job.proxydb.plot()

>>> job.configs["annualize_proxydb_months"] = [12, 1, 2]
>>> job.configs["annualize_proxydb_ptypes"] = ['coral']

Annualizing ProxyDatabase: 100%|██████████| 104/104 [00:05<00:00, 18.90it/s]

>>> 104 records remaining
>>> job.proxydb updated

../_images/notebooks_lmr-assim-no-calib_9_3.png

Load climate model priors#

[6]:

job.load_clim(
    tag='prior',
    path_dict={
        'tas': 'iCESM_past1000historical/tas',
    },
    anom_period=(1951, 1980),
    verbose=True,
)

>>> job.configs["prior_path"] = {'tas': './data/tas_sfc_Amon_iCESM_past1000historical_085001-200512.nc'}
>>> job.configs["prior_anom_period"] = (1951, 1980)
>>> job.configs["prior_lat_name"] = lat
>>> job.configs["prior_lon_name"] = lon
>>> job.configs["prior_time_name"] = time
>>> prior variables ['tas'] loaded
>>> job.prior created

Forward modeling without calling the PSMs#

This is the key step of this tutorial. Below we loop over each proxy record and perform forward modeling steps:

Get the climate data at the nearest grid point of the record.
Annualize the climate data to get the same seasonality as the record.
Create the pobj.pseudo object for each record, with pobj.pseudo.time and pobj.pseudo.value specified. The .value can be generated based on uncalibrated but known relationships (e.g., linear with a certain slope).
Add the tag “calibrated” so that the workflow will not complain.
Set the observation error based on estimation.

[24]:

for pid, pobj in job.proxydb.records.items():
    pobj.get_clim(job.prior['tas'], tag='model')
    clim_ann = pobj.clim['model.tas'].annualize(months=[12, 1, 2])
    pobj.pseudo = pobj.copy()
    pobj.pseudo.time = clim_ann.da.time.values
    pobj.pseudo.value = clim_ann.da.values # can be arbitrary relationship
    pobj.pseudo.seasonality = [12, 1, 2]
    pobj.tags.add('calibrated')
    pobj.R = 0.5 # obs err; should be estimated ahead

job.proxydb.nrec_tags("calibrated") # confirm that all the records have been "calibrated"

Annualize/seasonalize the climate model prior#

[12]:

job.annualize_clim(tag='prior', verbose=True, months=[12, 1, 2])
fig, ax = job.prior['tas'][-1].plot()

>>> job.configs["prior_annualize_months"] = [12, 1, 2]
>>> Processing tas ...
>>> job.prior updated

../_images/notebooks_lmr-assim-no-calib_15_1.png

Crop the prior field#

[14]:

job.crop_clim(tag='prior', lat_min=-35, lat_max=35, verbose=True)
fig, ax = job.prior['tas'][-1].plot()

>>> job.configs["prior_lat_min"] = -35
>>> job.configs["prior_lat_max"] = 35
>>> job.configs["prior_lon_min"] = 0
>>> job.configs["prior_lon_max"] = 360
>>> Processing tas ...

../_images/notebooks_lmr-assim-no-calib_19_1.png

Data assimilation#

[26]:

job.run_da_mc(
    save_dirpath='./recons/lmr-no-calib',
    recon_period=(1801, 2000),
    recon_seeds=list(range(1, 2)), verbose=True)

>>> job.configs["recon_period"] = (1801, 2000)
>>> job.configs["recon_seeds"] = [1]
>>> job.configs["save_dirpath"] = ./recons/lmr-no-calib
>>> seed: 1 | max: 1

KF updating: 100%|██████████| 200/200 [00:13<00:00, 15.07it/s]

>>> Reconstructed fields saved to: ./recons/lmr-no-calib/job_r01_recon.nc
>>> DONE! Total time spent: 0.24 mins.

Validation#

[27]:

res = cfr.ReconRes('./recons/lmr-no-calib', verbose=True)

>>> res.paths:
['./recons/lmr-no-calib/job_r01_recon.nc']

[28]:

res.load(['nino3.4', 'tas'], verbose=True)

>>> ReconRes.recons["nino3.4"] created
>>> ReconRes.da["nino3.4"] created
>>> ReconRes.recons["tas"] created
>>> ReconRes.da["tas"] created

[29]:

import scipy.io as sio
import numpy as np

data = sio.loadmat('./data/BC09_NINO34.mat')
start_yr = 1874
end_yr = 2000
nyr = end_yr-start_yr+1
nino34 = np.zeros(nyr)
for i in range(nyr):
    nino34[i] = np.mean(data['nino34'][12+i*12:14+i*12])

ref_value = nino34
ref_time = np.arange(start_yr, end_yr+1)
ref_name = 'BC09'

[30]:

fig, ax = res.recons['nino3.4'].compare(
    ref_time, ref_value, ref_name, timespan=(1874, 1900)).plot_qs()
ax.set_xlim(1800, 2000)
ax.set_ylabel('NINO3.4 [K]')

[30]:

Text(0, 0.5, 'NINO3.4 [K]')

../_images/notebooks_lmr-assim-no-calib_26_1.png

The above reconstruction shows that, without proper calibration of the PSMs, the posterior can become weired due to the inaccurate estimate of the observation error and the inaccurate forward modeling relationship.