<jupyter_start><jupyter_text>Compute the annual temperature anomaly from Temperature12kThis notebooks computes GMST using Generalized Additive Models (GAM). The main output is a netCDF file `annual-with-single-seas-gam-band-means.nc` with the temperature anomaly with reference to the 5000-3000 year BP period for the six latitudinal bands -90 to -60, -60 to -30, -30 to 0, 0 to 30, 30 to 60 and 60 to 90 degrees North.This notebook is divided into four sections.1. [Read the LiPD files](lipds)2. [Transform them into a combined dataframe](pandas)3. [Compute the grid cell anomalies](cell-anomalies)4. [Compute the band means](band-anomalies)Note that you have to set the [version for the database](db_version) manually.<jupyter_code>import pandas as pd
import numpy as np
import os.path as osp
import os
from urllib import request
import zipfile
import pickle
from openpyxl import load_workbook
import matplotlib.pyplot as plt
import warnings
import xarray as xr
import multiprocessing as mp
import psyplot.data as psyd
import gam_ensemble as ge
import pyleogrid as pg
from distributed import Client
import psyplot.project as psy
from dask_jobqueue import SLURMCluster
from matplotlib.animation import FuncAnimation
import contextlib
from sklearn.neighbors import BallTree
from latlon_utils import get_climate
import lipd<jupyter_output><empty_output><jupyter_text>Read the LiPD files In this section we read the temperature from the LiPD files. This involves some [filtering](filter) to only keep the necessary temperature series. Note however, that samples with NaN are filtered [down below](filter-nan).Read in the LipD data from http://lipdverse.org/globalHolocene/current_versionYou should set the latest version here manually!<jupyter_code>%matplotlib inline
db_version = '1_0_0'
lipd_dir = f'../lipds_{db_version}'<jupyter_output><empty_output><jupyter_text>Download the LiPD file if the input directory does not yet exist<jupyter_code>%%time
if not osp.exists(lipd_dir):
    os.makedirs(lipd_dir)
    zipped = f'Temp12k_v{db_version}_LiPD.zip'
    uri = f'https://www1.ncdc.noaa.gov/pub/data/paleo/reconstructions/climate12k/temperature/version{db_version.replace("_", ".")}/{zipped}'
    target = osp.join(lipd_dir, zipped)
    print('downloading ' + uri)
    request.urlretrieve(uri, target)
    with zipfile.ZipFile(target) as f:
        f.extractall(lipd_dir)
@contextlib.contextmanager
def remember_cwd():
    """Context manager to switch back to the current working directory

    Usage::

        with remember_cwd():
            os.chdir('test')
            print(os.getcwd())  # test
        print(os.getcwd())      # test/.."""
    curdir = os.getcwd()
    try:
        yield
    except:
        raise
    finally:
        os.chdir(curdir)<jupyter_output><empty_output><jupyter_text>Load the data from the LiPD files it is has not yet been serialized as a pickled object.<jupyter_code>%%time
pickle_file = f'{lipd_dir}/lipds.pkl'
if osp.exists(pickle_file):
    with open(pickle_file, 'rb') as f:
        data = pickle.load(f)
else:
    with remember_cwd():
        os.chdir(f'{lipd_dir}/')
        data = lipd.readLipd('.')
    with open(pickle_file, 'wb') as f:
        pickle.dump(data, f)<jupyter_output>CPU times: user 137 ms, sys: 55.8 ms, total: 193 ms
