The BAYSPAR TEX₈₆ calibration, in Python¶

Warning

This project is under heavy development.

baysparpy is an Open Source Python package for TEX86 calibration.

The package is based on the original BAYSPAR (BAYesian SPAtially-varying Regression) software for MATLAB.

Documentation¶

Quick Overview¶

BAYSPAR (BAYesian SPAtially-varying Regression) is a Bayesian calibration model for the TEX₈₆ temperature proxy.

There are two calibrations — one for sea-surface temperatures (SSTs) and one for subsurface (0-200 m) temperatures (subT) — and two methods to infer or predict past temperatures. The “standard” version of the model draws calibration parameters from the model gridpoint nearest to the core site to predict temperatures. The “Deep-Time” or “analog” version searches the calibration dataset for coretop TEX₈₆ values similar to those in the time series (within a user-set search tolerance) and uses the parameters from those analog locations to predict SSTs. The SST data in the calibration model are the statistical mean data from the World Ocean Atlas 2009.

For further details, refer to the original publication in Geochimica et Cosmochimica Acta, and the updated calibration publication in Scientific Data.

The original BAYSPAR MATLAB code is also available online.

What’s New¶

v0.0.4¶

Enhancements¶

Minor improvements to documentation.

Bug fixes¶

v0.0.3¶

Enhancements¶

Many minor changes to documentation.

v0.0.2¶

Enhancements¶

