v0.0.5

README.rst

@@ -1,5 +1,5 @@
-baysparpy
-=========
+baysparpy3
+==========
 
 .. image:: https://travis-ci.org/brews/baysparpy.svg?branch=master
  :target: https://travis-ci.org/brews/baysparpy
@@ -9,10 +9,16 @@ An Open Source Python package for TEX86 calibration.
 
 This package is based on the original BAYSPAR (BAYesian SPAtially-varying Regression) for MATLAB (https://github.com/jesstierney/BAYSPAR).
 
+This package is the updated version of baysparpy with new module(s) such as `predict_tex_analog`
+based on the original python version of BAYSPAR (https://github.com/brews/baysparpy).
 
 Quick example
 -------------
 
+Standard prediction
+~~~~~~~~~~~~~~~~~~~
+
+
 First, load key packages and an example dataset:
 
 .. code-block:: python
@@ -37,6 +43,81 @@ To see actual numbers from the prediction, directly parse ``prediction.ensemble`
 
 You can also plot your prediction with ``bsr.predictplot()`` or ``bsr.densityplot()``.
 
+
+Analog prediction
+~~~~~~~~~~~~~~~~~
+
+
+Begin by loading data for a “Deep-Time” analog prediction:
+
+This dataset is a TEX86 record from record from Wilson Lake, New Jersey (`Zachos et al. 2006 <https://doi.org/10.1130/G22522.1>`_). The file has two columns giving depth (m) and TEX86:
+
+
+.. code-block:: python
+
+ example_file = bsr.get_example_data('wilsonlake.csv')
+ d = np.genfromtxt(example_file, delimiter=',', names=True)
+
+
+We can run the analog prediction of SST with ``predict_seatemp_analog()``.
+We can also examine our prediction as though it were a standard prediction. For example, we can plot a time series of the prediction:
+
+.. code-block:: python
+
+ search_tolerance = np.std(d['tex86'], ddof=1) * 2
+ prediction = bsr.predict_seatemp_analog(d['tex86'], temptype='sst',prior_mean=30, prior_std=20,search_tol=search_tolerance,nens=500)
+ ax = bsr.predictplot(prediction, x=d['depth'], xlabel='Depth (m)', ylabel='SST (°C)')
+ ax.grid()
+ ax.legend()
+
+
+Forward prediction
+~~~~~~~~~~~~~~~~~~
+
+
+For this example, we make inferences about TEX86 from SST data using a forward-model prediction. We start by creating a SST series spanning from 0 - 40 °C. 
+
+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:
+
+.. code-block:: python
+
+ sst = np.arange(0, 41)
+ 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:
+
+.. code-block:: python
+
+ ax = bsr.predictplot(prediction, x=sst,xlabel='SST (°C)',ylabel=r'TEX$_{86}$')
+ ax.grid()
+ ax.legend()
+
+
+Analog forward prediction
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+This tool will calculate forwarded TEX using given SST data. Here is an example:
+
+.. code-block:: python
+
+ sst = np.arange(0, 41)
+ prediction = bsr.predict_tex_analog(sst, temptype = 'sst', search_tol = 5., nens=8000)
+ ax = bsr.predictplot(prediction, x=sst,xlabel='SST (°C)',ylabel=r'TEX$_{86}$')
+ ax.grid()
+ ax.legend()
+
+
+First, we make inferences about deep-time TEX86 from SST data using a forward-model analog prediction. We start by creating a SST series spanning from 0 - 40 °C.
+
+And then plug the SST data into ``predict_tex_analog()`` along with additional information (search tolerance is 5 °C).
+
+We can use the output of the forward prediction to parse and plot.
+
+Read More
+~~~~~~~~~
+
+
 For further details, examples, and additional prediction functions, see the online documentation (https://baysparpy.readthedocs.io).
 
 
@@ -47,7 +128,7 @@ To install **baysparpy** with pip, run:
 
 .. code-block:: bash
 
- $ pip install baysparpy
+ $ pip install baysparpy3
 
 To install with conda, run:
 
@@ -62,7 +143,7 @@ Support and development
 
 - Documentation is available online (https://baysparpy.readthedocs.io).
 
-- Please feel free to report bugs and issues or view the source code on GitHub (https://github.com/brews/baysparpy).
+- Please feel free to report bugs and issues or view the source code on GitHub (https://github.com/mingsongli/baysparpy).
 
 
 License

bayspar/__init__.py

@@ -1,3 +1,4 @@
-from bayspar.predict import predict_tex, predict_seatemp, predict_seatemp_analog
+name = 'baysparpy3'
+from bayspar.predict import predict_tex, predict_tex_analog, predict_seatemp, predict_seatemp_analog
 from bayspar.utils import get_example_data
 from bayspar.plot import predictplot, analogmap, densityplot

bayspar/modelparams/__init__.py

@@ -1 +1 @@
-from bayspar.modelparams.core import get_draws
+from bayspar.modelparams.core import get_draws, get_draws_analog

bayspar/modelparams/core.py

@@ -1,3 +1,17 @@
+"""
+Core
+ Originator: Steven Brewster Malevich
+ University of Arizona Department of Geosciences
+
+ Revisions: Mingsong Li
+ Penn State Geosciences
+ Date: Sept 23, 2019
+
+ Purpose: add get_draws_analog module
+ Get Draws instance for a draw type
+ for the analog model
+
+"""
 import os.path
 from copy import deepcopy
 from pkgutil import get_data
@@ -19,24 +33,25 @@ def get_matlab_resource(resource, package='bayspar', **kwargs):
 
 def read_draws(flstr, drawtype):
  """Grab single squeezed array from package resources
+ example:
+ alpha_samples_comp=read_draws('alpha_samples_comp.mat', 'sst')
  """
  drawtype = drawtype.lower()
-
- var_template = 'modelparams/Output_SpatAg_{0}/{1}'
-
  varstr = TRANSLATE_VAR[drawtype]
- varstr_full = os.path.splitext(flstr)[0]
+
+ var_template = 'modelparams/Output_SpatAg_{0}/{1}'
+ # resource_str = 'modelparams/Output_SpatAg_{0}/{1}'.format('SST', 'alpha_samples_comp.mat')
  resource_str = var_template.format(varstr, flstr)
+ varstr_full = os.path.splitext(flstr)[0]
  var = get_matlab_resource(resource_str, squeeze_me=True)
  return var[varstr_full]
 
-
 @attr.s
 class Draws:
  """Spatially-aware modelparams draws
  """
- alpha_samples_comp = attr.ib()
- beta_samples_comp = attr.ib()
+ alpha_samples_comp = attr.ib() # true data may be either alpha_samples_comp.mat or alpha_samples.mat
+ beta_samples_comp = attr.ib() # true data may be either beta_samples_comp.mat or beta_samples.mat
  tau2_samples = attr.ib()
  locs_comp = attr.ib()
  _half_grid_space = attr.ib(default=10)
@@ -66,7 +81,7 @@ def find_alphabeta_near(self, lat, lon):
  alpha_select = self.alpha_samples_comp[idx].squeeze()
  beta_select = self.beta_samples_comp[idx].squeeze()
  return alpha_select, beta_select
-
+ 
 
 class BadLatlonError(Exception):
  """Raised when latitude or longitude is outside (-90, 90) or (-180, 180)
@@ -79,7 +94,7 @@ class BadLatlonError(Exception):
  def __init__(self, latlon):
  self.latlon = latlon
 
-
+# for modern tex forward
 draws_sst = Draws(alpha_samples_comp=read_draws('alpha_samples_comp.mat', 'sst'),
  beta_samples_comp=read_draws('beta_samples_comp.mat', 'sst'),
  tau2_samples=read_draws('tau2_samples.mat', 'sst'),
@@ -91,11 +106,30 @@ def __init__(self, latlon):
  tau2_samples=read_draws('tau2_samples.mat', 'subt'),
  locs_comp=read_draws('Locs_Comp.mat', 'subt'))
 
-
 def get_draws(drawtype):
  """Get Draws instance for a draw type
  """
  if drawtype == 'sst':
  return deepcopy(draws_sst)
  elif drawtype == 'subt':
  return deepcopy(draws_subt)
+
+# for analog tex forward model
+draws_sst_analog = Draws(alpha_samples_comp=read_draws('alpha_samples.mat', 'sst'),
+ beta_samples_comp=read_draws('beta_samples.mat', 'sst'),
+ tau2_samples=read_draws('tau2_samples.mat', 'sst'),
+ locs_comp=read_draws('Locs_Comp.mat', 'sst'))
+
+draws_subt_analog = Draws(alpha_samples_comp=read_draws('alpha_samples.mat', 'subt'),
+ beta_samples_comp=read_draws('beta_samples.mat', 'subt'),
+ tau2_samples=read_draws('tau2_samples.mat', 'subt'),
+ locs_comp=read_draws('Locs_Comp.mat', 'subt'))
+
+def get_draws_analog(drawtype):
+ """Get Draws instance for a draw type
+ for the analog model"""
+ if drawtype == 'sst':
+ return deepcopy(draws_sst_analog)
+ elif drawtype == 'subt':
+ return deepcopy(draws_subt_analog)
+

bayspar/observations/core.py

@@ -1,3 +1,17 @@
+"""
+Core
+ Originator: Steven Brewster Malevich
+ University of Arizona Department of Geosciences
+
+ Revisions: Mingsong Li
+ Penn State Geosciences
+ Date: Sept 23, 2019
+
+ Purpose: add find_t_within_tolerance module in TexObs class
+ : Find mean modern sea temperature observations
+ that are within ± tolerance
+
+"""
 from copy import deepcopy
 from pkgutil import get_data
 from io import BytesIO
@@ -47,7 +61,10 @@ def read_tex(obstype):
  obs_stack = var['Obs_Stack'].squeeze().item()
  inds_stack = var['Inds_Stack'].squeeze().item()
  obslocs = np.array([x[...].squeeze() for x in var['ObsLocs'].item().ravel() if x.any()])
- return locs, obs_stack, inds_stack, obslocs
+ # return locs, obs_stack, inds_stack, obslocs
+ # M.Li
+ target_stack = var['Target_Stack'].squeeze().item()
+ return locs, obs_stack, inds_stack, obslocs, target_stack
 
 
 def chord_distance(latlon1, latlon2):
@@ -177,7 +194,55 @@ class TexObs:
  obs_stack = attr.ib()
  inds_stack = attr.ib()
  obslocs = attr.ib()
-
+ # M.Li added
+ target_stack = attr.ib()
+
+ #M.Li added
+ def find_t_within_tolerance(self,t,tolerance):
+ """Find mean modern sea temperature observations that are within ± tolerance
+ cycle through the alpha/beta grid cells
+
+ Parameters
+ ----------
+ t ： float
+ Mean sea temperature value.
+ tolerance : float
+ Value added and subtracted from 't' to get upper and lower tolerance bounds.
+
+ Returns
+ ----------
+ latlon_match : list
+ An n-length list of latlon tuples where matches were found.
+ vals_match : ndarray
+ A 1d array (n) of corresponding TEX86 averages from each match.
+ inder_g : ndarray
+ A 1d array (n) of location in the big grids that has sea temp within
+ tolerance
+
+ Originator: Mingsong Li, Penn State Geosciences
+ Date: 25 September, 2019
+ """
+ n_bg = len(self.locs) # number of big grids
+ upper_bound = t + tolerance # upper bound of sea temperature range
+ lower_bound = t - tolerance # lower bound of sea temperature range
+ inder_g =[]
+ vals_match = []
+ # search eachi big grid
+ for kk in range(1, n_bg+1):
+ # find the mean SST obs corresponding to this index location
+ # 
+ vals = self.target_stack[self.inds_stack == kk]
+ #
+ vals_mean = vals.mean()
+ # if the mean of vals in the big grids within tolerance, add it to inder_g
+ if lower_bound <= vals_mean <= upper_bound:
+ idx = kk-1 # Back to 0-based idx
+ inder_g.append(idx)
+ vals_match.append(vals_mean)
+
+ latlon_match = [tuple(x) for x in self.locs[inder_g, ::-1].tolist()]
+ return latlon_match, np.array(vals_match), np.array(inder_g)
+
  def find_within_tolerance(self, x, tolerance):
  """Find mean TEX86 observations that are within ± tolerance from x