Wall time: 192 ms<jupyter_text>Extract the individual time-series from the LiPD files<jupyter_code>%%time
all_series = lipd.extractTs(data)<jupyter_output>extracting paleoData...
extracting: GodziszewskieLake.Miotk.1986
extracting: Hudson.Clegg.2011
extracting: ParkPond1.Lynch.1998
extracting: ElGygytgynCrater.Schwamborn.2006
extracting: Dalmutladdo.Bjune.2004
extracting: GeoB10043_3.Setiawan.2015
extracting: MD99_2341.Eynaud.2009
extracting: CN01.Cheung.2017
extracting: HumberPond5.McAndrews.1989
extracting: TauldintreBrazi.Toth.2015
extracting: LagunaSaladaChiprana.EPD
extracting: PrazRodet.Shotyk.1997
extracting: SaladaPequena.EPD
extracting: Leviathan.Lachniet.2014
extracting: BloodPond.Marsicek.2013
extracting: Khatanga-12Lake.Syrykh.2017
extracting: rainbow.Clegg.2011
extracting: M77_2_024_5.Salvatteci.2019
extracting: GeoB7926_2.Kim.2012
extracting: LagoMoreno.Markgraf.2002
extracting: PolandSpring.Jacobson.Ja
extracting: MSM05-712.Muller.2013
extracting: Gammelheimvatnet.Seppa.2009
extracting: WCA3BGumboLimboMarsh.Willard.2001
extracting: NussbaumerSeen.Rsch.1995
extracting: DevonIslandGlacier.McAndrews.1984
extracting: klotjarne[...]<jupyter_text>Filter all the series to extract the calibrated temperature series<jupyter_code>filtered_ts_temp12k = lipd.filterTs(all_series,'paleoData_inCompilation == Temp12k')
temperatures = lipd.filterTs(filtered_ts_temp12k,'paleoData_units == degC')<jupyter_output>Found 1319 matches from 5088 columns
Found 1276 matches from 1319 columns<jupyter_text>Conversion to pandasNow that we have the time-series, we transform them to a [pandas.DataFrame](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.html). This allows us to use groupby and other database measures. For this purpose, we first define a few meta data attributes that should be kept for every record.- `meta_cols` is the attribute in the dictionary of the individual time series. If the attribute is within a key or list, you can separate them by a single `.`.     `paleoData_interpretation.0.seasonalityGeneral`     for instance will be interpreted as    `d['paleoData_interpretation'][0]['seasonalityGeneral']`- `meta_names` is the name of the corresponding attribute in `meta_cols` how it should be saved in the dataframe<jupyter_code>meta_cols = ['geo_meanLon', 'geo_meanLat', 'dataSetName', 'paleoData_variableName', 'paleoData_TSid', 
             'paleoData_interpretation.0.seasonalityGeneral', 'paleoData_datum', 'paleoData_proxy']
meta_names = ['lon', 'lat', 'dataSetName', 'variableName', 'TSid', 'seasonality', 'datum', 'proxy']<jupyter_output><empty_output><jupyter_text>Then every dictionary in `temperatures` is converted to a pandas datafram<jupyter_code>frames = []
failures = []
for d in pg.utils.log_progress(temperatures):
    if 'age' in d:
        try:
            s = pd.Series(
                np.asarray(d['paleoData_values'], dtype=float),
                index=pd.Index(np.asarray(d['age'], dtype=float), name='age'),
                name='temperature')
        except ValueError:
            failures.append(d)
        else:
            # keep NaN of ages, these are filtered out later
            s = s[(~np.isnan(s.index)) & (s.index > -70 ) & (s.index < 14000)]
            # if the time series has (in average) more than one sample per
            # decade, we use decadal averages
            if s.size > 1470:
                s = s.groupby(lambda i: 5 + i - (i % 10)).mean()
                s.index.name = 'age'
            s.index.name = 'age'
            frame = s.to_frame().reset_index()
            for i, name in enumerate(meta_cols):
                sub = d
                vals = name.split('.')
                name = vals[-1]
                for v in vals[:-1]:
                    sub = sub[v if isinstance(sub, dict) else int(v)]
                frame[meta_names[i]] = sub.get(name, np.nan)
            frames.append(frame)