Prediction.location refactored to Prediction.latlon (Issue #8).
Removed unused Draws.alpha_samples_field and Draws.beta_samples_field to cut down package size. Also, we dropped the maximum MCMC parameter ensemble size to 10,000 to further reduce size (Issue #1).
Many assert errors replaced with full errors (Issue #7).
Improved pip and conda install instructions (Issue #6).

Bug fixes¶

Minor corrections to documentation and docstrings (Issue #10, Issue #5).
Trim burnin from tau2 MCMC params (Issue #3)

v0.0.1¶

This is the first release of baysparpy.

Examples¶

To start things off, import a few basic tools, including bayspar:

In [1]: import numpy as np

In [2]: import matplotlib.pyplot as plt

In [3]: import bayspar as bsr

There are a few options when it comes to prediction with bayspar. Below, we use example data - included with the package - to walk through each type of prediction.

Standard prediction¶

We can access the example data with get_example_data() and use the returned stream with pandas.read_csv() or numpy.genfromtxt():

In [4]: example_file = bsr.get_example_data('castaneda2010.csv')

In [5]: d = np.genfromtxt(example_file, delimiter=',', names=True)

This dataset (from Castañeda et al. 2010) has two columns giving sediment age (calendar years BP) and TEX₈₆.

In [6]: d['age'][:5]
Out[6]: array([142., 284., 375., 466., 558.])

In [7]: d['tex86'][:5]
Out[7]: array([0.662, 0.662, 0.652, 0.66 , 0.664])

We can make a “standard” prediction of sea-surface temperature (SST) with predict_seatemp():

In [8]: prediction = bsr.predict_seatemp(d['tex86'], lon=34.0733, lat=31.6517,
   ...:                                  prior_std=6, temptype='sst')
   ...: 

The TEX₈₆ data and site position are passed to the prediction function. We also give standard deviation for the prior SST distribution and specify that SST are the target variable. A Prediction instance is returned, which we can poked and prod with direct queries or using a few built-in functions:

In [9]: bsr.predictplot(prediction)
Out[9]: <AxesSubplot:>

You might have noticed that predictplot() returns an matplotlib.Axes instance. This means we can catch it and then modify the plot as needed:

In [10]: ax = bsr.predictplot(prediction, x=d['age'], xlabel='Age',
   ....:                      ylabel='SST (°C)')
   ....: 

In [11]: ax.grid()

In [12]: ax.legend()
Out[12]: <matplotlib.legend.Legend at 0x7f7242634550>

_images/predictplot_castaneda_pretty.png

In much the same way, we can view the distribution of the prediction prior and posterior with densityplot():

In [13]: ax = bsr.densityplot(prediction,
   ....:                      x=np.arange(1, 40, 0.1),
   ....:                      xlabel='SST (°C)')
   ....: 

In [14]: ax.grid(axis='x')

In [15]: ax.legend()
Out[15]: <matplotlib.legend.Legend at 0x7f724252de48>

_images/densityplot_castaneda_pretty.png

We can access the full prediction ensemble with:

In [16]: prediction.ensemble
Out[16]: 
array([[22.15034897, 19.93735788, 28.71824344, ..., 22.63209912,
        21.65818533, 22.98855481],
       [22.9722932 , 20.9287456 , 30.04947195, ..., 28.87735768,
        21.41092858, 23.03992783],
       [24.5673588 , 22.16406483, 31.4779479 , ..., 23.27718433,
        25.2469783 , 29.23222538],
       ...,
       [12.90251903, 15.76992949, 16.21296555, ..., 14.98565479,
        15.49510714, 14.08382262],
       [15.45817073, 17.24524088, 18.86186716, ..., 20.43865174,
        16.93475229, 14.28741749],
       [12.2875954 , 16.50004115, 16.58069012, ..., 20.16921891,
        20.21230458, 19.06554537]])

A useful summary is available using:

In [17]: prediction.percentile()[:10]  # Showing only the first few values.
Out[17]: 
array([[20.12036618, 24.52498535, 29.44182454],
       [20.07042247, 24.51395169, 29.43802039],
       [19.81838546, 23.89578648, 28.93577603],
       [20.06907822, 24.31681749, 29.52331259],
       [20.29399647, 24.61661976, 29.45651891],
       [20.140115  , 24.32071292, 29.38219064],
       [19.32560673, 23.58948729, 28.46959632],
       [20.6733432 , 24.86163052, 29.97853893],
       [19.40572198, 23.64455923, 28.52515245],
       [19.87691773, 24.19149754, 29.27673702]])

See the percentile() method for more options.

Analog prediction¶

Begin by loading data for a “Deep-Time” analog prediction:

In [18]: example_file = bsr.get_example_data('wilsonlake.csv')

In [19]: d = np.genfromtxt(example_file, delimiter=',', names=True)

This dataset is a TEX₈₆ record from record from Wilson Lake, New Jersey (Zachos et al. 2006). The file has two columns giving depth (m) and TEX₈₆:

In [20]: d['depth'][:5]
Out[20]: array([91.51, 92.95, 93.83, 95.96, 96.94])

In [21]: d['tex86'][:5]
Out[21]: array([0.79, 0.74, 0.77, 0.7 , 0.87])

We can run the analog prediction of SST with predict_seatemp_analog():

In [22]: search_tolerance = np.std(d['tex86'], ddof=1) * 2

In [23]: prediction = bsr.predict_seatemp_analog(d['tex86'], temptype='sst',
   ....:                                         prior_mean=30, prior_std=20,
   ....:                                         search_tol=search_tolerance,
   ....:                                         nens=500)
   ....: 

A Prediction instance is returned from the function. The arguments that we pass to the function indicate that the calibration is for SST. We also indicate a prior mean and standard deviation for the sea temperature inference. Note that we pass a search tolerance used to find analogous conditions across the global TEX₈₆ dataset. The nens argument is to reduce the size of the model parameter ensemble used for the inference - we’re using this option because otherwise this page of documentation would take too long to compile. This argument isn’t required if you’re following along with these examples on your own machine. By default, a progress bar is printed to the screen. This is an optional feature that can be switched off, for example, if you are processing many cores in batch.

Note

Analog predictions can be slow to calculate if many analog are selected. Be careful not to select too large a value for search_tol.

Specifically, for analog predictions, we can map the information used for the inference with analogmap():

In [24]: paleo_location = (38.2, -56.7)

In [25]: ax = bsr.analogmap(prediction, latlon=paleo_location)

In [26]: ax.legend()
Out[26]: <matplotlib.legend.Legend at 0x7f724247fcf8>

This maps the location of all available TEX₈₆ observations. We passed an optional paleo-location of the site to the mapping function. The map also indicates the grids that were identified as analogs for the prediction.

We can also examine our prediction as though it were a standard prediction. For example, we can plot a time series of the prediction:

In [27]: bsr.predictplot(prediction, x=d['depth'], xlabel='Depth (m)',
   ....:                 ylabel='SST (°C)')
   ....: 
Out[27]: <AxesSubplot:xlabel='Depth (m)', ylabel='SST (°C)'>

Forward prediction¶

For this example, we make inferences about TEX₈₆ from SST data using a forward-model prediction. We start by creating a SST series spanning from 0 - 40 °C:

In [28]: sst = np.arange(0, 41)

And now plug the SST data into predict_tex() along with additional information. In this case we’re using the same site location as in Standard prediction:

In [29]: prediction = bsr.predict_tex(sst, lon=34.0733, lat=31.6517, temptype='sst')

As might be expected, we can use the output of the forward prediction to parse and plot:

In [30]: ax = bsr.predictplot(prediction, x=sst,
   ....:                      xlabel='SST (°C)',
   ....:                      ylabel=r'TEX$_{86}$')
   ....: 

In [31]: ax.grid()

In [32]: ax.legend()
Out[32]: <matplotlib.legend.Legend at 0x7f724229e358>

Installation¶

Requirements¶

Python >= 3.4
numpy
scipy
matplotlib
cartopy
tqdm
attrs

Instructions¶

conda¶

To install baysparpy with conda, run:

$ conda install baysparpy -c sbmalev

This is usually the easiest way to install the package.

pip¶

To install with pip, run:

$ pip install baysparpy

and follow the on-screen prompts. We do not recommend installing with pip unless you are brave, or already have a copy of cartopy installed.

Testing¶

To run the testing suite after installation, first install py.test. Run:

$ pytest --pyargs bayspar

API reference¶

This page is summary of the baysparpy API.

Prediction¶

predict.Prediction(temptype, ensemble[, ...])

MCMC prediction

Creating a prediction¶

`predict_seatemp`(tex, lat, lon, prior_std, ...)	Predict sea temperature with TEX86
`predict_seatemp_analog`(tex, prior_std, ...)	Predict sea temperature with TEX86, using the analog method
`predict_tex`(seatemp, lat, lon, temptype[, nens])	Predict TEX86 from sea temperature

Plots & Graphics¶

`predictplot`(prediction[, ylabel, x, xlabel, ...])	Lineplot of prediction with uncertainty estimate
`densityplot`(prediction[, x, xlabel, ax])	Plot density of prediction prior and posterior
`analogmap`(prediction[, latlon, ax])	Map analog prediction with grids used for analog and TEX86 sites

Support and development¶

Please feel free to report bugs and issues or view the source code on GitHub.

License¶

baysparpy is available under the Open Source GPLv3.

The BAYSPAR TEX86 calibration, in Python¶

Documentation¶

Quick Overview¶

What’s New¶

v0.0.4¶

Enhancements¶

Bug fixes¶

v0.0.3¶

Enhancements¶

v0.0.2¶

Enhancements¶

Bug fixes¶

v0.0.1¶

Examples¶

Standard prediction¶

Analog prediction¶

Forward prediction¶

Installation¶

Requirements¶

Instructions¶

conda¶

pip¶

Testing¶

API reference¶

Prediction¶

Creating a prediction¶

Plots & Graphics¶

See also¶

Support and development¶

License¶

The BAYSPAR TEX₈₆ calibration, in Python¶