<jupyter_start><jupyter_text>Constrain the ensembleUsing:- temperature- ocean heat uptake- CO2 concentration<jupyter_code>from climateforcing.utils import check_and_download
from matplotlib.lines import Line2D
import matplotlib.pyplot as pl
import pandas as pd
import numpy as np
SAMPLES = 100000
pl.rcParams['figure.figsize'] = (9/2.54, 9/2.54)
pl.rcParams['font.size'] = 9
pl.rcParams['font.family'] = 'Arial'
pl.rcParams['xtick.direction'] = 'out'
pl.rcParams['xtick.minor.visible'] = True
pl.rcParams['ytick.minor.visible'] = True
pl.rcParams['ytick.right'] = True
pl.rcParams['xtick.top'] = True
pl.rcParams['figure.dpi'] = 300<jupyter_output><empty_output><jupyter_text>Get observational data<jupyter_code># AR6 temperature dataset
check_and_download(
    "https://github.com/chrisroadmap/ar6/raw/main/data_input/observations/AR6%20FGD%20assessment%20time%20series%20-%20GMST%20and%20GSAT.xlsx",
    "../data_input/AR6_FGD_assessment_time_series_-_GMST_and_GSAT.xlsx"
)
xl = pd.read_excel('../data_input/AR6_FGD_assessment_time_series_-_GMST_and_GSAT.xlsx', skiprows=1, skipfooter=28, engine="openpyxl")
temp_gmst=xl['4-set mean'].values
temp_year=xl['Unnamed: 0'].values
pl.plot(temp_year, temp_gmst)

temp_gmst[-1]
# OHC dataset
check_and_download(
    "https://github.com/chrisroadmap/ar6/raw/main/data_input/observations/AR6_OHC_ensemble_FGDprelim.csv",
    "../data_input/AR6_OHC_ensemble_FGDprelim.csv"
)
OHC_df = pd.read_csv("../data_input/AR6_OHC_ensemble_FGDprelim.csv", skiprows=1)
OHUobs = OHC_df['Central Estimate Full-depth'].values
OHUyears = OHC_df['Year'].values
OHU90 = OHC_df['Full-depth Uncertainty (1-sigma)'].values
# rcmip concentrations files
check_and_download(
    'https://rcmip-protocols-au.s3-ap-southeast-2.amazonaws.com/v5.1.0/rcmip-concentrations-annual-means-v5-1-0.csv',
    '../data_input_large/rcmip-concentrations-annual-means-v5-1-0.csv'
)
conc_all = pd.read_csv('../data_input_large/rcmip-concentrations-annual-means-v5-1-0.csv')
co2_obs = conc_all[(conc_all['Scenario']=='ssp245')&(conc_all['Region']=='World')&(conc_all.Variable.str.endswith('CO2'))].loc[:,'1750':'2018'].values.squeeze()<jupyter_output><empty_output><jupyter_text>Constraining starts here<jupyter_code>C_CO2 = np.load('../data_output_large/C_CO2_unconstrained.npy')
T = np.load('../data_output_large/T_unconstrained.npy')
OHU = np.load('../data_output_large/OHU_unconstrained.npy')
def rmse(obs, mod):
    return np.sqrt(np.sum((obs-mod)**2)/len(obs))
## simple criterion: RMSE of temperature 1850-2014 < 0.135
rmse_temp = np.zeros((SAMPLES))
for i in range(SAMPLES):
    rmse_temp[i] = rmse(temp_gmst[:165], (T[100:265,i]-np.mean(T[100:151,i], axis=0)))
accept_temp=(rmse_temp<0.135)
print(np.sum(accept_temp))
valid_temp = np.arange(SAMPLES, dtype=int)[accept_temp]
accept_ohu = np.logical_and(
    ((OHUobs[-1]-OHUobs[0]) - np.sqrt(OHU90[0]**2+OHU90[-1]**2)) < 0.90*(OHU[268,:]-OHU[221,:])*1e-21,
    0.90*(OHU[268,:]-OHU[221,:])*1e-21 < ((OHUobs[-1]-OHUobs[0]) + np.sqrt(OHU90[0]**2+OHU90[-1]**2))
)
valid_ohu = np.arange(SAMPLES, dtype=int)[accept_ohu]
print(np.sum(accept_ohu))
accept_co2 = np.logical_and( (co2_obs[264] - 0.36) < C_CO2[264], C_CO2[264] < (co2_obs[264] + 0.36) )
print(np.sum(accept_co2))
accept_ssp245 = accept_temp * accept_ohu * accept_co2
print(np.sum(accept_ssp245))
# Also get SSP2-4.5
ssp245 = pd.read_csv('../data_output/ssp245.csv', index_col=0)
OFFSET = 0
PERIOD = slice(100,151)
ssp245 = ssp245.sub(ssp245.loc[:,'1850-07-01':'1900-07-01'].mean(axis=1), axis=0) + OFFSET
pl.plot(np.arange(1850, 2101), ssp245.T, color='red', lw=0.1)

pl.fill_between(
    np.arange(1850, 2101),
    np.percentile(T[100:,:]-np.mean(T[PERIOD,:], axis=0) + OFFSET, 5, axis=1),
    np.percentile(T[100:,:]-np.mean(T[PERIOD,:], axis=0) + OFFSET, 95, axis=1), 
    color='lightgray', lw=1, alpha=0.5,
    label='Unconstrained 90% range'
)
pl.fill_between(
    np.arange(1850, 2101),
    np.percentile(T[100:,accept_ssp245]-np.mean(T[PERIOD,accept_ssp245], axis=0) + OFFSET, 5, axis=1),
    np.percentile(T[100:,accept_ssp245]-np.mean(T[PERIOD,accept_ssp245], axis=0) + OFFSET, 95, axis=1), 
    color='gray', lw=1, alpha=0.5,
    label='Constrained 90% range'
)
pl.plot(
    np.arange(1850, 2101), 
    np.percentile(T[100:,accept_ssp245]-np.mean(T[PERIOD,accept_ssp245], axis=0) + OFFSET, 50, axis=1),
    color='k',
    label='Best estimate projection'
)
pl.plot(np.arange(1850, 2021), temp_gmst, color='blue', label='Observations', lw=1.0)

line = Line2D([0], [0], label='CMIP6 models', color='red', lw=0.5)

pl.title('SSP2-4.5 temperature projections')
pl.ylabel('°C warming since 1850-1900')
pl.xlim(1850, 2100)
pl.ylim(-1, 4.5)

handles, labels = pl.gca().get_legend_handles_labels()
handles.extend([line])
pl.legend(handles=handles, loc='upper left')

pl.tight_layout()
pl.savefig('../figures/constraining.png')
pl.savefig('../figures/constraining.pdf')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18925518925518925, 'fraction_numerical': 0.0757020757020757, 'mean_word_length': 6.11287988422576, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 13, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '284558', 'license': 'MIT', 'n_tokens_mistral': 2361, 'n_tokens_neox': 2103, 'n_words': 260, 'path': 'notebooks/03_constraining.ipynb', 'repo_name': 'chrisroadmap/cb-emulators', 'size': '4914'}