print(f"Failed to read {len(failures)} time series")<jupyter_output>Starting at 2020-04-17 20:49:06.189962
|████████████████████████████████████████████████████████████████████████████████| 100.0% Time left: 0.000 minutes
Finished at 2020-04-17 20:49:11.346282. Time needed: 0.086 minutes
Failed to read 0 time series<jupyter_text>Now we concatenate everything into one large dataframe<jupyter_code>df = pd.concat(frames, ignore_index=True)
df.head()<jupyter_output><empty_output><jupyter_text>Estimation of uncertainties The age uncertainty is estimated such, that it is 50 years at 0 BP and 250 years at 12k BP<jupyter_code>df['age_unc'] = 50 + df.age * (250 - 50) / 12000<jupyter_output><empty_output><jupyter_text>The temperature uncertainty is estimated per proxy and per season<jupyter_code>uncertainties = pd.read_excel(
    '../data/proxy-uncertainties.xlsx', sheet_name='Summary', 
    index_col='proxy', usecols=[0, 1, 2, 3])
uncertainties
uncertainty_map = {
    'dinocyst': 'other microfossils/dinocyst',
    'pollen': 'pollen', 
    'chironomid': 'chironomid',
    'alkenone': 'alkenone',
    'GDGT': 'GDGT (MBT/CBT as well as BrGDGT fractional abundance)',
    'Mg/Ca': 'MgCa',
    'macrofossils': 'macrofossils',
    'bubble frequency': 'bubble frequency',
    'borehole': 'borehole',
    'C15 fatty alcohols': 'C15 fatty alcohols',
    'd18O': 'd18O',
    'radiolaria': 'other microfossils/radiolaria',
    'foraminifera': 'other microfossils/foraminifera',
    'dD': 'dD',
    'd13C': 'd13C',
    'diatom': 'other microfossils/diatoms',
    'isotope diffusion': 'isotope diffusion',
    'TRW': 'TRW',
    'hybrid-ice': 'hybrid-ice',
    'deterium excess': 'deterium excess',
    'gas': 'gas',
    'chlorophyll': 'chlorophyll',
    'melt layer': 'melt layer',
    '3-OH-Fatty Acids': '3-OH-Fatty Acids',
    'TEX86': 'GDGT (Tex86)',
    'LDI': 'LDI',
    'long chain diol': 'long chain diol'
}
uncertainty_map_inv = {val: key for key, val in uncertainty_map.items()}
df['seas'] = np.nan
seasonality = df.seasonality.fillna('')
df.loc[seasonality.str.startswith('summer'), 'seas'] = 'summer'
df.loc[seasonality.str.startswith('winter'), 'seas'] = 'winter'
df.loc[seasonality.str.startswith('annual'), 'seas'] = 'annual'
last = uncertainties.isnull().all(axis=1)
uncertainties = uncertainties.iloc[:np.argmax(last.values)]
# default uncertainty: 75th percentile of all values without d18O
df['temp_unc'] = 1.975877719
for proxy, row in uncertainties.iterrows():
    for seas, val in row[row.notnull().values].items():
        mask = (df.proxy == uncertainty_map_inv[proxy]) & (df.seas == seas)
        df.loc[mask.values, 'temp_unc'] = val
import seaborn as sns
ax = sns.stripplot('proxy', 'temp_unc', data=df)
for t in ax.get_xticklabels():
    t.set_rotation(90)<jupyter_output><empty_output><jupyter_text>Now we filter out samples that have unrealistic temperatures, or no seasonality .<jupyter_code>df.loc[~df.temperature.between(-200, 200).values, 'temperature'] = np.nan
for col in ['seasonality', 'temperature']:
    print(f'NaN for {col}: {df[col].isnull().sum()}')
print(df.shape)
df.dropna(subset=['seasonality', 'temperature'], inplace=True)
print(df.shape)<jupyter_output>(121457, 13)
(110816, 13)<jupyter_text>Computation of the grid cell anomalies Now that the input data is filtered and ready to use, we compute the anomalies per grid cell. For this purpose, we use the `GAMEnsemble` class defined in [gam_ensemble.py](gam_ensemble.py).We use all annual records, and, if there is only a summer or winter record available in the LiPD dataset, we use this one as well.<jupyter_code>with psyd.open_dataset('../data/equal-area-grid-N4000.nc') as _ds:
    output_ds = _ds.load()
area = output_ds['cell_area']
output_ds['time'] = time = xr.Variable(
    ('time', ), np.arange(12000, -50, -100),
    attrs = {'long_name': 'Age', 'units': 'years since 1950-01-01'})
output_ds['temperature'] = xr.Variable(
    ('time', area.dims[0]), np.full((time.size, area.size), np.nan),
    encoding = dict(coordinates='clat clon'))

ensemble = ge.GAMEnsemble(
    df[(df.seasonality == 'annual') | 
       (df.seasonality == 'summerOnly') | 
       (df.seasonality == 'winterOnly')].set_index('age'), 
    target=output_ds.psy['temperature'],
    climate='temperature', ds_id='TSid')<jupyter_output><empty_output><jupyter_text>Generate the ensembleFor age, we use a gibbs sampling algorithm to generate an ensemble of 500 ages per sample. For temperature, we use a normal distribution defined by the temperature of the sample and the uncertainty [estimated above](temp-unc).<jupyter_code>%%time
ds = ensemble.input_data
ages = ensemble.sample_ages(size=500, use_dask=False)
ds['age_ensemble'] = ages.load()
ds['temperature_ensemble'] = ages.rename('temperature_ensemble').copy(
    data=np.random.normal(ds.temperature, ds.temp_unc,
                          size=ages.shape))<jupyter_output>CPU times: user 47.2 s, sys: 241 ms, total: 47.5 s
Wall time: 47.5 s<jupyter_text>We save the input data to share it with others<jupyter_code>output_dir = f'../output_{db_version}'
!mkdir -p {output_dir}
to_keep = [
    'dataSetName',
    'variableName',
    'seasonality',
    'datum',
    'proxy',
    'age_unc',
    'seas',
    'temp_unc',
    'age_ensemble',
    'temperature_ensemble'
]

key = 'annual-with-single-seas'

ds = ensemble.input_data[to_keep]
encoding = {v: dict(zlib=True, complevel=4)
            for v in ds.variables}
for var in ds.variables:
    if np.issubdtype(ds[var].dtype, np.number):
        encoding[var]['least_significant_digit'] = 5
ofile = f'{output_dir}/{key}_input.nc'
ds.rename_dims(age='sampleId').to_netcdf(
    ofile, encoding=encoding)<jupyter_output><empty_output><jupyter_text>Computation of grid cell anomaliesThe heavy part comes now, we compute the grid cell anomalies and save them to `annual-with-single-seas-gam-ref-5000-3000-bp.nc`.<jupyter_code>n_workers = 64

cluster = SLURMCluster(cores=1, dashboard_address=":33913", walltime="03:00:00")
cluster.scale(n_workers)

try:
    with Client(cluster) as client:
        with warnings.catch_warnings() as c:
            warnings.filterwarnings('ignore', '\s*indexing past lexsort depth')
            output = ensemble.predict(
                size=500, client=client, use_gam=True, return_counts=100, max_time_diff=14000,
                modern_young=3000, modern_old=5000)
            output['time'].values *= -1
            output.to_netcdf(f'{output_dir}/{key}-gam-ref-5000-3000-bp.nc')
finally:
    cluster.close()<jupyter_output>Starting at 2020-04-17 20:50:23.437711
|████████████████████████████████████████████████████████████████████████████████| 100.0% Time left: 0.042 minutes
Finished at 2020-04-17 21:04:09.345651. Time needed: 13.765 minutes<jupyter_text>We only want to use the gam prediction if- there has been more than 100 samples between 5k and 3k BP in the ensemble- there were at least 100 samples in the surrounding 300 years of the time of prediction<jupyter_code># valid = output.nsamples.sel(time=slice(-5000, -3000)).sum(dim='time')
valid = output.nsamples.where((output.time > -5000) & (output.time < -3000), drop=True).sum(dim='time')
surrounding = output.nsamples.rolling(time=3).sum()
output = output.where(valid > 100)
output = output.where(surrounding > 100)<jupyter_output><empty_output><jupyter_text>Computation of latitudinal bandsIn the last part of this analysis we compute the mean of the individual latitudinal bands. We are using an equal-area grid here, as such, we do not have to worry about the latitudinal weighting.<jupyter_code>bands = []
with warnings.catch_warnings() as c:
    warnings.filterwarnings('ignore', '\s*Mean of empty slice')
    for lat in pg.utils.log_progress(np.arange(-90, 90, 30)):
        idx = np.where(output.clat - (output.clat % 30) == lat)[0]
        sel = output.isel(cell=idx)
        band = sel.mean('cell')
        
        # compute the sum for the number of samples and the cell_area
        band['nsamples'] = sel.nsamples.sum('cell')
        band['cell_area'] = sel.cell_area.sum('cell')
        band['lat'] = ((), lat + 15)
        bands.append(band)
bands = xr.concat(bands, dim='lat')<jupyter_output>Starting at 2020-04-17 21:04:21.660716
|████████████████████████████████████████████████████████████████████████████████| 100.0% Time left: 0.011 minutes
Finished at 2020-04-17 21:04:25.294907. Time needed: 0.061 minutes<jupyter_text>Compute the standard deviation from the ensemble<jupyter_code>bands['temperature_std'] = bands.temperature_samples.std('ens')<jupyter_output>/gpfs/home/sommerp/miniconda3/envs/temperature12k/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1667: RuntimeWarning: Degrees of freedom <= 0 for slice.
  keepdims=keepdims)<jupyter_text>Plot the anomalies of the individual bands<jupyter_code>sp = bands.psy.plot.lineplot(
    name=[[['temperature', 'temperature_std']]],
    lat=list(range(6)), legendlabels='%(lat)s',
    xlabel='Time [years since 1950-01-01]',
    ylabel='Temperature anomaly to 5000-3000 yr cal BP [°C]',
    xlim='minmax', 
    legend=dict(loc='upper center', ncol=6, bbox_to_anchor=(0.5, 1.3),
                title='Center of latitudinal band [°N]'))

plt.gcf().set_dpi(300)<jupyter_output><empty_output><jupyter_text>Compute the samples per band that have been used. There is an nsamples variable already in the file, but this one is based on the entire ensemble.<jupyter_code>nsamples = ensemble.input_data.age.groupby_bins('lat', range(-90, 91, 30)).apply(
    lambda da: da.groupby_bins(da, np.arange(-50, 12051, 100)).apply(
        lambda da: xr.DataArray(np.unique(da.TSid).size)))
nsamples['lat_bins'] = ('lat_bins', [i.mid for i in nsamples.lat_bins.values])
nsamples['age_bins'] = ('age_bins', [i.mid for i in nsamples.age_bins.values])
nsamples.to_dataset(name='nsamples').psy.plot.lineplot(
    lat_bins=list(range(6)), legendlabels='%(lat_bins)s',
    ylabel='Number of individual records', xlabel='Age yr cal BP',
    xlim='minmax', 
    legend=dict(loc='upper center', ncol=6, bbox_to_anchor=(0.5, 1.3),
                    title='Center of latitudinal band [°N]'))
plt.gcf().set_dpi(250)
plt.gca().invert_xaxis()
bands = bands.rename({'nsamples': 'nsamples_ens'})
bands['nsamples'] = nsamples.rename({'lat_bins': 'lat', 'age_bins': 'age'})
bands.to_netcdf(f'{output_dir}/annual-with-single-seas-gam-band-means.nc',
                encoding={var: {'complevel': 4, 'zlib': True} 
                          for var in bands.variables})<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13805169883734056, 'fraction_numerical': 0.037419573315272604, 'mean_word_length': 4.375910194174757, 'pattern_counts': {'":': 1, '<': 79, '<?xml version=': 0, '>': 81, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '88466', 'license': 'MIT', 'n_tokens_mistral': 6774, 'n_tokens_neox': 5997, 'n_words': 1505, 'path': 'ScientificDataAnalysis/GAM_frozen/scripts/compute-annual-gmst.ipynb', 'repo_name': 'subond/Temperature12k', 'size': '17718'}