Pyleoclim User API

Pyleoclim, like a lot of other Python packages, follows an object-oriented design. It sounds fancy, but it really is quite simple. What this means for you is that we’ve gone through the trouble of coding up a lot of timeseries analysis methods that apply in various situations - so you don’t have to worry about that. These situations are described in classes, the beauty of which is called “inheritance” (see link above). Basically, it allows to define methods that will automatically apply to your dataset, as long as you put your data within one of those classes. A major advantage of object-oriented design is that you, the user, can harness the power of Pyleoclim methods in very few lines of code through the user API without ever having to get your hands dirty with our code (unless you want to, of course). The flipside is that any user would do well to understand Pyleoclim classes, what they are intended for, and what methods they support.

The following describes the various classes that undergird the Pyleoclim edifice.

Series (pyleoclim.Series)

class pyleoclim.core.series.Series(time, value, time_unit=None, time_name=None, value_name=None, value_unit=None, label=None, importedFrom=None, archiveType=None, control_archiveType=False, log=None, keep_log=False, sort_ts='ascending', dropna=True, verbose=True, clean_ts=False, auto_time_params=None)[source]

The Series class describes the most basic objects in Pyleoclim. A Series is a simple dictionary that contains 3 things:

value, an array of real-valued numbers;
time, a coordinate axis at which those values were obtained ;
optionally, some metadata about both axes, like units, labels and origin.

How to create, manipulate and use such objects is described in PyleoTutorials.

Parameters:

time (list or numpy.array) – time axis (prograde or retrograde)
value (list of numpy.array) – values of the dependent variable (y)
time_unit (string) – Units for the time vector (e.g., ‘ky BP’). Default is None. in which case ‘years CE’ are assumed
time_name (string) – Name of the time vector (e.g., ‘Time’,’Age’). Default is None. This is used to label the time axis on plots
value_name (string) – Name of the value vector (e.g., ‘temperature’) Default is None
value_unit (string) – Units for the value vector (e.g., ‘deg C’) Default is None
label (string) – Name of the time series (e.g., ‘NINO 3.4’) Default is None
log (dict) – Dictionary of tuples documentating the various transformations applied to the object
keep_log (bool) – Whether to keep a log of applied transformations. False by default
importedFrom (string) – source of the dataset. If it came from a LiPD file, this could be the datasetID property
archiveType (string) – climate archive, one of ‘Borehole’, ‘Coral’, ‘FluvialSediment’, ‘GlacierIce’, ‘GroundIce’, ‘LakeSediment’, ‘MarineSediment’, ‘Midden’, ‘MolluskShell’, ‘Peat’, ‘Sclerosponge’, ‘Shoreline’, ‘Speleothem’, ‘TerrestrialSediment’, ‘Wood’ Reference: https://lipdverse.org/vocabulary/archivetype/
control_archiveType ([True, False]) – Whether to standardize the name of the archiveType agains the vocabulary from: https://lipdverse.org/vocabulary/paleodata_proxy/. If set to True, will only allow for these terms and automatically convert known synonyms to the standardized name. Only standardized variable names will be automatically assigned a color scheme. Default is False.
dropna (bool) – Whether to drop NaNs from the series to prevent downstream functions from choking on them defaults to True
sort_ts (str) – Direction of sorting over the time coordinate; ‘ascending’ or ‘descending’ Defaults to ‘ascending’
verbose (bool) – If True, will print warning messages if there are any
clean_ts (boolean flag) – set to True to remove the NaNs and make time axis strictly prograde with duplicated timestamps reduced by averaging the values Default is None (marked for deprecation)
auto_time_params (bool) – If True, uses tsbase.disambiguate_time_metadata to ensure that time_name and time_unit are usable by Pyleoclim. This may override the provided metadata. If False, the provided time_name and time_unit are used. This may break some functionalities (e.g. common_time and convert_time_unit), so use at your own risk. If not provided, code will set to True for internal consistency.

Examples

Import the Southern Oscillation Index (SOI) and display a quick synopsis:

import pyleoclim as pyleo
soi = pyleo.utils.load_dataset('SOI')
soi.view()

	SOI [mb]
Time [year C.E.]
1866.002894	-0.62
1866.087769	-0.12
1866.164431	-0.62
1866.249306	-0.65
1866.331443	0.04
...	...
2024.752349	0.35
2024.837224	0.61
2024.919362	1.09
2025.004237	0.26
2025.089112	0.21

1910 rows × 1 columns

annualize(months=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12], min_res=0.25, frac_req_months=0.6666666666666666)[source]

Annualize subannual data by averaging values within specified months for each year.

This method converts subannual time series data to annual resolution by computing weighted averages of values from specified months using the custom_year_averages utility function.

Parameters:

months (list of int, optional) – List of months to include in the annual aggregation. Must be consecutive months either within a calendar year or spanning across years (e.g., [10, 11, 12, 1, 2, 3] for Oct-Mar). Default is a calendar average, [1 … 12].
min_res (float, optional) – Minimum resolution threshold (in years) to determine if data is subannual. Default is 0.25 (quarterly or finer resolution required).
frac_req_months (float, optional) – Minimum fraction of requested months that must be present for a year to be included. Default is 2/3 (67%). For example, if requesting 12 months, at least 8 months must be present. If requesting 3 months (e.g., DJF), at least 2 months must be present. Years that do not meet this threshold are dropped from the resulting series. Set to 1.0 to require all months, or lower values (e.g., 1/3) for more lenient inclusion.

Returns:

A new Series object with annualized data.

Return type:

Series

Examples

Calendar average

soi = pyleo.utils.load_dataset('SOI')
dt = soi.resolution().describe()['median']
print(f"The series' resolution is {dt*12:.0f} month")

soi_a = soi.annualize()
dta = soi_a.resolution().describe()['median']
print(f"The series' resolution is now {dta:.0f} year")
fig, ax = soi.plot(title='Jan-Dec averaging')
soi_a.plot(marker='o',ax=ax)

The series' resolution is 1 month

The series' resolution is now 1 year

<Axes: title={'center': 'Jan-Dec averaging'}, xlabel='Time [year C.E.]', ylabel='SOI [mb]'>

JJA average (straightforward)

soi_jja = soi.annualize(months=[6, 7 , 8])
fig, ax = soi.plot(title='JJA averaging')
soi_jja.plot(marker='o',ax=ax, label='JJA average')

<Axes: title={'center': 'JJA averaging'}, xlabel='Time [year C.E.]', ylabel='SOI [mb]'>

DJF average : straddles a year; handles it gracefully

soi_djf = soi.annualize(months=[12, 1 , 2])
fig, ax = soi.plot(title='DJF averaging')
soi_djf.plot(marker='o',ax=ax, label='DJF average')

<Axes: title={'center': 'DJF averaging'}, xlabel='Time [year C.E.]', ylabel='SOI [mb]'>

Varying the fraction of required months

AprMar = [4,5,6,7,8,9,10,11,12,1,2,3]
soi_am_default = soi.annualize(months=AprMar)
soi_am_stringent = soi.annualize(months=AprMar,frac_req_months=1.0)

fig, ax = soi.plot(title='Apr-Mar averaging', xlim = [2000, 2026])
soi_am_default.plot(marker='o',ax=ax, label='Apr-Mar, $f=2/3$')
soi_am_stringent.plot(marker='o',ax=ax, label='Apr-Mar, $f=1.0$')

<Axes: title={'center': 'Apr-Mar averaging'}, xlabel='Time [year C.E.]', ylabel='SOI [mb]'>

The last year is incomplete, so insisting on complete coverage results in dropping it

Notes

This method requires subannual data (median resolution <= min_res)
Months must be consecutive to define a clear averaging period
Uses weighted averaging to account for potentially uneven temporal spacing

bin(keep_log=False, **kwargs)[source]

Bin values in a time series

Parameters:

keep_log (Boolean) – if True, adds this step and its parameters to the series log.
kwargs – Arguments for binning function. See pyleoclim.utils.tsutils.bin for details

Returns:

new – An binned Series object

Return type:

Series

See also

pyleoclim.utils.tsutils.bin: bin the series values into evenly-spaced time bins

causality(target_series, method='liang', timespan=None, settings=None, common_time_kwargs=None)[source]

Perform causality analysis with the target timeseries. Specifically, whether there is information in the target series that influenced the original series.: If the two series have different time axes, they are first placed on a common timescale (in ascending order).

Parameters:

target_series (Series) – A pyleoclim Series object on which to compute causality
method ({'liang', 'granger'}) – The causality method to use.
timespan (tuple) – The time interval over which to perform the calculation
settings (dict) – Parameters associated with the causality methods. Note that each method has different parameters. See individual methods for details
common_time_kwargs (dict) – Parameters for the method MultipleSeries.common_time(). Will use interpolation by default.

Returns:

res – Dictionary containing the results of the the causality analysis. See indivudal methods for details

Return type:

dict

See also

pyleoclim.utils.causality.liang_causality: Liang causality
pyleoclim.utils.causality.granger_causality: Granger causality

Examples

Liang causality

ts_nino=pyleo.utils.load_dataset('NINO3')
ts_air=pyleo.utils.load_dataset('AIR')

We use the specific params below to lighten computations; you may drop settings for real work

liang_N2A = ts_air.causality(ts_nino, settings={'nsim': 20, 'signif_test': 'isopersist'})
print(liang_N2A)
liang_A2N = ts_nino.causality(ts_air, settings={'nsim': 20, 'signif_test': 'isopersist'})
print(liang_A2N)

liang_N2A['T21']/liang_A2N['T21']

{'T21': np.float64(0.01644548028647295), 'tau21': np.float64(0.011968992003953967), 'Z': np.float64(1.3740071244963796), 'dH1_star': np.float64(-0.6359251528278479), 'dH1_noise': np.float64(0.3521058551681981), 'signif_qs': [0.005, 0.025, 0.05, 0.95, 0.975, 0.995], 'T21_noise': array([-4.32589932e-05, -4.32589932e-05, -3.23075012e-05,  2.96208614e-03,
        3.24962027e-03,  3.24962027e-03]), 'tau21_noise': array([-3.19603765e-05, -3.19603765e-05, -2.40302600e-05,  2.22858615e-03,
        2.50778462e-03,  2.50778462e-03])}

{'T21': np.float64(0.005840218794917535), 'tau21': np.float64(0.0473182615992069), 'Z': np.float64(0.12342420447279116), 'dH1_star': np.float64(-0.5094709112672597), 'dH1_noise': np.float64(0.44321082713353344), 'signif_qs': [0.005, 0.025, 0.05, 0.95, 0.975, 0.995], 'T21_noise': array([-0.00020216, -0.00020216, -0.00015725,  0.00160315,  0.00231519,
        0.00231519]), 'tau21_noise': array([-0.00158658, -0.00158658, -0.00125493,  0.01069281,  0.01450644,
        0.01450644])}

np.float64(2.815901400951737)

Both information flows (T21) are positive, but the flow from NINO3 to AIR is about 3x as large as the other way around, suggesting that NINO3 influences AIR much more than the other way around, which conforms to physical intuition.

To implement Granger causality, simply specfiy the method:

granger_A2N = ts_nino.causality(ts_air, method='granger')
granger_N2A = ts_air.causality(ts_nino, method='granger')

Granger Causality
number of lags (no zero) 1
ssr based F test:         F=20.8492 , p=0.0000  , df_denom=1592, df_num=1
ssr based chi2 test:   chi2=20.8885 , p=0.0000  , df=1
likelihood ratio test: chi2=20.7529 , p=0.0000  , df=1
parameter F test:         F=20.8492 , p=0.0000  , df_denom=1592, df_num=1

Granger Causality
number of lags (no zero) 1
ssr based F test:         F=18.6927 , p=0.0000  , df_denom=1592, df_num=1
ssr based chi2 test:   chi2=18.7280 , p=0.0000  , df=1
likelihood ratio test: chi2=18.6189 , p=0.0000  , df=1
parameter F test:         F=18.6927 , p=0.0000  , df_denom=1592, df_num=1

Note that the output is fundamentally different for the two methods. Granger causality cannot discriminate between NINO3 -> AIR or AIR -> NINO3, in this case. This is not unusual, and one reason why it is no longer in wide use.

center(timespan=None, keep_log=False)[source]

Centers the series (i.e. renove its estimated mean)

Parameters:

timespan (tuple or list) – The timespan over which the mean must be estimated. In the form [a, b], where a, b are two points along the series’ time axis.
keep_log (Boolean) – if True, adds the previous mean and method parameters to the series log.

Returns:

new – The centered series object

Return type:

Series

clean(verbose=False, keep_log=False)[source]

Clean up the timeseries by removing NaNs and sort with increasing time points

Parameters:

verbose (bool) – If True, will print warning messages if there are any
keep_log (Boolean) – if True, adds this step and its parameters to the series log.

Returns:

new – Series object with removed NaNs and sorting

Return type:

Series

compare(ts, row_clr='palegoldenrod')[source]

Graphically compare two Series objects

Parameters:

ts (pyleo.Series) – The target series

Returns:

fig (Matplotlib figure)
df (pandas dataframe) – dataframe containing metadata ; identical rows are highlighted in Jupyter Notebook

Examples

Compare the SOI and NINO3 series

GISP2 = pyleo.utils.load_dataset('GISP2')
EDC_dD = pyleo.utils.datasets.load_dataset('EDC-dD')

fig, df = GISP2.compare(EDC_dD)
df  # this will display the metadata as a dataframe in a notebook

	GISP2	EPICA Dome C dD
archiveType	GlacierIce	GlacierIce
control_archiveType	nan	False
elevation	nan	3233
importedFrom	https://www.ncei.noaa.gov/access/paleo-search/study/17796	https://www.ncei.noaa.gov/pub/data/paleo/icecore/antarctica/epica_domec/edc3deuttemp2007.txt
lat	nan	-75.101100
log	nan	None
lon	nan	123.347800
observationType	nan	hydrogen isotopes
sensorType	nan	ice sheet
time_name	Age	Age
time_unit	yr BP	y BP
value_name	$\delta^{18} \mathrm{O}$	$\delta \mathrm{D}$
value_unit	‰	‰

convert_time_unit(time_unit='ky BP', keep_log=False)[source]

Convert the time units of the Series object

Parameters:

time_unit (str) –
the target time unit, possible input: {

’year’, ‘years’, ‘yr’, ‘yrs’, ‘CE’, ‘AD’, ‘y BP’, ‘yr BP’, ‘yrs BP’, ‘year BP’, ‘years BP’, ‘ky BP’, ‘kyr BP’, ‘kyrs BP’, ‘ka BP’, ‘ka’, ‘my BP’, ‘myr BP’, ‘myrs BP’, ‘ma BP’, ‘ma’,

}
keep_log (Boolean) – if True, adds this step and its parameter to the series log.

Examples

ts = pyleo.utils.load_dataset('SOI')
tsBP = ts.convert_time_unit(time_unit='yrs BP')
print('Original timeseries:')
print('time unit:', ts.time_unit)
print('time:', ts.time[:10])
print()
print('Converted timeseries:')
print('time unit:', tsBP.time_unit)
print('time:', tsBP.time[:10])

Original timeseries:
time unit: year C.E.
time: [1866.00289418 1866.08776937 1866.16443083 1866.24930601 1866.33144329
 1866.41631848 1866.49845576 1866.58333094 1866.66820613 1866.75034341]

Converted timeseries:
time unit: yrs BP
time: [-75.08716159 -75.0022864  -74.91741122 -74.83527394 -74.75039875
 -74.66826147 -74.58338629 -74.4985111  -74.41637382 -74.33149863]

copy()[source]

Make a copy of the Series object

Returns:: Series – A copy of the Series object
Return type:: Series

correlation(target_series, alpha=0.05, statistic='pearsonr', method='phaseran', number=1000, timespan=None, settings=None, seed=None, common_time_kwargs=None, mute_pbar=False)[source]

Estimates the correlation and its associated significance between two time series (not necessarily IID) as per [1]

The significance of the correlation is assessed using one of the following methods:

‘ttest’: T-test adjusted for effective sample size, see [1]
‘ar1sim’: AR(1) modeling of x and y (Monte Carlo method)
‘CN’: colored noise (power-law spectrum) modeling of x and y (Monte Carlo method)
‘phaseran’: phase randomization of original inputs. (Monte Carlo method, default), see [2]
‘built-in’: uses built-in method from scipy (function of the statistic used)

Note: Up to version v0.14.0. ar1sim was called “isopersistent”, phaseran was called “isospectral”

The T-test is a parametric test, hence computationally cheap, but can only be performed in ideal circumstances (Gaussian, identically distributed data and moderate autocorrelation). The others are non-parametric, but their computational requirements scale with the number of simulations.

The choice of significance test and associated number of Monte-Carlo simulations are passed through the settings parameter.

Parameters:

target_series (Series) – A pyleoclim Series object
alpha (float) – The significance level (default: 0.05)
statistic (str) –

statistic being evaluated. Can use any of the SciPy-supported ones:
https://docs.scipy.org/doc/scipy/reference/stats.html#association-correlation-tests Currently supported: [‘pearsonr’,’spearmanr’,’pointbiserialr’,’kendalltau’,’weightedtau’]

Default: ‘pearsonr’.
method (str, {'ttest','built-in','ar1sim','phaseran','CN'}) – method for significance testing. Default is ‘phaseran’ - ‘ttest’ implements the T-test with degrees of freedom adjusted for autocorrelation, as done in [1] - ‘built-in’ uses the p-value that ships with the SciPy function. - ‘ar1sim’ (formerly ‘isopersistent’) tests against an ensemble of AR(1) series fitted to the originals - ‘CN’ tests against an ensemble of colored noise series (power-law spectra) fitted to the originals - ‘phaseran’ (formerly ‘isospectral’) tests against phase-randomized surrogates (aka the method of Ebisuzaki [2]) The old options ‘isopersistent’ and ‘isospectral’ still work, but trigger a deprecation warning. Note that ‘weightedtau’ does not have a known distribution, so the ‘built-in’ method returns an error in that case.
number (int) – the number of surrogates series used for simulation, as appropriate for the method. Ignored if the method is ‘ttest’ or ‘built-in’
timespan (tuple) – The time interval over which to perform the calculation
settings (dict) – optional parameters for the measures of association. See https://docs.scipy.org/doc/scipy/reference/stats.html#association-correlation-tests
common_time_kwargs (dict) – Parameters for the method MultipleSeries.common_time(). Will use interpolation by default.
seed (float or int) – random seed for isopersistent and isospectral methods
mute_pbar (bool, optional) – Mute the progress bar. The default is False.

Returns:

corr – the result object, containing

rfloat
association metric (typically, correlation coefficient)
r_critfloat
critical value of the statistic (unless ‘method’ = ‘built-in’, in which case this quantity is not computed, and output as np.nan)
pfloat
the p-value
signifbool
true if significant; false otherwise Note that signif = True if and only if p <= alpha.
alphafloat
the significance level

Return type:

pyleoclim.Corr

See also

pyleoclim.utils.correlation.association: SciPy measures of association between variables
pyleoclim.series.surrogates: parametric and non-parametric surrogates of any Series object
pyleoclim.multipleseries.common_time: Aligning time axes

References

_ [1] Hu, J., J. Emile-Geay, and J. Partin (2017), Correlation-based interpretations of paleoclimate data – where statistics meet past climates, Earth and Planetary Science Letters, 459, 362–371, doi:10.1016/j.epsl.2016.11.048.

_ [2] Ebisuzaki, W. (1997), A method to estimate the statistical significance of a correla- tion when the data are serially correlated, Journal of Climate, 10(9), 2147–2153, doi:10.1175/1520- 0442(1997)010¡2147:AMTETS¿2.0.CO;2.

Examples

Let us compute the correlation between the Nino3.4 index and the Deasonalized All Indian Rainfall Index. For expediency, we limit the Monte Carlo tests to 20 surrogates, which is not sufficient to obtain accurate p-values. The default number, 1000, is more respectable, though to avoid p-hacking, we recommend bumping it even higher if at all possible.

ts_air = pyleo.utils.load_dataset('AIR')
ts_nino = pyleo.utils.load_dataset('NINO3')

# with  20 surrogates and the default significance method, 'phaseran'`
corr_res = ts_nino.correlation(ts_air, number=20)
print(corr_res)

# using a simple t-test with DOFs adjusted for autocorrelation (only with Person's r')
corr_res = ts_nino.correlation(ts_air, method='ttest')
print(corr_res)

# using  "isopersistent" surrogates (AR(1) simulation)
# and setting an arbitrary random seed to fix the result
corr_res = ts_nino.correlation(ts_air, method = 'ar1sim', number=20, seed=2333)
print(corr_res)

# changing the statistic works with all surrogate-based significance methods (but not the 'ttest' method, which only applies to a 'pearsonr' statistic)
corr_res = ts_nino.correlation(ts_air, statistic='spearmanr', method = 'phaseran', number=20)
print(corr_res)

# To use the built-in signficance test:
corr_res2s = ts_nino.correlation(ts_air, statistic='spearmanr', method = 'built-in')
print(corr_res2s)

# To modify the method, use `settings`. For instance, to specify a 1-sided test instead of a 2-sided one for Spearman's R:
corr_res1s = ts_nino.correlation(ts_air, statistic='spearmanr', method = 'built-in', settings={'alternative':'less'})
print(corr_res1s.p-corr_res2s.p) # shows the difference in p-values

  statistic    critical value  p-value    signif. (α: 0.05)
-----------  ----------------  ---------  -------------------
  -0.152394        -0.0478127  < 1e-6     True

  statistic    critical value  p-value    signif. (α: 0.05)
-----------  ----------------  ---------  -------------------
  -0.152394        -0.0557396  < 1e-2     True

  statistic    critical value  p-value    signif. (α: 0.05)
-----------  ----------------  ---------  -------------------
  -0.152394        -0.0722729  < 1e-6     True

  statistic    critical value  p-value    signif. (α: 0.05)
-----------  ----------------  ---------  -------------------
 -0.0931086        -0.0487597  < 1e-6     True

  statistic    critical value  p-value    signif. (α: 0.05)
-----------  ----------------  ---------  -------------------
 -0.0931086               nan  < 1e-3     True

-9.771046201450402e-05

property datetime_index

Convert time to pandas DatetimeIndex.

Note: conversion will happen using time_unit, and will assume:

the number of seconds per year is calculated using UDUNITS, see http://cfconventions.org/cf-conventions/cf-conventions#time-coordinate
time refers to the Gregorian calendar. If using a different calendar, then please make sure to do any conversions before hand.

detrend(method='emd', keep_log=False, preserve_mean=False, **kwargs)[source]

Detrend Series object

Parameters:

method (str, optional) –
The method for detrending. The default is ‘emd’. Options include:
- ”linear”: the result of a n ordinary least-squares stright line fit to y is subtracted.
- ”constant”: only the mean of data is subtracted.
- ”savitzky-golay”, y is filtered using the Savitzky-Golay filters and the resulting filtered series is subtracted from y.
- ”emd” (default): Empirical mode decomposition. The last mode is assumed to be the trend and removed from the series
keep_log (boolean) – if True, adds the removed trend and method parameters to the series log.
preserve_mean (boolean) – if True, ensures that the mean of the series is preserved despite the detrending
kwargs (dict) – Relevant arguments for each of the methods.

Returns:

new – Detrended Series object in “value”, with new field “trend” added

Return type:

Series

See also

pyleoclim.utils.tsutils.detrend: detrending wrapper functions

Examples

lr04 = pyleo.utils.load_dataset('LR04')
fig, ax = lr04.plot(invert_yaxis=True)
ts_emd = lr04.detrend(method='emd',preserve_mean=True)
ts_emd.plot(label=lr04.label+', EMD detrend',ax=ax)

<Axes: xlabel='Age [ky BP]', ylabel='$\\delta^{18} \\mathrm{O}$ [‰]'>

equals(ts, index_tol=5, value_tol=1e-05)[source]

Test whether two objects contain the same elements (values and datetime_index) A printout is returned if metadata are different, but the statement is considered True as long as data match.

Parameters:

ts (Series object) – The target series for the comparison
index_tol (int, default 5) – tolerance on difference in datetime indices (in dtype units, which are seconds by default)
value_tol (float, default 1e-5) – tolerance on difference in values (in %)

Returns:

same_data (bool) – Truth value of the proposition “the two series have the same data”.
same_metadata (bool) – Truth value of the proposition “the two series have the same metadata”.

Examples

soi = pyleo.utils.load_dataset('SOI')
NINO3 = pyleo.utils.load_dataset('NINO3')
soi.equals(NINO3)

The two series have different lengths, left: 1910 vs right: 1596
Metadata are different:
value_unit property -- left: mb, right: $^{\circ}$C
value_name property -- left: SOI, right: NINO3
label property -- left: Southern Oscillation Index, right: NINO3 SST

(False, False)

fill_na(timespan=None, dt=1, keep_log=False)[source]

Fill NaNs into the timespan

Parameters:

timespan (tuple or list) – The list of time points for slicing, whose length must be 2. For example, if timespan = [a, b], then the sliced output includes one segment [a, b]. If None, will use the start point and end point of the original timeseries
dt (float) – The time spacing to fill the NaNs; default is 1.
keep_log (Boolean) – if True, adds this step and its parameters to the series log.

Returns:

new – The sliced Series object.

Return type:

Series

filter(cutoff_freq=None, cutoff_scale=None, method='butterworth', keep_log=False, **kwargs)[source]

Filtering methods for Series objects using four possible methods:

By default, this method implements a lowpass filter, though it can easily be turned into a bandpass or high-pass filter (see examples below).

Parameters:

method (str, {'savitzky-golay', 'butterworth', 'firwin', 'lanczos'}) – the filtering method: - ‘butterworth’: a Butterworth filter (default = 3rd order) - ‘savitzky-golay’: Savitzky-Golay filter - ‘firwin’: finite impulse response filter design using the window method, with default window as Hamming - ‘lanczos’: Lanczos zero-phase filter
cutoff_freq (float or list) – The cutoff frequency only works with the Butterworth method. If a float, it is interpreted as a low-frequency cutoff (lowpass). If a list, it is interpreted as a frequency band (f1, f2), with f1 < f2 (bandpass). Note that only the Butterworth option (default) currently supports bandpass filtering.
cutoff_scale (float or list) – cutoff_freq = 1 / cutoff_scale The cutoff scale only works with the Butterworth method and when cutoff_freq is None. If a float, it is interpreted as a low-frequency (high-scale) cutoff (lowpass). If a list, it is interpreted as a frequency band (f1, f2), with f1 < f2 (bandpass).
keep_log (Boolean) – if True, adds this step and its parameters to the series log.
kwargs (dict) – a dictionary of the keyword arguments for the filtering method, see pyleoclim.utils.filter.savitzky_golay, pyleoclim.utils.filter.butterworth, pyleoclim.utils.filter.lanczos and pyleoclim.utils.filter.firwin for the details

Returns:

new

Return type:

Series

See also

pyleoclim.utils.filter.butterworth: Butterworth method
pyleoclim.utils.filter.savitzky_golay: Savitzky-Golay method
pyleoclim.utils.filter.firwin: FIR filter design using the window method
pyleoclim.utils.filter.lanczos: lowpass filter via Lanczos resampling

Examples

In the example below, we generate a signal as the sum of two signals with frequency 10 Hz and 20 Hz, respectively. Then we apply a low-pass filter with a cutoff frequency at 15 Hz, and compare the output to the signal of 10 Hz. After that, we apply a band-pass filter with the band 15-25 Hz, and compare the outcome to the signal of 20 Hz.

Generating the test data

import numpy as np
t = np.linspace(0, 1, 1000)
sig1 = np.sin(2*np.pi*10*t)
sig2 = np.sin(2*np.pi*20*t)
sig = sig1 + sig2
ts1 = pyleo.Series(time=t, value=sig1)
ts2 = pyleo.Series(time=t, value=sig2)
ts = pyleo.Series(time=t, value=sig)

fig, ax = ts.plot(label='mix')
ts1.plot(ax=ax, label='10 Hz')
ts2.plot(ax=ax, label='20 Hz')
ax.legend(loc='upper left', bbox_to_anchor=(0, 1.1), ncol=3)

Time axis values sorted in ascending order
Time axis values sorted in ascending order
Time axis values sorted in ascending order

<matplotlib.legend.Legend at 0x7150390427e0>

Applying a low-pass filter

fig, ax = ts.plot(label='mix')
ts.filter(cutoff_freq=15).plot(ax=ax, label='After 15 Hz low-pass filter')
ts1.plot(ax=ax, label='10 Hz')
ax.legend(loc='upper left', bbox_to_anchor=(0, 1.1), ncol=3)

<matplotlib.legend.Legend at 0x7150390425a0>

Applying a band-pass filter

fig, ax = ts.plot(label='mix')
ts.filter(cutoff_freq=[15, 25]).plot(ax=ax, label='After 15-25 Hz band-pass filter')
ts2.plot(ax=ax, label='20 Hz')
ax.legend(loc='upper left', bbox_to_anchor=(0, 1.1), ncol=3)

<matplotlib.legend.Legend at 0x715038f11970>

Above is using the default Butterworth filtering. To use FIR filtering with a window like Hanning is also simple:

fig, ax = ts.plot(label='mix')
ts.filter(cutoff_freq=[15, 25], method='firwin', window='hann').plot(ax=ax, label='After 15-25 Hz band-pass filter')
ts2.plot(ax=ax, label='20 Hz')
ax.legend(loc='upper left', bbox_to_anchor=(0, 1.1), ncol=3)

filtering with windown hann

<matplotlib.legend.Legend at 0x715038f655b0>

Applying a high-pass filter

fig, ax = ts.plot(label='mix')
ts_low  = ts.filter(cutoff_freq=15)
ts_high = ts.copy()
ts_high.value = ts.value - ts_low.value # subtract low-pass filtered series from original one
ts_high.plot(label='High-pass filter @ 15Hz',ax=ax)
ax.legend(loc='upper left', bbox_to_anchor=(0, 1.1), ncol=3)

<matplotlib.legend.Legend at 0x715038fae180>

flip(axis='value', keep_log=False)[source]

Flips the Series along one or both axes

Parameters:

axis (str, optional) – The axis along which the Series will be flipped. The default is ‘value’. Other acceptable options are ‘time’ or ‘both’.
keep_log (Boolean) – if True, adds this transformation to the series log.

Returns:

new – The flipped series object

Return type:

Series

Examples

ts = pyleo.utils.load_dataset('SOI')
tsf = ts.flip(keep_log=True)

fig, ax = tsf.plot()
tsf.log

({0: 'flip', 'applied': True, 'axis': 'value'},)

classmethod from_csv(path)[source]

Read in Series object from CSV file. Expects a metadata header dealineated by ‘###’ lines, as written by the Series.to_csv() method.

Parameters:

filename (str) – name of the file, e.g. ‘myrecord.csv’
path (str) – directory of the file. Default: current directory, ‘.’

Returns:

pyleoclim Series object containing data and metadata.

Return type:

Series

See also

pyleoclim.Series.to_csv

classmethod from_json(path)[source]

Creates a pyleoclim.Series from a JSON file

The keys in the JSON file must correspond to the parameter associated with a Series object

Parameters:: path (str) – Path to the JSON file
Returns:: ts – A Pyleoclim Series object.
Return type:: pyleoclim.core.series.Series

classmethod from_pandas(ser, metadata, verbose=True)[source]

Class method to create a Series object from a pandas Series.

Parameters:

ser (pandas.Series) – The pandas Series object to convert. The index must be a DatetimeIndex.
metadata (dict) – A dictionary containing metadata for the Series. If ‘time_name’ or ‘value_name’ are not provided, defaults to the name of the index and the name of the pandas Series, respectively.
verbose (bool) – If True (default), will print warning messages if there are any

Returns:

The created Series object.

Return type:

Series

Raises:

ValueError – If the index of the pandas Series is not a DatetimeIndex.

Notes

This method first checks if the index of the pandas Series is a DatetimeIndex. If it is, it converts the index to seconds if it’s not already in that unit. Then it converts the datetime index to a time array using the ‘time_unit’ and ‘time_name’ from the metadata.

The method then returns a new Series object with the converted time and values, and the provided metadata.

gaussianize(keep_log=False)[source]

Gaussianizes the timeseries (i.e. maps its values to a standard normal)

Returns:

new (Series) – The Gaussianized series object
keep_log (Boolean) – if True, adds this transformation to the series log.

References

Emile-Geay, J., and M. Tingley (2016), Inferring climate variability from nonlinear proxies: application to palaeo-enso studies, Climate of the Past, 12 (1), 31–50, doi:10.5194/cp- 12-31-2016.

gkernel(step_style='max', keep_log=False, step_type=None, **kwargs)[source]

Coarse-grain a Series object via a Gaussian kernel.

Like .bin() this technique is conservative and uses the max space between points as the default spacing. Unlike .bin(), gkernel() uses a gaussian kernel to calculate the weighted average of the time series over these intervals.

Note that if the series being examined has very low resolution sections with few points, you may need to tune the parameter for the kernel e-folding scale (h).

Parameters:

step_style (str) – type of timestep: ‘mean’, ‘median’, or ‘max’ of the time increments
keep_log (Boolean) – if True, adds the step type and its keyword arguments to the series log.
kwargs – Arguments for kernel function. See pyleoclim.utils.tsutils.gkernel for details

Returns:

new – The coarse-grained Series object

Return type:

Series

See also

pyleoclim.utils.tsutils.gkernel: application of a Gaussian kernel

global_coherence(target_series=None, coh=None, method='cwt', wavelet_kwargs=None)[source]

Compute the global coherence between two Series objects using wavelet analysis

Built on top of pyleoclim.core.series.wavelet_coherence.

Parameters:

target_series (pyleo.Series) – The target series to compare with
coh (pyleo.core.coherence.Coherence) – A coherence object. If None, will compute the coherence from scratch. If passed, will take precedence over target_series.
method (str; {'ww','cwt'}) – The method to use for the wavelet analysis. Default is ‘cwt’, which only works if the series share an evenly spaced time axis.
wavelet_kwargs (dict) – Keyword arguments to pass to the pyleoclim.core.series.wavelet_coherence

Returns:

coh

Return type:

pyleo.core.globalcoherence.GlobalCoherence

See also

pyleoclim.core.globalcoherence.GlobalCoherence: global Coherence object
pyleoclim.core.series.wavelet_coherence: Wavelet coherence analysis

Examples

soi = pyleo.utils.load_dataset('SOI')
nino3 = pyleo.utils.load_dataset('NINO3')

gcoh = soi.global_coherence(nino3)
gcoh.plot()

(<Figure size 800x800 with 2 Axes>,
 <Axes: xlabel='Period [year]', ylabel='PSD'>)

histplot(figsize=[10, 4], title=None, savefig_settings=None, ax=None, ylabel='KDE', vertical=False, edgecolor='w', **plot_kwargs)[source]

Plot the distribution of the timeseries values

Parameters:

figsize (list) – a list of two integers indicating the figure size
title (str) – the title for the figure
savefig_settings (dict) –
the dictionary of arguments for plt.savefig(); some notes below:
- ”path” must be specified; it can be any existed or non-existed path, with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”}
ax (matplotlib.axis, optional) – A matplotlib axis
ylabel (str) – Label for the count axis
vertical ({True,False}) – Whether to flip the plot vertically
edgecolor (matplotlib.color) – The color of the edges of the bar
plot_kwargs (dict) – Plotting arguments for seaborn histplot: https://seaborn.pydata.org/generated/seaborn.histplot.html

See also

pyleoclim.utils.plotting.savefig: saving figure in Pyleoclim

Examples

Distribution of the SOI record

ts = pyleo.utils.load_dataset('SOI')
fig, ax = ts.histplot()

interp(method='linear', keep_log=False, **kwargs)[source]

Interpolate a Series object onto a new time axis

Parameters:

method ({‘linear’, ‘nearest’, ‘zero’, ‘slinear’, ‘quadratic’, ‘cubic’, ‘previous’, ‘next’}) – where ‘zero’, ‘slinear’, ‘quadratic’ and ‘cubic’ refer to a spline interpolation of zeroth, first, second or third order; ‘previous’ and ‘next’ simply return the previous or next value of the point) or as an integer specifying the order of the spline interpolator to use. Default is ‘linear’.
keep_log (Boolean) – if True, adds the method name and its parameters to the series log.
kwargs – Arguments specific to each interpolation function. See pyleoclim.utils.tsutils.interp for details

Returns:

new – An interpolated Series object

Return type:

Series

See also

pyleoclim.utils.tsutils.interp: interpolation function

is_evenly_spaced(tol=0.001)[source]

Check if the Series time axis is evenly-spaced, within tolerance

Parameters:: tol (float) – tolerance. If time increments are all within tolerance, the series is declared evenly-spaced. default = 1e-3
Returns:: res
Return type:: bool

make_labels()[source]

Initialization of plot labels based on Series metadata

Returns:

time_header (str) – Label for the time axis
value_header (str) – Label for the value axis

outliers(method='kmeans', remove=True, settings=None, fig_outliers=True, figsize_outliers=[10, 4], plotoutliers_kwargs=None, savefigoutliers_settings=None, fig_clusters=True, figsize_clusters=[10, 4], plotclusters_kwargs=None, savefigclusters_settings=None, keep_log=False)[source]

Remove outliers from timeseries data. The method employs clustering to identify clusters in the data, using the k-means and DBSCAN algorithms from scikit-learn. Points falling a certain distance from the cluster (either away from the centroid for k-means or in a area of low density for DBSCAN) are considered outliers. The silhouette score is used to optimize parameter values.

A tutorial explaining how to use this method and set the parameters is available at https://github.com/LinkedEarth/PyleoTutorials/blob/main/notebooks/L2_outliers_detection.ipynb.

Parameters:

method (str, {'kmeans','DBSCAN'}, optional) – The clustering method to use. The default is ‘kmeans’.
remove (bool, optional) – If True, removes the outliers. The default is True.
settings (dict, optional) – Specific arguments for the clustering functions. The default is None.
fig_outliers (bool, optional) – Whether to display the timeseries showing the outliers. The default is True.
figsize_outliers (list, optional) – The dimensions of the outliers figure. The default is [10,4].
plotoutliers_kwargs (dict, optional) – Arguments for the plot displaying the outliers. The default is None.
savefigoutliers_settings (dict, optional) –
Saving options for the outlier plot. The default is None. - “path” must be specified; it can be any existed or non-existed path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”}
fig_clusters (bool, optional) – Whether to display the clusters. The default is True.
figsize_clusters (list, optional) – The dimensions of the cluster figures. The default is [10,4].
plotclusters_kwargs (dict, optional) – Arguments for the cluster plot. The default is None.
savefigclusters_settings (dict, optional) –
Saving options for the cluster plot. The default is None. - “path” must be specified; it can be any existed or non-existed path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”}
keep_log (Boolean) – if True, adds the previous method parameters to the series log.

Returns:

ts – A new Series object without outliers if remove is True. Otherwise, returns the original timeseries

Return type:

Series

See also

pyleoclim.utils.tsutils.detect_outliers_DBSCAN: Outlier detection using the DBSCAN method
pyleoclim.utils.tsutils.detect_outliers_kmeans: Outlier detection using the kmeans method
pyleoclim.utils.tsutils.remove_outliers: Remove outliers from the series

Examples

LR04 = pyleo.utils.load_dataset('LR04')
LR_out = LR04.detrend().standardize().outliers(method='kmeans')

To set the number of clusters:

LR_out = LR04.detrend().standardize().outliers(method='kmeans', settings={'nbr_clusters':2})

The log contains diagnostic information, to access it, set the keep_log parameter to True:

LR_out = LR04.detrend().standardize().outliers(method='kmeans', settings={'nbr_clusters':2}, keep_log=True)

overlap(ts)[source]

Convenience method to check the degree of overlap between two Series objects

Parameters:: ts (pyleo.Series) – The target series
Returns:: ovrlp – length of overlap (negative means overlap deficit) in units of the original object
Return type:: float

See also

pyleoclim.utils.tsbase.overlap: compute length of overlap

pandas_method(method)[source]

Apply a pandas method to the Series object

Parameters:: method (str) – The name of the pandas method to apply to the Series object
Returns:: A new Series object with the result of the method applied to the original Series object
Return type:: Series

plot(figsize=[10, 4], marker=None, markersize=None, color=None, linestyle=None, linewidth=None, xlim=None, ylim=None, label=None, xlabel=None, ylabel=None, title=None, zorder=None, legend=True, plot_kwargs=None, lgd_kwargs=None, alpha=None, savefig_settings=None, ax=None, invert_xaxis=False, invert_yaxis=False)[source]

Plot the timeseries

Parameters:

figsize (list) – a list of two integers indicating the figure size
marker (str) – e.g., ‘o’ for dots See [matplotlib.markers](https://matplotlib.org/stable/api/markers_api.html) for details
markersize (float) – the size of the marker
color (str, list) – the color for the line plot e.g., ‘r’ for red See [matplotlib colors](https://matplotlib.org/stable/gallery/color/color_demo.html) for details
linestyle (str) – e.g., ‘–’ for dashed line See [matplotlib.linestyles](https://matplotlib.org/stable/gallery/lines_bars_and_markers/linestyles.html) for details
linewidth (float) – the width of the line
label (str) – the label for the line
xlabel (str) – the label for the x-axis
ylabel (str) – the label for the y-axis
title (str) – the title for the figure
zorder (int) – The default drawing order for all lines on the plot
legend ({True, False}) – plot legend or not
invert_xaxis (bool, optional) – if True, the x-axis of the plot will be inverted
invert_yaxis (bool, optional) – same for the y-axis
plot_kwargs (dict) – the dictionary of keyword arguments for ax.plot() See [matplotlib.pyplot.plot](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.plot.html) for details
lgd_kwargs (dict) – the dictionary of keyword arguments for ax.legend() See [matplotlib.pyplot.legend](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.legend.html) for details
alpha (float) – Transparency setting
savefig_settings (dict) –
the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existed or non-existed path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”}
ax (matplotlib.axis, optional) – the axis object from matplotlib See [matplotlib.axes](https://matplotlib.org/api/axes_api.html) for details.

Returns:

fig (matplotlib.figure) – the figure object from matplotlib See [matplotlib.pyplot.figure](https://matplotlib.org/stable/api/figure_api.html) for details.
ax (matplotlib.axis) – the axis object from matplotlib See [matplotlib.axes](https://matplotlib.org/stable/api/axes_api.html) for details.

Notes

When ax is passed, the return will be ax only; otherwise, both fig and ax will be returned.

See also

pyleoclim.utils.plotting.savefig: saving a figure in Pyleoclim

Examples

Plot the SOI record

ts = pyleo.utils.load_dataset('SOI')
fig, ax = ts.plot()

Change the line color

fig, ax = ts.plot(color='r')

Save the figure. Two options available, only one is needed:

Within the plotting command
After the figure has been generated

fig, ax = ts.plot(color='k', savefig_settings={'path': 'ts_plot3.png'}); pyleo.closefig(fig)
pyleo.savefig(fig,path='ts_plot3.png')

Figure saved at: "ts_plot3.png"
Figure saved at: "ts_plot3.png"

Plot the CENOGRID d18O timeseries with added Geologic Time Scale annotation

import pyleoclim as pyleo

ts_18 = pyleo.utils.load_dataset('cenogrid_d18O')
fig, ax = ts_18.plot(figsize=(10, 4),linewidth=0.5)
ax.invert_yaxis() # d18O is traditionally inverted

fig, ax = pyleo.add_GTS(fig, ax, ranks=['Period', 'Epoch'], location='above')

Using rdflib to load GTS information from ICS

resample(rule, keep_log=False, **kwargs)[source]

Analogue to pandas.Series.resample.

This is a convenience method: doing

ser.resample(‘YS’).mean()

will do the same thing as

ser.pandas_method(lambda x: x.resample(‘YS’).mean())

but will also accept some extra resampling rules, such as ‘Ga’ (see below).

NOTE: this feature may break until [this pandas bug](https://github.com/pandas-dev/pandas/issues/57427) is fixed.

Parameters:

rule (str) –
The offset string or object representing target conversion. Can also accept pyleoclim units, such as ‘ka’ (1000 years), ‘Ma’ (1 million years), and ‘Ga’ (1 billion years).

Check the [pandas resample docs](https://pandas.pydata.org/docs/dev/reference/api/pandas.DataFrame.resample.html) for more details.
kwargs (dict) – Any other arguments which will be passed to pandas.Series.resample.

Returns:

Resampler object, not meant to be used to directly. Instead, an aggregation should be called on it, see examples below.

Return type:

SeriesResampler

Examples

ts = pyleo.utils.load_dataset('LR04')
ts5k = ts.resample('5ka').mean()
fig, ax = ts.plot(invert_yaxis='True',xlim=[0, 1000])
ts5k.plot(ax=ax,color='C1')

<Axes: xlabel='Age [ky BP]', ylabel='$\\delta^{18} \\mathrm{O}$ [‰]'>

resolution()[source]

Generate a resolution object

Increments are assigned to the preceding time value. E.g. for time_axis = [0,1,3], resolution.resolution = [1,2] resolution.time = [0,1]

Returns:: resolution – Resolution object
Return type:: Resolution

See also

pyleoclim.core.resolutions.Resolution

Examples

To create a resolution object, apply the .resolution() method to a Series object

ts = pyleo.utils.load_dataset('EDC-dD')
resolution = ts.resolution()

Several methods are then available:

Summary statistics can be obtained via .describe()

resolution.describe()

{'nobs': np.int64(5784),
 'minmax': (np.float64(8.244210000000002), np.float64(1364.0)),
 'mean': np.float64(138.5932963710235),
 'variance': np.float64(29806.73648249974),
 'skewness': np.float64(2.661861461835658),
 'kurtosis': np.float64(8.705801510819656),
 'median': np.float64(58.132250000006024)}

A simple plot can be created using .plot()

resolution.plot()

(<Figure size 1000x400 with 1 Axes>,
 <Axes: xlabel='Age [y BP]', ylabel='resolution [y BP]'>)

The distribution of resolution

resolution.histplot()

(<Figure size 1000x400 with 1 Axes>,
 <Axes: xlabel='resolution [y BP]', ylabel='KDE'>)

Or a dashboard combining plot() and histplot() side by side:

resolution.dashboard()

(<Figure size 1100x800 with 2 Axes>,
 {'res': <Axes: xlabel='Age [y BP]', ylabel='resolution [y BP]'>,
  'res_hist': <Axes: xlabel='Counts'>})

segment(factor=10, verbose=False)[source]

Gap detection

This function segments a timeseries into n number of parts following a gap: detection algorithm. The rule of gap detection is very simple: we define the intervals between time points as dts, then if dts[i] is larger than factor * dts[i-1], we think that the change of dts (or the gradient) is too large, and we regard it as a breaking point and divide the time series into two segments here

Parameters:

factor (float) – The factor that adjusts the threshold for gap detection
verbose (bool) – If True, will print warning messages if there is any

Returns:

res – If gaps were detected, returns the segments in a MultipleSeries object, else, returns the original timeseries.

Return type:

MultipleSeries or Series

sel(value=None, time=None, tolerance=0)[source]

Slice Series based on ‘value’ or ‘time’.

Parameters:

value (int, float, slice) – If int/float, then the Series will be sliced so that self.value is equal to value (+/- tolerance). If slice, then the Series will be sliced so self.value is between slice.start and slice.stop (+/- tolerance).
time (int, float, slice) – If int/float, then the Series will be sliced so that self.time is equal to time. (+/- tolerance) If slice of int/float, then the Series will be sliced so that self.time is between slice.start and slice.stop. If slice of datetime (or str containing datetime, such as ‘2020-01-01’), then the Series will be sliced so that self.datetime_index is between time.start and time.stop (+/- tolerance, which needs to be a timedelta).
tolerance (int, float, default 0.) – Used by value and time, see above.

Return type:

Copy of self, sliced according to value and time.

Examples

>>> ts = pyleo.Series(
...     time=np.array([1, 1.1, 2, 3]), value=np.array([4, .9, 6, 1]), time_unit='years BP'
... )
>>> ts.sel(value=1)
{'log': ({0: 'clean_ts', 'applied': True, 'verbose': False},
        {2: 'clean_ts', 'applied': True, 'verbose': False})}

None time [years BP] 3.0 1.0 Name: value, dtype: float64

If you also want to include the value 3.9, you could set tolerance to .1:

>>> ts.sel(value=1, tolerance=.1)
{'log': ({0: 'clean_ts', 'applied': True, 'verbose': False},
        {2: 'clean_ts', 'applied': True, 'verbose': False})}

None time [years BP] 1.1 0.9 3.0 1.0 Name: value, dtype: float64

You can also pass a slice to select a range of values:

>>> ts.sel(value=slice(4, 6))
{'log': ({0: 'clean_ts', 'applied': True, 'verbose': False},
        {2: 'clean_ts', 'applied': True, 'verbose': False})}

None time [years BP] 1.0 4.0 2.0 6.0 Name: value, dtype: float64

>>> ts.sel(value=slice(4, None))
{'log': ({0: 'clean_ts', 'applied': True, 'verbose': False},
        {2: 'clean_ts', 'applied': True, 'verbose': False})}

None time [years BP] 1.0 4.0 2.0 6.0 Name: value, dtype: float64

>>> ts.sel(value=slice(None, 4))
{'log': ({0: 'clean_ts', 'applied': True, 'verbose': False},
        {2: 'clean_ts', 'applied': True, 'verbose': False})}

None time [years BP] 1.0 4.0 1.1 0.9 3.0 1.0 Name: value, dtype: float64

Similarly, you filter using time instead of value.

slice(timespan)[source]

Slicing the timeseries with a timespan (tuple or list)

Parameters:: timespan (tuple or list) – The list of time points for slicing, whose length must be even. When there are n time points, the output Series includes n/2 segments. For example, if timespan = [a, b], then the sliced output includes one segment [a, b]; if timespan = [a, b, c, d], then the sliced output includes segment [a, b] and segment [c, d].
Returns:: new – The sliced Series object.
Return type:: Series

Examples

slice the SOI from 1972 to 1998

ts = pyleo.utils.load_dataset('SOI')
ts_slice = ts.slice([1972, 1998])
print("New time bounds:",ts_slice.time.min(),ts_slice.time.max())

New time bounds: 1972.0010513731472 1997.9181004754107

sort(verbose=False, ascending=True, keep_log=False)[source]

Ensure timeseries is set to a monotonically increasing axis.: If the time axis is prograde to begin with, no transformation is applied.

Parameters:

verbose (bool) – If True, will print warning messages if there is any
keep_log (Boolean) – if True, adds this step and its parameter to the series log.

Returns:

new – Series object with removed NaNs and sorting

Return type:

Series

spectral(method='lomb_scargle', freq=None, freq_kwargs=None, settings=None, label=None, scalogram=None, verbose=False)[source]

Perform spectral analysis on the timeseries

Parameters:

method (str;) – {‘wwz’, ‘mtm’, ‘lomb_scargle’, ‘welch’, ‘periodogram’, ‘cwt’} Default is Lomb-Scargle, because it can handle unevenly spaced series, and is fast. However, for evenly spaced data mtm would almost surely be a better choice.
freq (str or array, optional) – Information to produce the frequency vector. This can be ‘log’,’scale’, ‘nfft’, ‘lomb_scargle’ or ‘welch’ or a NumPy array. If a string, will use make_freq_vector with the specified frequency-generating method. If an array, this will be passed directly to the spectral method. If None (default), will use ‘log’ for WWZ and ‘lomb_scargle’ for Lomb-Scargle. This parameter is highly consequential for the WWZ and Lomb-Scargle methods, but is otherwise ignored, as other spectral methods generate their frequency vector internally.
freq_kwargs (dict) – Arguments for frequency vector
settings (dict) – Arguments for the specific spectral method
label (str) – Label for the PSD object
scalogram (pyleoclim.core.series.Series.Scalogram) – The return of the wavelet analysis; effective only when the method is ‘wwz’ or ‘cwt’
verbose (bool) – If True, will print warning messages if there is any

Returns:

psd – A PSD object

Return type:

PSD

See also

pyleoclim.utils.spectral.mtm: Spectral analysis using the Multitaper approach
pyleoclim.utils.spectral.lomb_scargle: Spectral analysis using the Lomb-Scargle method
pyleoclim.utils.spectral.welch: Spectral analysis using the Welch segement approach
pyleoclim.utils.spectral.periodogram: Spectral anaysis using the basic Fourier transform
pyleoclim.utils.spectral.wwz_psd: Spectral analysis using the Wavelet Weighted Z transform
pyleoclim.utils.spectral.cwt_psd: Spectral analysis using the continuous Wavelet Transform as implemented by Torrence and Compo
pyleoclim.utils.spectral.make_freq_vector: Functions to create the frequency vector
pyleoclim.utils.tsutils.detrend: Detrending function
pyleoclim.core.psds.PSD: PSD object
pyleoclim.core.psds.MultiplePSD: Multiple PSD object

Examples

Calculate the spectrum of SOI using the various methods and compute significance

ts = pyleo.utils.load_dataset('SOI')
ts_std = ts.standardize()

Lomb-Scargle

psd_ls = ts_std.spectral(method='lomb_scargle')
psd_ls_signif = psd_ls.signif_test(number=20) #in practice, need more AR(1) simulations
fig, ax = psd_ls_signif.plot(title='PSD using Lomb-Scargle method')

We may pass in method-specific arguments via “settings”, which is a dictionary. For instance, to adjust the number of overlapping segment for Lomb-Scargle, we may specify the method-specific argument “n50”; to adjust the frequency vector, we may modify the “freq” or modify the method-specific argument “freq”.

psd_LS_n50 = ts_std.spectral(method='lomb_scargle', settings={'n50': 4})  # c=1e-2 yields lower frequency resolution
psd_LS_freq = ts_std.spectral(method='lomb_scargle',freq=np.linspace(1/20, 1/0.2, 51))
psd_LS_log = ts_std.spectral(method='lomb_scargle', freq='log')  # with frequency vector generated using REDFIT method
fig, ax = psd_LS_n50.plot(
    title='Lomb-Scargle PSD with 4 overlapping segments',
    label='settings={"n50": 4}')
psd_ls.plot(ax=ax, label='settings={"n50": 3}', marker='o', alpha=.2)

fig, ax = psd_LS_freq.plot(
    title='Lomb-Scargle PSD with specified frequency vector',
    label='freq=np.linspace(1/20, 1/0.2, 51)', marker='o')
psd_ls.plot(ax=ax, label='Lomb-Scargle frequency vector', marker='o', alpha=.2)

fig, ax = psd_LS_log.plot(
    title='Lomb-Scargle PSD with different frequency vectors',
    label='log frequency method', marker='o')
psd_ls.plot(ax=ax, label='Lomb-Scargle frequency vector', marker='o', alpha=.2)

<Axes: title={'center': 'Lomb-Scargle PSD with different frequency vectors'}, xlabel='Period [year]', ylabel='PSD'>

You may notice the differences in the PSD curves regarding smoothness and the locations of the analyzed period points.

For other method-specific arguments, please look up the specific methods in the “See also” section.

WWZ

psd_wwz = ts_std.spectral(method='wwz')  # wwz is the default method
psd_wwz_signif = psd_wwz.signif_test(number=1)  # significance test; for real work, should use number=200 or even larger
fig, ax = psd_wwz_signif.plot(title='PSD using WWZ method')

We may take advantage of a pre-calculated scalogram using WWZ to accelerate the spectral analysis (although note that the default parameters for spectral and wavelet analysis using WWZ are different):

scal_wwz = ts_std.wavelet(method='wwz')  # wwz is the default method
psd_wwz_fast = ts_std.spectral(method='wwz', scalogram=scal_wwz)
fig, ax = psd_wwz_fast.plot(title='PSD using WWZ method w/ pre-calculated scalogram')

Periodogram

ts_interp = ts_std.interp()
psd_perio = ts_interp.spectral(method='periodogram')
psd_perio_signif = psd_perio.signif_test(number=20, method='ar1sim') #in practice, need more AR(1) simulations
fig, ax = psd_perio_signif.plot(title='PSD using Periodogram method')

Welch

psd_welch = ts_interp.spectral(method='welch')
psd_welch_signif = psd_welch.signif_test(number=20, method='ar1sim') #in practice, need more AR(1) simulations
fig, ax = psd_welch_signif.plot(title='PSD using Welch method')

MTM

psd_mtm = ts_interp.spectral(method='mtm', label='MTM, NW=4')
psd_mtm_signif = psd_mtm.signif_test(number=20, method='ar1sim') #in practice, need more AR(1) simulations
fig, ax = psd_mtm_signif.plot(title='PSD using the multitaper method')

By default, MTM uses a half-bandwidth of 4 times the fundamental (Rayleigh) frequency, i.e. NW = 4, which is the most conservative choice. NW runs from 2 to 4 in multiples of 1/2, and can be adjusted like so (note the sharper peaks and higher overall variance, which may not be desirable):

psd_mtm2 = ts_interp.spectral(method='mtm', settings={'NW':2}, label='MTM, NW=2')
fig, ax = psd_mtm2.plot(title='MTM with NW=2')

Continuous Wavelet Transform

ts_interp = ts_std.interp()
psd_cwt = ts_interp.spectral(method='cwt')
psd_cwt_signif = psd_cwt.signif_test(number=20)
fig, ax = psd_cwt_signif.plot(title='PSD using the CWT method')

ssa(M=None, nMC=0, f=0.3, trunc=None, var_thresh=80, online=True)[source]

Singular Spectrum Analysis

Nonparametric, orthogonal decomposition of timeseries into constituent oscillations. This implementation uses the method of [1], with applications presented in [2]. Optionally (MC>0), the significance of eigenvalues is assessed by Monte-Carlo simulations of an AR(1) model fit to X, using [3]. The method expects regular spacing, but is tolerant to missing values, up to a fraction 0<f<1 (see [4]).

Parameters:

M (int, optional) – window size. The default is None , which will use 10% of the length of the series.
MC (int, optional) – Number of iteration in the Monte-Carlo process. The default is 0.
f (float, optional) – maximum allowable fraction of missing values. The default is 0.3.
trunc (str) –
if present, truncates the expansion to a level K < M owing to one of 4 criteria:
1. ’kaiser’: variant of the Kaiser-Guttman rule, retaining eigenvalues larger than the median
2. ’mcssa’: Monte-Carlo SSA (use modes above the 95% quantile from an AR(1) process)
3. ’var’: first K modes that explain at least var_thresh % of the variance. Default is None, which bypasses truncation (K = M)
4. ’knee’: Wherever the “knee” of the screeplot occurs. (See kneed’s documentation)
The knee method is recommended as a first pass at identifying significant modes as it tends to be more robust than ‘kaiser’ or ‘var’, and faster than ‘mcssa’. While no truncation method is imposed by default, if the goal is to enhance the S/N ratio and reconstruct a smooth version of the attractor’s skeleton, then the knee-finding method is a good compromise between objectivity and efficiency.
var_thresh (float) – variance threshold for reconstruction (only impactful if trunc is set to ‘var’)
online (bool; {True,False}) –
Whether or not to conduct knee finding analysis online or offline. Only called when trunc = ‘knee’. Default is True See kneed’s documentation for details.

Returns:

res (object of the SsaRes class containing:)
eigvals ((M, ) array of eigenvalues)
eigvecs ((M, M) Matrix of temporal eigenvectors (T-EOFs))
PC ((N - M + 1, M) array of principal components (T-PCs))
RCmat ((N, M) array of reconstructed components)
RCseries ((N,) reconstructed series, with mean and variance restored (same type as original))
pctvar ((M, ) array of the fraction of variance (%) associated with each mode)
eigvals_q ((M, 2) array contaitning the 5% and 95% quantiles of the Monte-Carlo eigenvalue spectrum [ if nMC >0 ])

References

[1]

Vautard, R., and M. Ghil (1989), Singular spectrum analysis in nonlinear dynamics, with applications to paleoclimatic time series, Physica D, 35, 395–424.

[2]

Ghil, M., R. M. Allen, M. D. Dettinger, K. Ide, D. Kondrashov, M. E. Mann, A. Robertson, A. Saunders, Y. Tian, F. Varadi, and P. Yiou (2002), Advanced spectral methods for climatic time series, Rev. Geophys., 40(1), 1003–1052, doi:10.1029/2000RG000092.

[3] (1,2)

Allen, M. R., and L. A. Smith (1996), Monte Carlo SSA: Detecting irregular oscillations in the presence of coloured noise, J. Clim., 9, 3373–3404.

[4]

Schoellhamer, D. H. (2001), Singular spectrum analysis for time series with missing data, Geophysical Research Letters, 28(16), 3187–3190, doi:10.1029/2000GL012698.

See also

pyleoclim.core.utils.decomposition.ssa: Singular Spectrum Analysis utility
pyleoclim.core.ssares.SsaRes.modeplot: plot SSA modes
pyleoclim.core.ssares.SsaRes.screeplot: plot SSA eigenvalue spectrum

Examples

SSA with SOI

ts = pyleo.utils.load_dataset('NINO3')
fig, ax = ts.plot()
nino_ssa = ts.ssa(M=60)

The result of this operation yields an SsaRes object, which comes with dedicated methods. To inspect the eigenvalue spectrum, we use scree_plot().

fig, ax = nino_ssa.screeplot()

This highlights a few common properties of SSA:

the eigenvalues are in descending order
their uncertainties are proportional to the eigenvalues themselves
the eigenvalues tend to come in pairs : (1,2) (3,4), are all clustered
after a point (here i=15), the eigenvalues appear to reach a floor, explaining a very small amount of variance.

Summing the variance of the first 14 modes, we get:

print(nino_ssa.pctvar[:14].sum())

95.57806552633473

That is a typical result for a (paleo)climate timeseries; a few modes do the vast majority of the work. This means that we can focus our attention on these modes and capture most of the interesting behavior. To see this, let’s use the reconstructed components (RCs), and sum the RC matrix over the first 14 columns:

RCmat = nino_ssa.RCmat[:,:14]
RCk = (RCmat-RCmat.mean()).sum(axis=1) + ts.value.mean()
fig, ax = ts.plot(title='NINO3')
ax.plot(nino_ssa.orig.time,RCk,label='SSA reconstruction, 14 modes',color='orange')
ax.legend()

<matplotlib.legend.Legend at 0x715021847dd0>

Indeed, these first few modes capture the vast majority of the low-frequency behavior, including all the El Niño/La Niña events. What is left (the blue wiggles not captured in the orange curve) are high-frequency oscillations that might be considered “noise” from the standpoint of ENSO dynamics. This illustrates how SSA might be used for filtering a timeseries. One must be careful however:

there was not much rhyme or reason for picking 14 modes. Why not 5, or 39? All we have seen so far is that they gather >95% of the variance, which is by no means a magic number.

there is no guarantee that the first few modes will filter out high-frequency behavior, or at what frequency cutoff they will do so. If you need to cut out specific frequencies, you are better off doing it with a classical filter, like the butterworth filter implemented in Pyleoclim. However, in many instances the choice of a cutoff frequency is itself rather arbitrary. In such cases, SSA provides a principled alternative for generating a version of a timeseries that preserves dominant features and excludes the rest.

as with all orthgonal decompositions, summing over all RCs will recover the original signal within numerical precision.

Selecting meaningful modes in eigenproblems (e.g. EOF analysis) is more art than science. However, one technique stands out: Monte Carlo SSA, introduced by [3] to identify SSA modes that rise above what one would expect from “red noise”, specifically an AR(1) process). To run it, simply provide the parameter nMC, ideally large enough to get decent statistics. Here let’s use nMC = 1000. The result will be stored in the eigval_q array, which has the same length as eigval, and its two columns contain the 5% and 95% quantiles of the ensemble of MC-SSA eigenvalues.

nino_mcssa = ts.ssa(M = 60, nMC=1000)

Now let’s look at the result:

fig, ax = nino_mcssa.screeplot()
print('Indices of modes retained: '+ str(nino_mcssa.mode_idx))

Indices of modes retained: [0 1 5 6 7 8 9]

This suggests that modes 1-2 fall above the red noise benchmark, but so 5-9. To inspect mode 1 (index 0), just type:

fig, ax = nino_mcssa.modeplot(index=0)

To inspect the reconstructed series, simply do:

fig, ax = ts.plot()
nino_mcssa.RCseries.plot(ax=ax)

<Axes: xlabel='Time [year C.E.]', ylabel='NINO3 [$^{\\circ}$C]'>

For more details, see the PyleoTutorial <http://linked.earth/PyleoTutorials/notebooks/L2_singular_spectrum_analysis.html>_

standardize(keep_log=False, scale=1)[source]

Standardizes the series ((i.e. remove its estimated mean and divides by its estimated standard deviation)

Returns:

new (Series) – The standardized series object
keep_log (Boolean) – if True, adds the previous mean, standard deviation and method parameters to the series log.

stats()[source]

Compute basic statistics from a Series

Computes the mean, median, min, max, standard deviation, and interquartile range of a numpy array y, ignoring NaNs.

Returns:: res – Contains the mean, median, minimum value, maximum value, standard deviation, and interquartile range for the Series.
Return type:: dictionary

Examples

Compute basic statistics for the SOI series

ts = pyleo.utils.load_dataset('SOI')
ts.stats()

{'mean': np.float64(-0.045235602094240844),
 'median': np.float64(-0.02),
 'min': np.float64(-4.34),
 'max': np.float64(4.07),
 'std': np.float64(1.1298792605604073),
 'IQR': np.float64(1.52)}

stripes(figsize=[8, 1], cmap='RdBu_r', ref_period=None, sat=1.0, top_label=None, bottom_label=None, label_color='gray', label_size=None, xlim=None, xlabel=None, savefig_settings=None, ax=None, invert_xaxis=False, show_xaxis=False, x_offset=0.03)[source]

Represents the Series as an Ed Hawkins “stripes” pattern

Credit: https://matplotlib.org/matplotblog/posts/warming-stripes/

Parameters:

ref_period (array-like (2-elements)) – dates of the reference period, in the form “(first, last)”
figsize (list) – a list of two integers indicating the figure size (in inches)
cmap (str) – colormap name (https://matplotlib.org/stable/tutorials/colors/colormaps.html)
sat (float > 0) – Controls the saturation of the colormap normalization by scaling the vmin, vmax in https://matplotlib.org/stable/tutorials/colors/colormapnorms.html default = 0.9
xlim (list) – time axis limits
top_label (str) – the “title” label for the stripe
bottom_label (str) – the “ylabel” explaining which variable is being plotted
invert_xaxis (bool, optional) – if True, the x-axis of the plot will be inverted
x_offset (float) – value controlling the horizontal offset between stripes and labels (default = 0.05)
show_xaxis (bool) – flag indicating whether or not the x-axis should be shown (default = False)
savefig_settings (dict) –
the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existed or non-existed path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”}
ax (matplotlib.axis, optional) – the axis object from matplotlib See [matplotlib.axes](https://matplotlib.org/api/axes_api.html) for details.

Returns:

fig (matplotlib.figure) – the figure object from matplotlib See [matplotlib.pyplot.figure](https://matplotlib.org/stable/api/figure_api.html) for details.
ax (matplotlib.axis) – the axis object from matplotlib See [matplotlib.axes](https://matplotlib.org/stable/api/axes_api.html) for details.

Notes

When ax is passed, the return will be ax only; otherwise, both fig and ax will be returned.

See also

pyleoclim.utils.plotting.stripes: stripes representation of a timeseries
pyleoclim.utils.plotting.savefig: saving a figure in Pyleoclim

Examples

Plot the HadCRUT5 Global Mean Surface Temperature

gmst = pyleo.utils.load_dataset('HadCRUT5')
fig, ax = gmst.stripes(ref_period=(1971,2000))

For a more pastel tone, dial down saturation:

fig, ax = gmst.stripes(ref_period=(1971,2000), sat = 0.8)

To change the colormap:

fig, ax = gmst.stripes(ref_period=(1971,2000), cmap='Spectral_r')
fig, ax = gmst.stripes(ref_period=(1971,2000), cmap='magma_r')

To show the time axis:

fig, ax = gmst.stripes(ref_period=(1971,2000), show_xaxis=True)

summary_plot(psd, scalogram, figsize=[8, 10], title=None, time_lim=None, value_lim=None, period_lim=None, psd_lim=None, time_label=None, value_label=None, period_label=None, psd_label=None, ts_plot_kwargs=None, wavelet_plot_kwargs=None, psd_plot_kwargs=None, gridspec_kwargs=None, y_label_loc=None, legend=None, savefig_settings=None)[source]

Produce summary plot of timeseries.

Generate cohesive plot of timeseries alongside results of wavelet analysis and spectral analysis on said timeseries. Requires wavelet and spectral analysis to be conducted outside of plotting function, psd and scalogram must be passed as arguments.

Parameters:

psd (PSD) – the PSD object of a Series.
scalogram (Scalogram) – the Scalogram object of a Series. If the passed scalogram object contains stored signif_scals these will be plotted.
figsize (list) – a list of two integers indicating the figure size
title (str) – the title for the figure
time_lim (list or tuple) – the limitation of the time axis. This is for display purposes only, the scalogram and psd will still be calculated using the full time series.
value_lim (list or tuple) – the limitation of the value axis of the timeseries. This is for display purposes only, the scalogram and psd will still be calculated using the full time series.
period_lim (list or tuple) – the limitation of the period axis
psd_lim (list or tuple) – the limitation of the psd axis
time_label (str) – the label for the time axis
value_label (str) – the label for the value axis of the timeseries
period_label (str) – the label for the period axis
psd_label (str) – the label for the amplitude axis of PDS
legend (bool) – if set to True, a legend will be added to the open space above the psd plot
ts_plot_kwargs (dict) – arguments to be passed to the timeseries subplot, see Series.plot for details
wavelet_plot_kwargs (dict) – arguments to be passed to the scalogram plot, see pyleoclim.Scalogram.plot for details
psd_plot_kwargs (dict) –
arguments to be passed to the psd plot, see PSD.plot for details Certain psd plot settings are required by summary plot formatting. These include:
- ylabel
- legend
- tick parameters
These will be overriden by summary plot to prevent formatting errors
gridspec_kwargs (dict) –
arguments used to build the specifications for gridspec configuration The plot is constructed with six slots:
- slot [0] contains a subgridspec containing the timeseries and scalogram (shared x axis)
- slot [1] contains a subgridspec containing an empty slot and the PSD plot (shared y axis with scalogram)
- slot [2] and slot [3] are empty to allow ample room for xlabels for the scalogram and PSD plots
- slot [4] contains the scalogram color bar
- slot [5] is empty
It is possible to tune the size and spacing of the various slots:
- ’width_ratios’: list of two values describing the relative widths of the column containig the timeseries/scalogram/colorbar and the column containig the PSD plot (default: [6, 1])
- ’height_ratios’: list of three values describing the relative heights of the three timeseries, scalogram and colorbar (default: [2, 7, .35])
- ’hspace’: vertical space between timeseries and scalogram (default: 0, however if either the scalogram xlabel or the PSD xlabel contain ‘n’, .05)
- ’wspace’: lateral space between scalogram and psd plot (default: 0)
- ’cbspace’: vertical space between the scalogram and colorbar
y_label_loc (float) – Plot parameter to adjust horizontal location of y labels to avoid conflict with axis labels, default value is -0.15
savefig_settings (dict) –
the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existed or non-existed path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”}

See also

pyleoclim.core.series.Series.spectral: Spectral analysis for a timeseries
pyleoclim.core.series.Series.wavelet: Wavelet analysis for a timeseries
pyleoclim.utils.plotting.savefig: saving figure in Pyleoclim
pyleoclim.core.psds.PSD: PSD object
pyleoclim.core.psds.MultiplePSD: Multiple PSD object

Examples

Summary_plot with pre-generated psd and scalogram objects. Note that if the scalogram contains saved noise realizations these will be flexibly reused. See pyleo.Scalogram.signif_test() for details

series = pyleo.utils.load_dataset('NINO3')
psd = series.spectral(freq = 'welch')
scalogram = series.wavelet(freq = 'welch')

fig, ax = series.summary_plot(psd = psd,scalogram = scalogram)
pyleo.closefig(fig)

Summary_plot with pre-generated psd and scalogram objects from before and some plot modification arguments passed. Note that if the scalogram contains saved noise realizations these will be flexibly reused. See pyleo.Scalogram.signif_test() for details

psd = series.spectral(freq = 'welch').signif_test(number=20)
fig, ax = series.summary_plot(psd = psd,scalogram = scalogram, period_lim = [0,5], ts_plot_kwargs = {'color':'Purple','linewidth':.5}, psd_plot_kwargs = {'color':'Purple','linewidth':1.5})
pyleo.closefig(fig)

to_csv(metadata_header=True, path=None)[source]

Export Series to csv

Parameters:

metadata_header (boolean, optional) – DESCRIPTION. The default is True.
path (str, optional) – system path to save the file. Default is ‘.’

Return type:

None

See also

pyleoclim.Series.from_csv

Examples

LR04 = pyleo.utils.load_dataset('LR04')
LR04.to_csv()
lr04 = pyleo.Series.from_csv('LR04_benthic_stack.csv')
LR04.equals(lr04)

Series exported to LR04_benthic_stack.csv
Time axis values sorted in ascending order

(True, True)

to_json(path=None)[source]

Export the pyleoclim.Series object to a json file

Parameters:: path (string, optional) – The path to the file. The default is None, resulting in a file saved in the current working directory using the label for the dataset as filename if available or ‘series.json’ if label is not provided.
Return type:: None.

Examples

ts = pyleo.utils.load_dataset('SOI')
ts.to_json('soi.json')

to_pandas(paleo_style=False)[source]

Export to pandas Series

Parameters:: paleo_style (boolean, optional) – If True, will replace datetime with time and label columns with units . The default is False.
Returns:: ser
Return type:: pd.Series representation of the pyleo.Series object

view()[source]

Generates a DataFrame version of the Series object, suitable for viewing in a Jupyter Notebook

Return type:: pd.DataFrame

Examples

Plot the HadCRUT5 Global Mean Surface Temperature

ts = pyleo.utils.load_dataset('HadCRUT5')
ts.view()

	GMST [$^{\circ}$C]
Time [year C.E.]
1850.0	-0.417659
1851.0	-0.233350
1852.0	-0.229399
1853.0	-0.270354
1854.0	-0.291630
...	...
2018.0	0.762654
2019.0	0.891073
2020.0	0.922794
2021.0	0.761856
2022.0	0.801242

173 rows × 1 columns

wavelet(method='cwt', settings=None, freq=None, freq_kwargs=None, verbose=False)[source]

Perform wavelet analysis on a timeseries

Parameters:

method (str {wwz, cwt}) –

cwt - the continuous wavelet transform [1]
is appropriate for evenly-spaced series.

wwz - the weighted wavelet Z-transform [2]
is appropriate for unevenly-spaced series.

Default is cwt, returning an error if the Series is unevenly-spaced.
freq (str or array, optional) – Information to produce the frequency vector (highly consequential for the WWZ method) This can be ‘log’,’scale’, ‘nfft’, ‘lomb_scargle’, ‘welch’ or a NumPy array. If a string, will use make_freq_vector() with the specified frequency-generating method. If an array, this will be passed directly to the spectral method. If None (default), will use the ‘log’ method
freq_kwargs (dict) – Arguments for the frequency vector
settings (dict) – Arguments for the specific wavelet method
verbose (bool) – If True, will print warning messages if there are any

Returns:

scal

Return type:

Scalogram object

See also

pyleoclim.utils.wavelet.wwz: wwz function
pyleoclim.utils.wavelet.cwt: cwt function
pyleoclim.utils.spectral.make_freq_vector: Functions to create the frequency vector
pyleoclim.utils.tsutils.detrend: Detrending function
pyleoclim.core.series.Series.spectral: spectral analysis tools
pyleoclim.core.scalograms.Scalogram: Scalogram object
pyleoclim.core.scalograms.MultipleScalogram: Multiple Scalogram object

References

[1] Torrence, C. and G. P. Compo, 1998: A Practical Guide to Wavelet Analysis. Bull. Amer. Meteor. Soc., 79, 61-78. Python routines available at http://paos.colorado.edu/research/wavelets/

[2] Foster, G., 1996: Wavelets for period analysis of unevenly sampled time series. The Astronomical Journal, 112, 1709.

Examples

Wavelet analysis on the evenly-spaced NINO3 record. The CWT method will be applied by default.

ts = pyleo.utils.load_dataset('NINO3')
scal1 = ts.wavelet()
scal_signif = scal1.signif_test(number=20)  # for research-grade work, use number=200 or larger
fig, ax = scal_signif.plot()

If you wanted to invoke the WWZ method instead (here with no significance testing, to lower computational cost):

scal2 = ts.wavelet(method='wwz')
fig, ax = scal2.plot()

Notice that the two scalograms have different amplitudes, which are relative. Method-specific arguments may be passed via settings. For instance, if you wanted to change the default mother wavelet (‘MORLET’) to a derivative of a Gaussian (DOG), with degree 2 by default (“Mexican Hat wavelet”):

scal3 = ts.wavelet(settings = {'mother':'DOG'})
fig, ax = scal3.plot(title='CWT scalogram with DOG mother wavelet')

As for WWZ, note that, for computational efficiency, the time axis is coarse-grained by default to 50 time points, which explains in part the difference with the CWT scalogram.

If you need a custom axis, it (and other method-specific parameters) can also be passed via the settings dictionary:

tau = np.linspace(np.min(ts.time), np.max(ts.time), 60)
scal4 = ts.wavelet(method='wwz', settings={'tau':tau})
fig, ax = scal4.plot(title='WWZ scalogram with finer time axis')

wavelet_coherence(target_series, method='cwt', settings=None, freq=None, freq_kwargs=None, verbose=False, common_time_kwargs=None)[source]

Performs wavelet coherence analysis with the target timeseries

Parameters:

target_series (Series) – A pyleoclim Series object on which to perform the coherence analysis
method (str) – Possible methods {‘wwz’,’cwt’}. Default is ‘cwt’, which only works if the series share the same evenly-spaced time axis. ‘wwz’ is designed for unevenly-spaced data, but is far slower.
freq (str or array, optional) – Information to produce the frequency vector (highly consequential for the WWZ method) This can be ‘log’,’scale’, ‘nfft’, ‘lomb_scargle’, ‘welch’ or a NumPy array. If a string, will use make_freq_vector() with the specified frequency-generating method. If an array, this will be passed directly to the spectral method. If None (default), will use the ‘log’ method
freq_kwargs (dict) – Arguments for frequency vector
common_time_kwargs (dict) – Parameters for the method MultipleSeries.common_time(). Will use interpolation by default.
settings (dict) – Arguments for the specific wavelet method (e.g. decay constant for WWZ, mother wavelet for CWT) and common properties like standardize, detrend, gaussianize, pad, etc.
verbose (bool) – If True, will print warning messages, if any

Returns:

coh

Return type:

pyleoclim.core.coherence.Coherence

References

Grinsted, A., Moore, J. C. & Jevrejeva, S. Application of the cross wavelet transform and wavelet coherence to geophysical time series. Nonlin. Processes Geophys. 11, 561–566 (2004).

See also

pyleoclim.utils.spectral.make_freq_vector: Functions to create the frequency vector
pyleoclim.utils.tsutils.detrend: Detrending function
pyleoclim.core.multipleseries.MultipleSeries.common_time: put timeseries on common time axis
pyleoclim.core.series.Series.wavelet: wavelet analysis
pyleoclim.utils.wavelet.wwz_coherence: coherence using the wwz method
pyleoclim.utils.wavelet.cwt_coherence: coherence using the cwt method

Examples

Calculate the wavelet coherence of NINO3 and All India Rainfall with default arguments:

ts_air = pyleo.utils.load_dataset('AIR')
ts_nino = pyleo.utils.load_dataset('NINO3')
coh = ts_air.wavelet_coherence(ts_nino)
coh.plot()

(<Figure size 1000x800 with 2 Axes>,
 <Axes: xlabel='Time [year C.E.]', ylabel='Scale [yrs]'>)

Note that in this example both timeseries area already on a common, evenly-spaced time axis. If they are not (either because the data are unevenly spaced, or because the time axes are different in some other way), an error will be raised. To circumvent this error, you can either put the series on a common time axis (e.g. using common_time()) prior to applying CWT, or you can use the Weighted Wavelet Z-transform (WWZ) instead, as it is designed for unevenly-spaced data. However, it is usually far slower:

coh_wwz = ts_air.wavelet_coherence(ts_nino, method = 'wwz')
coh_wwz.plot()

(<Figure size 1000x800 with 2 Axes>,
 <Axes: xlabel='Time [year C.E.]', ylabel='Scale [yrs]'>)

As with wavelet analysis, both CWT and WWZ admit optional arguments through settings. For instance, one can adjust the resolution of the time axis on which coherence is evaluated:

coh_wwz = ts_air.wavelet_coherence(ts_nino, method = 'wwz', settings = {'ntau':20})
coh_wwz.plot()

(<Figure size 1000x800 with 2 Axes>,
 <Axes: xlabel='Time [year C.E.]', ylabel='Scale [yrs]'>)

The frequency (scale) axis can also be customized, e.g. to focus on scales from 1 to 20y, with 24 scales:

coh = ts_air.wavelet_coherence(ts_nino, freq_kwargs={'fmin':1/20,'fmax':1,'nf':24})
coh.plot()

(<Figure size 1000x800 with 2 Axes>,
 <Axes: xlabel='Time [year C.E.]', ylabel='Scale [yrs]'>)

Significance is assessed similarly to PSD or Scalogram objects:

cwt_sig = coh.signif_test(number=20, qs=[.9,.95]) # specifiying 2 significance thresholds does not take any more time.
# by default, the plot function will look for the closest quantile to 0.95, but it is easy to adjust:
cwt_sig.plot(signif_thresh = 0.9)

(<Figure size 1000x800 with 2 Axes>,
 <Axes: title={'center': 'Lines:90% threshold'}, xlabel='Time [year C.E.]', ylabel='Scale [yrs]'>)

Another plotting option, dashboard(), allows to visualize both timeseries as well as the wavelet transform coherency (WTC), which quantifies where two timeseries exhibit similar behavior in time-frequency space, and the cross-wavelet transform (XWT), which indicates regions of high common power.

cwt_sig.dashboard()

(<Figure size 900x1200 with 6 Axes>,
 {'ts1': <Axes: ylabel='AIR [mm/month]'>,
  'ts2': <Axes: xlabel='Time [year C.E.]', ylabel='NINO3 [$^{\\circ}$C]'>,
  'wtc': <Axes: ylabel='Scale [yrs]'>,
  'xwt': <Axes: xlabel='Time [year C.E.]', ylabel='Scale [yrs]'>})

Note: this design balances many considerations, and is not easily customizable.

GeoSeries (pyleoclim.GeoSeries)

class pyleoclim.core.geoseries.GeoSeries(time, value, lat, lon, elevation=None, time_unit=None, time_name=None, value_name=None, value_unit=None, label=None, importedFrom=None, archiveType=None, control_archiveType=False, sensorType=None, observationType=None, log=None, keep_log=False, verbose=True, depth=None, depth_name=None, depth_unit=None, sort_ts='ascending', dropna=True, clean_ts=False, auto_time_params=None)[source]

The GeoSeries class is a child of the Series class, and requires geolocation information (latitude, longitude). Elevation is optional, but can be used in mapping, if present. The class also allows for ancillary data and metadata, detailed below.

Parameters:

time (list or numpy.array) – independent variable (t)
value (list of numpy.array) – values of the dependent variable (y)
lat (float) – latitude N in decimal degrees. Must be in the range [-90;+90]
lon (float) – longitude East in decimal degrees. Must be in the range [-180;+360] No conversion is applied as mapping utilities convert to [-180,+180] internally
elevation (float) – elevation of the sample, in meters above sea level. Negative numbers indicate depth below global mean sea level, therefore.
time_unit (string) – Units for the time vector (e.g., ‘years’). Default is None
time_name (string) – Name of the time vector (e.g., ‘Time’,’Age’). Default is None. This is used to label the time axis on plots
value_name (string) – Name of the value vector (e.g., ‘temperature’) Default is None
value_unit (string) – Units for the value vector (e.g., ‘deg C’) Default is None
label (string) – Name of the time series (e.g., ‘Nino 3.4’) Default is None
log (dict) – Dictionary of tuples documentating the various transformations applied to the object
keep_log (bool) – Whether to keep a log of applied transformations. False by default
importedFrom (string) – source of the dataset. If it came from a LiPD file, this could be the datasetID property
archiveType (string) – climate archive, one of ‘Borehole’, ‘Coral’, ‘FluvialSediment’, ‘GlacierIce’, ‘GroundIce’, ‘LakeSediment’, ‘MarineSediment’, ‘Midden’, ‘MolluskShell’, ‘Peat’, ‘Sclerosponge’, ‘Shoreline’, ‘Speleothem’, ‘TerrestrialSediment’, ‘Wood’ Reference: https://lipdverse.org/vocabulary/archivetype/
control_archiveType ([True, False]) – Whether to standardize the name of the archiveType agains the vocabulary from: https://lipdverse.org/vocabulary/paleodata_proxy/. If set to True, will only allow for these terms and automatically convert known synonyms to the standardized name. Only standardized variable names will be automatically assigned a color scheme. Default is False.
sensorType (string) – sensor, e.g. a paleoclimate proxy sensor. This property can be used to differentiate between species of foraminifera
observationType (string) – observation type, e.g. a proxy observation. See https://lipdverse.org/vocabulary/paleodata_proxy/. Note: this is preferred terminology but not enforced
depth (array) – depth at which the values were collected
depth_name (string) – name of the field, e.g. ‘mid-depth’, ‘top-depth’, etc
depth_unit (string) – units of the depth axis, e.g. ‘cm’
dropna (bool) – Whether to drop NaNs from the series to prevent downstream functions from choking on them defaults to True
sort_ts (str) – Direction of sorting over the time coordinate; ‘ascending’ or ‘descending’ Defaults to ‘ascending’
verbose (bool) – If True, will print warning messages if there is any
clean_ts (boolean flag) – set to True to remove the NaNs and make time axis strictly prograde with duplicated timestamps reduced by averaging the values Default is None (marked for deprecation)
auto_time_params (bool,) – If True, uses tsbase.disambiguate_time_metadata to ensure that time_name and time_unit are usable by Pyleoclim. This may override the provided metadata. If False, the provided time_name and time_unit are used. This may break some functionalities (e.g. common_time and convert_time_unit), so use at your own risk. If not provided, code will set to True for internal consistency.

Examples

Import the EPICA Dome C deuterium record and display a quick synopsis:

ts = pyleo.utils.datasets.load_dataset('EDC-dD')
ts_interp = ts.convert_time_unit('kyr BP').interp(step=.5) # interpolate for a faster result
fig, ax = ts_interp.dashboard()

dashboard(figsize=[11, 8], gs=None, plt_kwargs=None, histplt_kwargs=None, spectral_kwargs=None, spectralsignif_kwargs=None, spectralfig_kwargs=None, map_kwargs=None, hue='archiveType', marker='archiveType', size=None, scatter_kwargs=None, gridspec_kwargs=None, savefig_settings=None)[source]

Create a dashboard of plots for the GeoSeries object

Parameters:

figsize (list or tuple, optional) – Figure size. The default is [11,8].
gs (matplotlib.gridspec object, optional) – Requires at least two rows and 4 columns. - top row, left: timeseries - top row, right: histogram - bottom left: map - bottom right: PSD See [matplotlib.gridspec.GridSpec](https://matplotlib.org/stable/tutorials/intermediate/gridspec.html) for details.
plt_kwargs (dict, optional) – Optional arguments for the timeseries plot. See Series.plot() or EnsembleSeries.plot_envelope(). The default is None.
histplt_kwargs (dict, optional) – Optional arguments for the distribution plot. See Series.histplot() or EnsembleSeries.plot_distplot(). The default is None.
spectral_kwargs (dict, optional) – Optional arguments for the spectral method. Default is to use Lomb-Scargle method. See Series.spectral() or EnsembleSeries.spectral(). The default is None.
spectralsignif_kwargs (dict, optional) – Optional arguments to estimate the significance of the power spectrum. See PSD.signif_test. Note that we currently do not support significance testing for ensembles. The default is None.
spectralfig_kwargs (dict, optional) – Optional arguments for the power spectrum figure. See PSD.plot() or MultiplePSD.plot_envelope(). The default is None.
map_kwargs (dict, optional) – Optional arguments for map configuration - projection: str; Optional value for map projection. Default ‘auto’. - proj_default: bool - lakes, land, ocean, rivers, borders, coastline, background: bool or dict; - lgd_kwargs: dict; Optional values for how the map legend is configured - gridspec_kwargs: dict; Optional values for adjusting the arrangement of the colorbar, map and legend in the map subplot - legend: bool; Whether to draw a legend on the figure. Default is True - colorbar: bool; Whether to draw a colorbar on the figure if the data associated with hue are numeric. Default is True The default is None.
hue (str, optional) – Variable associated with color coding for points plotted on map. May correspond to a continuous or categorical variable. The default is ‘archiveType’.
size (str, optional) – Variable associated with size. Must correspond to a continuous numeric variable. The default is None.
marker (string, optional) – Grouping variable that will produce points with different markers. Can have a numeric dtype but will always be treated as categorical. The default is ‘archiveType’.
scatter_kwargs (dict, optional) – Optional arguments configuring how data are plotted on a map. See description of scatter_kwargs in pyleoclim.utils.mapping.scatter_map
gridspec_kwargs (dict, optional) – Optional dictionary for configuring dashboard layout using gridspec For information about Gridspec configuration, refer to `Matplotlib documentation <https://matplotlib.org/3.5.0/api/_as_gen/matplotlib.gridspec.GridSpec.html#matplotlib.gridspec.GridSpec>_. The default is None.
savefig_settings (dict, optional) –
the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existed or non-existed path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”}.
The default is None.

Returns:

fig (matplotlib.figure) – The figure
ax (dict) – dictionary of matplotlib ax

See also

pyleoclim.core.series.Series.plot: plot a timeseries
pyleoclim.core.ensembleseries.EnsembleSeries.plot_envelope: Envelope plots for an ensemble
pyleoclim.core.series.Series.histplot: plot a distribution of the timeseries
pyleoclim.core.ensembleseries.EnsembleSeries.histplot: plot a distribution of the timeseries across ensembles
pyleoclim.core.series.Series.spectral: spectral analysis method.
pyleoclim.core.multipleseries.MultipleSeries.spectral: spectral analysis method for multiple series.
pyleoclim.core.psds.PSD.signif_test: significance test for timeseries analysis
pyleoclim.core.psds.PSD.plot: plot power spectrum
pyleoclim.core.psds.MulitplePSD.plot: plot envelope of power spectrum
pyleoclim.core.geoseries.GeoSeries.map: map location of dataset
pyleoclim.utils.mapping.scatter_map: Underlying mapping function for Pyleoclim

Examples

ts = pyleo.utils.datasets.load_dataset('EDC-dD')
ts_interp = ts.convert_time_unit('kyr BP').interp(step=.5) # interpolate for a faster result
fig, ax = ts_interp.dashboard()

classmethod from_json(path)[source]

Creates a pyleoclim.Series from a JSON file

The keys in the JSON file must correspond to the parameter associated with a GeoSeries object

Parameters:: path (str) – Path to the JSON file
Returns:: ts – A Pyleoclim Series object.
Return type:: pyleoclim.core.series.Series

map(projection='Orthographic', proj_default=True, background=True, borders=False, coastline=True, rivers=False, lakes=False, ocean=True, land=True, fig=None, gridspec_slot=None, title=None, figsize=None, marker='archiveType', hue='archiveType', size=None, edgecolor='w', markersize=None, scatter_kwargs=None, cmap=None, colorbar=False, gridspec_kwargs=None, legend=True, lgd_kwargs=None, savefig_settings=None)[source]

Map the location of the record

Parameters:

projection (str, optional) – The projection to use. The default is ‘Orthographic’.
proj_default (bool; {True, False}, optional) – Whether to use the Pyleoclim defaults for each projection type. The default is True.
background (bool, optional) – If True, uses a shaded relief background (only one available in Cartopy) Default is on (True).
borders (bool or dict, optional) – Draws the countries border. If a dictionary of formatting arguments is supplied (e.g. linewidth, alpha), will draw according to specifications. Defaults is off (False).
coastline (bool or dict, optional) – Draws the coastline. If a dictionary of formatting arguments is supplied (e.g. linewidth, alpha), will draw according to specifications. Defaults is on (True).
land (bool or dict, optional) – Colors land masses. If a dictionary of formatting arguments is supplied (e.g. color, alpha), will draw according to specifications. Default is off (True). Overriden if background=True.
ocean (bool or dict, optional) – Colors oceans. If a dictionary of formatting arguments is supplied (e.g. color, alpha), will draw according to specifications. Default is on (True). Overriden if background=True.
rivers (bool or dict, optional) – Draws major rivers. If a dictionary of formatting arguments is supplied (e.g. linewidth, alpha), will draw according to specifications. Default is off (False).
lakes (bool or dict, optional) – Draws major lakes. If a dictionary of formatting arguments is supplied (e.g. color, alpha), will draw according to specifications. Default is off (False).
fig (matplotlib.pyplot.figure, optional) – See matplotlib.pyplot.figure <https://matplotlib.org/3.5.0/api/_as_gen/matplotlib.pyplot.figure.html#matplotlib-pyplot-figure>_. The default is None.
gridspec_slot (Gridspec slot, optional) – If generating a map for a multi-plot, pass a gridspec slot. The default is None.
title (string) – Title for the plot. Default is None
size (string, optional) – Grouping variable that will produce points with different sizes. Expects to be associated with numeric values (e.g., ‘elevation’ has values). Any data without a value for the size variable will be filtered out. The default is None.
edgecolor (color (string) or list of rgba tuples, optional) – Color of marker edge. The default is ‘w’.
figsize (list or tuple, optional) – The size of the figure. The default is None.
marker (str, optional) – The marker type for each archive. The default is None. Uses plot_default
hue (str, optional) – Variable associated with color coding. The default is None. Uses plot_default.
markersize (float, optional) – Size of the marker. The default is None.
scatter_kwargs (dict, optional) – Parameters for the scatter plot. The default is None.
cmap (string or list, optional) – Matplotlib supported colormap id or list of colors for creating a colormap. See choosing a matplotlib colormap. The default is None.
colorbar (bool, optional) – Whether the draw a colorbar on the figure if the data associated with hue are numeric. Default is True.
gridspec_kwargs (dict, optional) – Function assumes the possibility of a colorbar, map, and legend. A list of floats associated with the keyword width_ratios will assume the first (index=0) is the relative width of the colorbar, the second to last (index = -2) is the relative width of the map, and the last (index = -1) is the relative width of the area for the legend. For information about Gridspec configuration, refer to Matplotlib documentation. The default is None.
legend (bool; {True, False}, optional) – Whether to plot the legend. The default is True.
lgd_kwargs (dict, optional) – Arguments for the legend. The default is None.
savefig_settings (dict, optional) –
the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existed or non-existed path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”}. The default is None.

Returns:

res

Return type:

fig,ax_d

See also

pyleoclim.utils.mapping.scatter_map: Underlying mapping function for Pyleoclim

Examples

ts = pyleo.utils.datasets.load_dataset('EDC-dD')
fig, ax = ts.map()

By default, the figure has no title. For a title built from the available labels:

fig, ax = ts.map(title=True)

For a custom title, and custom projection:

fig, ax = ts.map(title='Insert title here', projection='RotatedPole',
                 proj_default={'pole_longitude':0.0, 'pole_latitude':-90.0, 'central_rotated_longitude':45.0})

map_neighbors(mgs, radius=3000, projection='Orthographic', proj_default=True, background=True, borders=False, rivers=False, lakes=False, ocean=True, land=True, fig=None, gridspec_slot=None, title=None, figsize=None, marker='archiveType', hue='archiveType', size=None, edgecolor=None, markersize=None, scatter_kwargs=None, cmap=None, colorbar=False, gridspec_kwargs=None, legend=True, lgd_kwargs=None, savefig_settings=None)[source]

Map all records within a given radius of the object

Parameters:

mgs (MultipleGeoSeries) – object containing the series to be considered as neighbors
radius (float) – search radius for the record, in km. Default is 3000.
projection (str, optional) – The projection to use. The default is ‘Orthographic’.
proj_default (bool; {True, False}, optional) – Whether to use the Pyleoclim defaults for each projection type. The default is True.
background (bool; {True, False}, optional) – Whether to use a background. The default is True.
borders (bool; {True, False}, optional) – Draw borders. The default is False.
rivers (bool; {True, False}, optional) – Draw rivers. The default is False.
lakes (bool; {True, False}, optional) – Draw lakes. The default is False.
figsize (list or tuple, optional) – The size of the figure. The default is None.
marker (str, optional) – The marker type for each archive. The default is None. Uses plot_default
hue (str, optional) – Variable associated with color coding. The default is None. Uses plot_default.
markersize (float, optional) – Size of the marker. The default is None.
scatter_kwargs (dict, optional) – Parameters for the scatter plot. The default is None.
legend (bool; {True, False}, optional) – Whether to plot the legend. The default is True.
lgd_kwargs (dict, optional) – Arguments for the legend. The default is None.
savefig_settings (dict, optional) –
the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existed or non-existed path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”}. The default is None.
title (bool or str) – the title for the figure. If True or None, made automatically from the objects’ labels. Set to False for an empty title.

Returns:

res

Return type:

fig,ax_d

See also

pyleoclim.utils.mapping.map: Underlying mapping function for Pyleoclim

Examples

from pylipd.utils.dataset import load_dir
lipd = load_dir(name='Pages2k')
df = lipd.get_timeseries_essentials()
dfs = df.query("archiveType in ('Wood','Documents','Coral','Lake sediment') and paleoData_variableName not in ('year')")

# place in a MultipleGeoSeries object
ts_list = []
for _, row in dfs.iterrows():
    ts_list.append(pyleo.GeoSeries(time=row['time_values'],value=row['paleoData_values'],
                                    time_name=row['time_variableName'],value_name=row['paleoData_variableName'],
                                    time_unit=row['time_units'], value_unit=row['paleoData_units'],
                                    lat = row['geo_meanLat'], lon = row['geo_meanLon'],
                                    archiveType = row['archiveType'], verbose = False,
                                    label=row['dataSetName']+'_'+row['paleoData_variableName']))

mgs = pyleo.MultipleGeoSeries(series_list=ts_list,time_unit='years AD',label='Euro2k')
gs = ts_list[6] # extract one record as the target one
gs.map_neighbors(mgs, radius=4000)

Loading 16 LiPD files

Loaded..

(<Figure size 1800x700 with 2 Axes>,
 {'map': <GeoAxes: xlabel='lon', ylabel='lat'>, 'leg': <Axes: >})

By default, the figure has no title. For a title built from the available labels:

gs.map_neighbors(mgs, radius=4000, title=True)

(<Figure size 1800x700 with 2 Axes>,
 {'map': <GeoAxes: title={'center': 'Euro2k neighbors for Ocn-RedSea.Felis.2000_d18O within 4000 km'}, xlabel='lon', ylabel='lat'>,
  'leg': <Axes: >})

For a custom title:

gs.map_neighbors(mgs, radius=4000, title='Insert title here')

(<Figure size 1800x700 with 2 Axes>,
 {'map': <GeoAxes: title={'center': 'Insert title here'}, xlabel='lon', ylabel='lat'>,
  'leg': <Axes: >})

resample(rule, keep_log=False, **kwargs)[source]

Analogue to pandas.Series.resample.

This is a convenience method: doing

ser.resample(‘YS’).mean()

will do the same thing as

ser.pandas_method(lambda x: x.resample(‘YS’).mean())

but will also accept some extra resampling rules, such as ‘Ga’ (see below).

NOTE: this feature may break until [this pandas bug](https://github.com/pandas-dev/pandas/issues/57427) is fixed.

Parameters:

rule (str) –
The offset string or object representing target conversion. Can also accept pyleoclim units, such as ‘ka’ (1000 years), ‘Ma’ (1 million years), and ‘Ga’ (1 billion years).

Check the [pandas resample docs](https://pandas.pydata.org/docs/dev/reference/api/pandas.DataFrame.resample.html) for more details.
kwargs (dict) – Any other arguments which will be passed to pandas.Series.resample.

Returns:

Resampler object, not meant to be used to directly. Instead, an aggregation should be called on it, see examples below.

Return type:

SeriesResampler

Examples

ts = pyleo.utils.load_dataset('EDC-dD').convert_time_unit('ky BP')
ts5k = ts.resample('1ka').mean()
fig, ax = ts.plot()
ts5k.plot(ax=ax,color='C1')

<Axes: xlabel='Age [ka]', ylabel='$\\delta \\mathrm{D}$ [‰]'>

segment(factor=10, verbose=False)[source]

Gap detection

This function segments a timeseries into n parts following a gap- detection algorithm. The rule of gap detection is very simple:: we define the intervals between time points as dts, then if dts[i] is larger than factor * dts[i-1], we think that the change of dts (or the gradient) is too large, and we regard it as a breaking point and divide the time series into two segments here

Parameters:

factor (float) – The factor that adjusts the threshold for gap detection
verbose (bool) – If True, will print warning messages if there is any

Returns:

res – If gaps were detected, returns the segments in a MultipleGeoSeries object, else, returns the original timeseries.

Return type:

MultiplegGeoSeries or GeoSeries

Examples

import numpy as np
gs = pyleo.utils.datasets.load_dataset('EDC-dD')
gs.value[4000:5000] = np.nan # cut a large gap in the middle
mgs = gs.segment()
mgs.plot()

(<Figure size 1000x400 with 1 Axes>,
 <Axes: xlabel='Age [y BP]', ylabel='$\\delta \\mathrm{D}$ [‰]'>)

MultipleSeries (pyleoclim.MultipleSeries)

class pyleoclim.core.multipleseries.MultipleSeries(series_list, time_unit=None, label=None, name=None)[source]

MultipleSeries object.

This object handles a collection of the type Series and can be created from a list of such objects. MultipleSeries should be used when the need to run analysis on multiple records arises, such as running principal component analysis. Some of the methods automatically transform the time axis prior to analysis to ensure consistency.

Parameters:

series_list (list) – a list of pyleoclim.Series objects
time_unit (str) – The target time unit for every series in the list. If None, then no conversion will be applied; Otherwise, the time unit of every series in the list will be converted to the target.
label (str) – label of the collection of timeseries (e.g. ‘PAGES 2k ice cores’)

Examples

soi = pyleo.utils.load_dataset('SOI')
nino = pyleo.utils.load_dataset('NINO3')
ms = soi & nino
ms.label = 'ENSO'
ms

            Southern Oscillation Index  NINO3 SST
datetime                                         
1866-01-01                       -0.62        NaN
1866-02-01                       -0.12        NaN
1866-03-01                       -0.62        NaN
1866-04-01                       -0.65        NaN
1866-05-01                        0.04        NaN
...                                ...        ...
2024-10-01                        0.35        NaN
2024-11-01                        0.61        NaN
2024-12-01                        1.09        NaN
2025-01-01                        0.26        NaN
2025-02-01                        0.21        NaN

[3506 rows x 2 columns]

append(ts, inplace=False)[source]

Append timeseries ts to MultipleSeries object

Parameters:: ts (pyleoclim.Series) – The pyleoclim Series object to be appended to the MultipleSeries object
Returns:: ms – The augmented object, comprising the old one plus ts
Return type:: MultipleSeries

See also

pyleoclim.core.series.Series: A Pyleoclim Series object

Examples

soi = pyleo.utils.load_dataset('SOI')
NINO3 = pyleo.utils.load_dataset('NINO3')
ms = pyleo.MultipleSeries([soi], label = 'ENSO')
ms.append(NINO3)

The two series have different lengths, left: 1910 vs right: 1596
Metadata are different:
value_unit property -- left: mb, right: $^{\circ}$C
value_name property -- left: SOI, right: NINO3
label property -- left: Southern Oscillation Index, right: NINO3 SST

            Southern Oscillation Index  NINO3 SST
datetime                                         
1866-01-01                       -0.62        NaN
1866-02-01                       -0.12        NaN
1866-03-01                       -0.62        NaN
1866-04-01                       -0.65        NaN
1866-05-01                        0.04        NaN
...                                ...        ...
2024-10-01                        0.35        NaN
2024-11-01                        0.61        NaN
2024-12-01                        1.09        NaN
2025-01-01                        0.26        NaN
2025-02-01                        0.21        NaN

[3506 rows x 2 columns]

bin(**kwargs)[source]

Aligns the time axes of a MultipleSeries object, via binning.

This is critical for workflows that need to assume a common time axis for the group of series under consideration.

The common time axis is characterized by the following parameters:

start : the latest start date of the bunch (maximin of the minima)

stop : the earliest stop date of the bunch (minimum of the maxima)

step : The representative spacing between consecutive values (mean of the median spacings)

This is a special case of the common_time function.

Parameters:: kwargs (dict) – Arguments for the binning function. See pyleoclim.utils.tsutils.bin
Returns:: ms – The MultipleSeries objects with all series aligned to the same time axis.
Return type:: MultipleSeries

See also

pyleoclim.core.multipleseries.MultipleSeries.common_time: Base function on which this operates
pyleoclim.utils.tsutils.bin: Underlying binning function
pyleoclim.core.series.Series.bin: Bin function for Series object

Examples

soi = pyleo.utils.load_dataset('SOI')
nino = pyleo.utils.load_dataset('NINO3')
ms = soi & nino
msbin = ms.bin()

common_time(method='interp', step=None, start=None, stop=None, step_style=None, time_axis=None, **kwargs)[source]

Aligns the time axes of a MultipleSeries object

The alignment is achieved via binning, interpolation, or Gaussian kernel. Alignment is critical for workflows that need to assume a common time axis for the group of series under consideration.

The common time axis is characterized by the following parameters:

start : the latest start date of the bunch (maximun of the minima)

stop : the earliest stop date of the bunch (minimum of the maxima)

step : The representative spacing between consecutive values

Optional arguments for binning, Gaussian kernel (gkernel) interpolation are those of the underling functions.

If any of the time axes are retrograde, this step makes them prograde.

Parameters:

method (string; {'bin','interp','gkernel'}) – either ‘bin’, ‘interp’ [default] or ‘gkernel’
step (float) – common step for all time axes. Default is None and inferred from the timeseries spacing
start (float) – starting point of the common time axis. Default is None and inferred as the max of the min of the time axes for the timeseries.
stop (float) – end point of the common time axis. Default is None and inferred as the min of the max of the time axes for the timeseries.
step_style (string; {'median', 'mean', 'mode', 'max'}) – Method to obtain a representative step among all Series (using tsutils.increments). Default value is None, so that it will be chosen according to the method: ‘max’ for bin and gkernel, ‘mean’ for interp.
time_axis (array) – Time axis onto which all the series will be aligned. Will override step,start,stop, and step_style if they are passed.
kwargs (dict) – keyword arguments (dictionary) of the bin, gkernel or interp methods

Returns:

ms – The MultipleSeries objects with all series aligned to the same time axis.

Return type:

MultipleSeries

Notes

start, stop, step, and step_style are interpreted differently depending on the method used. Interp uses these to specify the time_axis onto which interpolation will be applied. Bin and gkernel use these to specify the bin_edges which define the “buckets” used for the respective methods.

See also

pyleoclim.utils.tsutils.bin: put timeseries values into bins of equal size (possibly leaving NaNs in).
pyleoclim.utils.tsutils.gkernel: coarse-graining using a Gaussian kernel
pyleoclim.utils.tsutils.interp: interpolation onto a regular grid (default = linear interpolation)
pyleoclim.utils.tsutils.increments: infer grid properties

Examples

from pyleoclim.utils.tsmodel import colored_noise

# create 2 incompletely sampled series
ns = 2 ; nt = 200; n_del = 20
serieslist = []

for j in range(ns):
    t = np.arange(nt)
    v = colored_noise(alpha=1, t=t)
    deleted_idx = np.random.choice(range(np.size(t)), n_del, replace=False)
    tu =  np.delete(t, deleted_idx)
    vu =  np.delete(v, deleted_idx)
    ts = pyleo.Series(time = tu, value = vu, label = 'series {}'.format(j+1), verbose=False)
    serieslist.append(ts)

# create MS object from the list
ms = pyleo.MultipleSeries(serieslist)

import matplotlib.pyplot as plt
fig, ax = plt.subplots(2,2,sharex=True,sharey=True, figsize=(10,8))
ax = ax.flatten()
# apply common_time with default parameters
msc = ms.common_time()
msc.plot(title='linear interpolation',ax=ax[0], legend=False)

# apply common_time with binning
msc = ms.common_time(method='bin')
msc.plot(title='Binning',ax=ax[1], legend=False)

# apply common_time with gkernel
msc = ms.common_time(method='gkernel')
msc.plot(title=r'Gaussian kernel ($h=3$)',ax=ax[2],legend=False)

# apply common_time with gkernel and a large bandwidth
msc = ms.common_time(method='gkernel', h=.5)
msc.plot(title=r'Gaussian kernel ($h=.5$)',ax=ax[3],legend=False)
fig.tight_layout()
# Optional close fig after plotting

convert_time_unit(time_unit=None)[source]

Convert the time units of the object

Parameters:

time_unit (str) –

the target time unit, possible input: {

’year’, ‘years’, ‘yr’, ‘yrs’, ‘y BP’, ‘yr BP’, ‘yrs BP’, ‘year BP’, ‘years BP’, ‘ky BP’, ‘kyr BP’, ‘kyrs BP’, ‘ka BP’, ‘ka’, ‘my BP’, ‘myr BP’, ‘myrs BP’, ‘ma BP’, ‘ma’,

}

Examples

soi = pyleo.utils.load_dataset('SOI')
nino = pyleo.utils.load_dataset('NINO3')
ms = soi & nino
new_ms = ms.convert_time_unit('yr BP')
print('Original timeseries:')
print('time unit:', ms.time_unit)
print()
print('Converted timeseries:')
print('time unit:', new_ms.time_unit)

Original timeseries:
time unit: None

Converted timeseries:
time unit: yr BP

copy()[source]

Copy the object

Returns:: ms – The copied version of the pyleoclim.MultipleSeries object
Return type:: MultipleSeries

Examples

soi = pyleo.utils.load_dataset('SOI')
nino = pyleo.utils.load_dataset('NINO3')
ms = soi & nino
ms_copy = ms.copy()

correlation(target=None, timespan=None, alpha=0.05, statistic='pearsonr', method='phaseran', number=1000, settings=None, fdr_kwargs=None, common_time_kwargs=None, mute_pbar=False, seed=None)[source]

Calculate the correlation between a MultipleSeries and a target Series

This function recursively applies Series.correlation() to members of the MultipleSeries object.

The significance of the correlation is assessed using one of the following methods:

‘ttest’: T-test adjusted for effective sample size, see [1]
‘ar1sim’: AR(1) modeling of x and y (Monte Carlo method)
‘CN’: colored noise (power-law spectrum) modeling of x and y (Monte Carlo method)
‘phaseran’: phase randomization of original inputs. (Monte Carlo method, default), see [2]
‘built-in’: uses built-in method from scipy (function of the statistic used)

Note: Up to version v0.14.0. ar1sim was called “isopersistent”, phaseran was called “isospectral”

The T-test is a parametric test, hence computationally cheap, but can only be performed in ideal circumstances. The others are non-parametric, but their computational requirements scale with the number of simulations.

The False Discovery Rate method is applied to the assessment of significance when plotting the result.

If the computation fails, a diagnostic message returns the index and label of the incriminated series.

Parameters:

target (pyleoclim.Series, optional) – The Series against which to take the correlation. If the target Series is not specified, then the 1st member of MultipleSeries will be used as the target
timespan (tuple, optional) – The time interval over which to perform the calculation
alpha (float) – The significance level (0.05 by default)
statistic (str) –

statistic being evaluated. Can use any of the SciPy-supported ones:
https://docs.scipy.org/doc/scipy/reference/stats.html#association-correlation-tests Currently supported: [‘pearsonr’,’spearmanr’,’pointbiserialr’,’kendalltau’,’weightedtau’]

Default: ‘pearsonr’.
method (str, {'ttest','built-in','ar1sim','phaseran','CN'}) – method for significance testing. Default is ‘phaseran’ - ‘ttest’ implements the T-test with degrees of freedom adjusted for autocorrelation, as done in [1] - ‘built-in’ uses the p-value that ships with the SciPy function. - ‘ar1sim’ (formerly ‘isopersistent’) tests against an ensemble of AR(1) series fitted to the originals - ‘CN’ tests against an ensemble of colored noise series (power-law spectra) fitted to the originals - ‘phaseran’ (formerly ‘isospectral’) tests against phase-randomized surrogates (aka the method of Ebisuzaki [2]) The old options ‘isopersistent’ and ‘isospectral’ still work, but trigger a deprecation warning. Note that ‘weightedtau’ does not have a known distribution, so the ‘built-in’ method returns an error in that case.
settings (dict) – Parameters for the correlation function (per scipy)
number (int) – the number of simulations (default: 1000)
fdr_kwargs (dict) – Parameters for the FDR function
common_time_kwargs (dict) – Parameters for the method MultipleSeries.common_time()
mute_pbar (bool; {True,False}) – If True, the progressbar will be muted. Default is False.
seed (float or int) – random seed for isopersistent and isospectral methods

Returns:

corr – the result object, containing the following: - statistic r (array of real numbers) - p-values pval (array of real numbers) - signif (array of booleans) - alpha (significance level)

Return type:

CorrEns

See also

pyleoclim.core.series.Series.correlation: Series-level correlation
pyleoclim.utils.correlation.association: SciPy measures of association between variables
pyleoclim.series.surrogates: parametric and non-parametric surrogates of any Series object
pyleoclim.multipleseries.common_time: Aligning time axes
pyleoclim.utils.correlation.fdr: FDR function
pyleoclim.core.correns.CorrEns: the correlation ensemble object

Examples

from pyleoclim.utils.tsmodel import colored_noise
import numpy as np

nt = 100
t0 = np.arange(nt)
v0 = colored_noise(alpha=1, t=t0)
noise = np.random.normal(loc=0, scale=1, size=nt)

ts0 = pyleo.Series(time=t0, value=v0, verbose=False)
ts1 = pyleo.Series(time=t0, value=v0+noise, verbose=False)
ts2 = pyleo.Series(time=t0, value=v0+2*noise, verbose=False)
ts3 = pyleo.Series(time=t0, value=v0+1/2*noise, verbose=False)

ts_list = [ts1, ts2, ts3]

ms = pyleo.MultipleSeries(ts_list)
ts_target = ts0

Correlation between the MultipleSeries object and a target Series. We also set an arbitrary random seed to ensure reproducibility:

corr_res = ms.correlation(ts_target, number=20, seed=2333)
print(corr_res)

Looping over 3 Series in collection

  correlation  p-value    signif. w/o FDR (α: 0.05)    signif. w/ FDR (α: 0.05)
-------------  ---------  ---------------------------  --------------------------
     0.720717  < 1e-6     True                         True
     0.460458  < 1e-6     True                         True
     0.901376  < 1e-6     True                         True
Ensemble size: 3

Correlation among the series of the MultipleSeries object

corr_res = ms.correlation(number=20, seed=2333)
print(corr_res)

Looping over 3 Series in collection

  correlation  p-value    signif. w/o FDR (α: 0.05)    signif. w/ FDR (α: 0.05)
-------------  ---------  ---------------------------  --------------------------
     1         < 1e-6     True                         True
     0.947227  < 1e-6     True                         True
     0.949831  < 1e-6     True                         True
Ensemble size: 3

detrend(method='emd', **kwargs)[source]

Detrend timeseries

Parameters:

method (str, optional) –
The method for detrending. The default is ‘emd’. Options include:
- linear: the result of a linear least-squares fit to y is subtracted from y.
- constant: only the mean of data is subtrated.
- ’savitzky-golay’, y is filtered using the Savitzky-Golay filters and the resulting filtered series is subtracted from y.
- ’emd’ (default): Empirical mode decomposition. The last mode is assumed to be the trend and removed from the series
**kwargs (dict) – Relevant arguments for each of the methods.

Returns:

ms – The detrended timeseries

Return type:

MultipleSeries

See also

pyleoclim.core.series.Series.detrend: Detrending for a single series
pyleoclim.utils.tsutils.detrend: Detrending function

equal_lengths()[source]

Test whether all series in object have equal length

Returns:

flag (bool) – Whether or not the Series in the pyleo.MultipleSeries object are of equal length
lengths (list) – List of the lengths of the series in object

See also

pyleoclim.core.multipleseries.MultipleSeries.common_time: Aligns the time axes of a MultipleSeries object

Examples

soi = pyleo.utils.load_dataset('SOI')
nino = pyleo.utils.load_dataset('NINO3')
ms = soi & nino
flag, lengths = ms.equal_lengths()
print(flag)

False

filter(cutoff_freq=None, cutoff_scale=None, method='butterworth', **kwargs)[source]

Filtering the timeseries in the MultipleSeries object

Parameters:

method (str; {'savitzky-golay', 'butterworth', 'firwin', 'lanczos'}) – The filtering method - ‘butterworth’: the Butterworth method (default) - ‘savitzky-golay’: the Savitzky-Golay method - ‘firwin’: FIR filter design using the window method, with default window as Hamming - ‘lanczos’: lowpass filter via Lanczos resampling
cutoff_freq (float or list) – The cutoff frequency only works with the Butterworth method. If a float, it is interpreted as a low-frequency cutoff (lowpass). If a list, it is interpreted as a frequency band (f1, f2), with f1 < f2 (bandpass).
cutoff_scale (float or list) – cutoff_freq = 1 / cutoff_scale The cutoff scale only works with the Butterworth method and when cutoff_freq is None. If a float, it is interpreted as a low-frequency (high-scale) cutoff (lowpass). If a list, it is interpreted as a frequency band (f1, f2), with f1 < f2 (bandpass).
kwargs (dict) – A dictionary of the keyword arguments for the filtering method, See pyleoclim.utils.filter.savitzky_golay, pyleoclim.utils.filter.butterworth, pyleoclim.utils.filter.firwin, and pyleoclim.utils.filter.lanczos for the details

Returns:

ms

Return type:

MultipleSeries

See also

pyleoclim.series.Series.filter: filtering for Series objects
pyleoclim.utils.filter.butterworth: Butterworth method
pyleoclim.utils.filter.savitzky_golay: Savitzky-Golay method
pyleoclim.utils.filter.firwin: FIR filter design using the window method
pyleoclim.utils.filter.lanczos: lowpass filter via Lanczos resampling

Examples

air = pyleo.utils.load_dataset('AIR')
nino = pyleo.utils.load_dataset('NINO3')
ms = air & nino
ms_filter = ms.filter(method='lanczos',cutoff_scale=20)

flip(axis='value')[source]

Flips the Series along one or both axes

Parameters:

axis (str, optional) – The axis along which the Series will be flipped. The default is ‘value’. Other acceptable options are ‘time’ or ‘both’.

Returns:

ms (MultipleSeries) – The flipped object

Examples
.. jupyter-execute:: – soi = pyleo.utils.load_dataset(‘SOI’) nino = pyleo.utils.load_dataset(‘NINO3’) ms = soi & nino ms.name = ‘ENSO’ fig, ax = ms.flip().stackplot()

Note that labels have been updated to reflect the flip

classmethod from_json(path)[source]

Creates a pyleoclim.MulitpleSeries from a JSON file

The keys in the JSON file must correspond to the parameter associated with MulitpleSeries and Series objects

Parameters:: path (str) – Path to the JSON file
Returns:: ts – A Pyleoclim MultipleSeries object.
Return type:: pyleoclim.core.series.MulitplesSeries

gkernel(**kwargs)[source]

Aligns the time axes of a MultipleSeries object, via Gaussian kernel.

This is critical for workflows that need to assume a common time axis for the group of series under consideration.

The common time axis is characterized by the following parameters:

start : the latest start date of the bunch (maximin of the minima)

stop : the earliest stop date of the bunch (minimum of the maxima)

step : The representative spacing between consecutive values (mean of the median spacings)

This is a special case of the common_time function.

Parameters:: kwargs (dict) – Arguments for gkernel. See pyleoclim.utils.tsutils.gkernel for details.
Returns:: ms – The MultipleSeries objects with all series aligned to the same time axis.
Return type:: MultipleSeries

See also

pyleoclim.core.multipleseries.MultipleSeries.common_time: Base function on which this operates
pyleoclim.utils.tsutils.gkernel: Underlying kernel module

Examples

soi = pyleo.utils.load_dataset('SOI')
nino = pyleo.utils.load_dataset('NINO3')
ms = soi & nino
msk = ms.gkernel()

increments(step_style='median', verbose=False)[source]

Extract grid properties (start, stop, step) of all the Series objects in a collection.

Parameters:

step_style (str; {'median','mean','mode','max'}) –
Method to obtain a representative step if x is not evenly spaced. Valid entries: ‘median’ [default], ‘mean’, ‘mode’ or ‘max’. The “mode” is the most frequent entry in a dataset, and may be a good choice if the timeseries is nearly equally spaced but for a few gaps.

”max” is a conservative choice, appropriate for binning methods and Gaussian kernel coarse-graining
verbose (bool) – If True, will print out warning messages when they appear

Returns:

increments – n x 3 array, where n is the number of series,

index 0 is the earliest time among all Series
index 1 is the latest time among all Series
index 2 is the step, chosen according to step_style

Return type:

numpy.array

See also

pyleoclim.utils.tsutils.increments: underlying array-level utility

Examples

soi = pyleo.utils.load_dataset('SOI')
nino = pyleo.utils.load_dataset('NINO3')
ms = soi & nino
increments = ms.increments()

interp(**kwargs)[source]

Aligns the time axes of a MultipleSeries object, via interpolation.

This is critical for workflows that need to assume a common time axis for the group of series under consideration.

The common time axis is characterized by the following parameters:

start : the latest start date of the bunch (maximin of the minima)

stop : the earliest stop date of the bunch (minimum of the maxima)

step : The representative spacing between consecutive values (mean of the median spacings)

This is a special case of the common_time function.

Parameters:: kwargs (keyword arguments (dictionary) for the interpolation method) –
Returns:: ms – The MultipleSeries objects with all series aligned to the same time axis.
Return type:: MultipleSeries

See also

pyleoclim.core.multipleseries.MultipleSeries.common_time: Base function on which this operates
pyleoclim.utils.tsutils.interp: Underlying interpolation function
pyleoclim.core.series.Series.interp: Interpolation function for Series object

Examples

soi = pyleo.utils.load_dataset('SOI')
nino = pyleo.utils.load_dataset('NINO3')
ms = soi & nino
msinterp = ms.interp()

pca(weights=None, name=None, missing='fill-em', tol_em=0.005, max_em_iter=100, **pca_kwargs)[source]

Principal Component Analysis (Empirical Orthogonal Functions)

Decomposition of MultipleSeries in terms of orthogonal basis functions. Tolerant to missing values, infilled by an EM algorithm.

Do make sure the time axes are aligned, however! (e.g. use common_time())

Algorithm from statsmodels: https://www.statsmodels.org/stable/generated/statsmodels.multivariate.pca.PCA.html

Parameters:

weights (ndarray, optional) – Series weights to use after transforming data according to standardize or demean when computing the principal components.
missing ({str, None}) –
Method for missing data. Choices are:
- ’drop-row’ - drop rows with missing values.
- ’drop-col’ - drop columns with missing values.
- ’drop-min’ - drop either rows or columns, choosing by data retention.
- ’fill-em’ - use EM algorithm to fill missing value. ncomp should be set to the number of factors required.
- None raises if data contains NaN values.
tol_em (float) – Tolerance to use when checking for convergence of the EM algorithm.
max_em_iter (int) – Maximum iterations for the EM algorithm.

Returns:

res – Resulting pyleoclim.MultivariateDecomp object

Return type:

MultivariateDecomp

See also

pyleoclim.utils.tsutils.eff_sample_size: Effective Sample Size of timeseries y
pyleoclim.core.multivardecomp.MultivariateDecomp: The spatial decomposition object
pyleoclim.core.mulitpleseries.MulitpleSeries.common_time: align time axes

Examples

soi = pyleo.utils.load_dataset('SOI')
nino = pyleo.utils.load_dataset('NINO3')
ms = (soi & nino).common_time()
ms.label = ms.series_list[0].label
res = ms.pca() # carry out PCA

fig1, ax1 = res.screeplot() # plot the eigenvalue spectrum
fig2, ax2 = res.modeplot() # plot the first mode

The provided eigenvalue array has only one dimension. UQ defaults to NB82

plot(figsize=[10, 4], marker=None, markersize=None, linestyle=None, linewidth=None, colors=None, cmap='tab10', norm=None, xlabel=None, ylabel=None, title=None, time_unit=None, legend=True, plot_kwargs=None, lgd_kwargs=None, savefig_settings=None, ax=None, invert_xaxis=False)[source]

Plot multiple timeseries on the same axis

Parameters:

figsize (list, optional) – Size of the figure. The default is [10, 4].
marker (str, optional) – Marker type. The default is None.
markersize (float, optional) – Marker size. The default is None.
linestyle (str, optional) – Line style. The default is None.
linewidth (float, optional) – The width of the line. The default is None.
colors (a list of, or one, Python supported color code (a string of hex code or a tuple of rgba values)) – Colors for plotting. If None, the plotting will cycle the ‘tab10’ colormap; if only one color is specified, then all curves will be plotted with that single color; if a list of colors are specified, then the plotting will cycle that color list.
cmap (str) – The colormap to use when “colors” is None.
norm (matplotlib.colors.Normalize) – The normalization for the colormap. If None, a linear normalization will be used.
xlabel (str, optional) – x-axis label. The default is None.
ylabel (str, optional) – y-axis label. The default is None.
title (str, optional) – Title. The default is None.
time_unit (str) –
the target time unit, possible input: {

’year’, ‘years’, ‘yr’, ‘yrs’, ‘y BP’, ‘yr BP’, ‘yrs BP’, ‘year BP’, ‘years BP’, ‘ky BP’, ‘kyr BP’, ‘kyrs BP’, ‘ka BP’, ‘ka’, ‘my BP’, ‘myr BP’, ‘myrs BP’, ‘ma BP’, ‘ma’,

} default is None, in which case the code picks the most common time unit in the collection. If no unambiguous winner can be found, the unit of the first series in the collection is used.
legend (bool, optional) – Whether the show the legend. The default is True.
plot_kwargs (dict, optional) – Plot parameters. The default is None.
lgd_kwargs (dict, optional) – Legend parameters. The default is None.
savefig_settings (dictionary, optional) –
the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existing or non-existing path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”} The default is None.
ax (matplotlib.ax, optional) – The matplotlib axis onto which to return the figure. The default is None.
invert_xaxis (bool, optional) – if True, the x-axis of the plot will be inverted

Returns:

fig (matplotlib.figure) – the figure object from matplotlib See [matplotlib.pyplot.figure](https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.figure.html) for details.
ax (matplotlib.axis) – the axis object from matplotlib See [matplotlib.axes](https://matplotlib.org/api/axes_api.html) for details.

See also

pyleoclim.utils.plotting.savefig: Saving figure in Pyleoclim

Examples

soi = pyleo.utils.load_dataset('SOI')
nino = pyleo.utils.load_dataset('NINO3')
ms = soi & nino
ms.name = 'ENSO'
fig, ax = ms.plot()

remove(label)[source]

Remove Series based on given label.

Modifies the MultipleSeries, does not return anything.

resolution(time_unit=None, verbose=True, statistic='median')[source]

Generate a MultipleResolution object

Increments are assigned to the preceding time value. E.g. for time_axis = [0,1,3], resolution.resolution = [1,2] resolution.time = [0,1]. Note that the MultipleResolution class requires a shared time unit. If the time_unit parameter is not passed, a time unit will be automatically determined.

Returns:

multipleresolution (pyleoclim.MultipleResolution) – MultipleResolution object
time_unit (str) – Time unit to convert objects to. See pyleo.Series.convert_time_unit for options.
verbose (bool) – Whether or not to print messages warning the user about automated decisions.
statistic (str; {‘median’,’mean’,None}) – If a recognized statistic is passed, this function will simply output that statistic applied to the resolution of each series in the MulitipleSeries object. Options are ‘mean’ or ‘median’. If statistic is None, then the function will return a new MultipleResolution class with plotting capabilities.

Examples

To create a resolution object, apply the .resolution() method to a Series object with statistic=None.

co2ts = pyleo.utils.load_dataset('AACO2')
edc = pyleo.utils.load_dataset('EDC-dD')
ms = edc & co2ts # create MS object
ms_resolution = ms.resolution(statistic=None)

Time unit not found, attempting conversion.
Converted to ky BP

Several methods are then available:

Summary statistics can be obtained via .describe()

ms_resolution.describe()

{'EPICA Dome C dD': {'nobs': np.int64(5784),
  'minmax': (np.float64(0.008244210009855486), np.float64(1.3640000000207237)),
  'mean': np.float64(0.13859329637102363),
  'variance': np.float64(0.029806736482500852),
  'skewness': np.float64(2.6618614618357794),
  'kurtosis': np.float64(8.705801510816693),
  'median': np.float64(0.05813225000042621)},
 'EPICA Dome C CO2': {'nobs': np.int64(1900),
  'minmax': (np.float64(1.9999983537945243e-05),
   np.float64(4.171250000015277)),
  'mean': np.float64(0.4240631052631468),
  'variance': np.float64(0.24737134999156402),
  'skewness': np.float64(2.0625788668360423),
  'kurtosis': np.float64(6.7967879720624484),
  'median': np.float64(0.27533999998908243)}}

A simple plot can be created using .summary_plot()

ms_resolution.summary_plot()

(<Figure size 1000x800 with 1 Axes>, <Axes: xlabel='Resolution [ky BP]'>)

sel(value=None, time=None, tolerance=0)[source]

Slice MulitpleSeries based on ‘value’ or ‘time’. See examples in pyleoclim.series.Series for usage.

Parameters:

value (int, float, slice) – If int/float, then the Series will be sliced so that self.value is equal to value (+/- tolerance). If slice, then the Series will be sliced so self.value is between slice.start and slice.stop (+/- tolerance).
time (int, float, slice) – If int/float, then the Series will be sliced so that self.time is equal to time. (+/- tolerance) If slice of int/float, then the Series will be sliced so that self.time is between slice.start and slice.stop. If slice of datetime (or str containing datetime, such as ‘2020-01-01’), then the Series will be sliced so that self.datetime_index is between time.start and time.stop (+/- tolerance, which needs to be a timedelta).
tolerance (int, float, default 0.) – Used by value and time, see above.

Returns:

ms_new – Copy of self, sliced according to value and time.

Return type:

pyleoclim.mulitpleseries.MultipleSeries

See also

pyleoclim.series.Series.sel: Slicing a series by value and time.

spectral(method='lomb_scargle', freq=None, settings=None, mute_pbar=False, freq_kwargs=None, label=None, verbose=False, scalogram_list=None)[source]

Perform spectral analysis on the timeseries

Parameters:

method (str; {'wwz', 'mtm', 'lomb_scargle', 'welch', 'periodogram', 'cwt'}) –
freq (str or array, optional) – Information to produce the frequency vector. This can be ‘log’,’scale’, ‘nfft’, ‘lomb_scargle’ or ‘welch’ or a NumPy array. If a string, will use make_freq_vector with the specified frequency-generating method. If an array, this will be passed directly to the spectral method. If None (default), will use ‘log’ for WWZ and ‘lomb_scargle’ for Lomb-Scargle. This parameter is highly consequential for the WWZ and Lomb-Scargle methods, but is otherwise ignored, as other spectral methods generate their frequency vector internally.
freq_kwargs (dict) – Arguments for frequency vector
settings (dict) – Arguments for the specific spectral method
label (str) – Label for the PSD object
verbose (bool) – If True, will print warning messages if there is any
mute_pbar (bool) – Mute the progress bar. Default is False.
scalogram_list (pyleoclim.MultipleScalogram) – Multiple scalogram object containing pre-computed scalograms to use when calculating spectra, only works with wwz or cwt

Returns:

psd – A Multiple PSD object

Return type:

MultiplePSD

See also

pyleoclim.utils.spectral.mtm: Spectral analysis using the Multitaper approach
pyleoclim.utils.spectral.lomb_scargle: Spectral analysis using the Lomb-Scargle method
pyleoclim.utils.spectral.welch: Spectral analysis using the Welch segement approach
pyleoclim.utils.spectral.periodogram: Spectral anaysis using the basic Fourier transform
pyleoclim.utils.spectral.wwz_psd: Spectral analysis using the Wavelet Weighted Z transform
pyleoclim.utils.spectral.cwt_psd: Spectral analysis using the continuous Wavelet Transform as implemented by Torrence and Compo
pyleoclim.utils.wavelet.make_freq_vector: Functions to create the frequency vector
pyleoclim.utils.tsutils.detrend: Detrending function
pyleoclim.core.series.Series.spectral: Spectral analysis for a single timeseries
pyleoclim.core.PSD.PSD: PSD object
pyleoclim.core.psds.MultiplePSD: Multiple PSD object

Examples

air = pyleo.utils.load_dataset('AIR')
nino = pyleo.utils.load_dataset('NINO3')
ms = air & nino
ms_psd = ms.spectral(method='mtm')
ms_psd.plot()

(<Figure size 1000x400 with 1 Axes>,
 <Axes: xlabel='Period [year]', ylabel='PSD'>)

stackplot(figsize=None, savefig_settings=None, time_unit=None, xlim=None, fill_between_alpha=0.2, colors=None, cmap='tab10', norm=None, labels='auto', ylabel_fontsize=8, spine_lw=1.5, grid_lw=0.5, label_x_loc=-0.15, v_shift_factor=0.75, linewidth=1.5, yticks_minor=False, xticks_minor=False, ylims='auto', plot_kwargs=None)[source]

Stack plot of multiple series

Time units are harmonized prior to plotting. Note that the plotting style is uniquely designed for this one and cannot be properly reset with pyleoclim.set_style().

Parameters:

figsize (list) – Size of the figure.
savefig_settings (dictionary) –
the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existing or non-existing path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”} The default is None.
time_unit (str) –
the target time unit, possible inputs: {

’year’, ‘years’, ‘yr’, ‘yrs’, ‘y BP’, ‘yr BP’, ‘yrs BP’, ‘year BP’, ‘years BP’, ‘ky BP’, ‘kyr BP’, ‘kyrs BP’, ‘ka BP’, ‘ka’, ‘my BP’, ‘myr BP’, ‘myrs BP’, ‘ma BP’, ‘ma’,

} default is None, in which case the code picks the most common time unit in the collection. If no discernible winner can be found, the unit of the first series in the collection is used.
xlim (list) – The x-axis limit.
fill_between_alpha (float) – The transparency for the fill_between shades.
colors (a list of, or one, Python supported color code (a string of hex code or a tuple of rgba values)) – Colors for plotting. If None, the plotting will cycle the ‘tab10’ colormap; if only one color is specified, then all curves will be plotted with that single color; if a list of colors are specified, then the plotting will cycle that color list.
cmap (str) – The colormap to use when “colors” is None.
norm (matplotlib.colors.Normalize like) – The normalization for the colormap. If None, a linear normalization will be used.
labels (None, 'auto' or list) – If None, doesn’t add labels to the subplots If ‘auto’, uses the labels passed during the creation of pyleoclim.Series If list, pass a list of strings for each labels. Default is ‘auto’
spine_lw (float) – The linewidth for the spines of the axes.
grid_lw (float) – The linewidth for the gridlines.
label_x_loc (float) – The x location for the label of each curve.
v_shift_factor (float) – The factor for the vertical shift of each axis. The default value 3/4 means the top of the next axis will be located at 3/4 of the height of the previous one.
linewidth (float) – The linewidth for the curves.
ylabel_fontsize (int) – Size for ylabel font. Default is 8, to avoid crowding.
yticks_minor (bool) – Whether the y axes should contain minor ticks (use sparingly!). Default: False
xticks_minor (bool) – Whether the x axis should contain minor ticks. Default: False
ylims (str {'spacious', 'auto'}) – Method for determining the limits of the y axes. Default is ‘spacious’, which is mean +/- 4 x std ‘auto’ activates the Matplotlib default
plot_kwargs (dict or list of dict) –
Arguments to further customize the plot from matplotlib.pyplot.plot.
- Dictionary: Arguments will be applied to all lines in the stackplots
- List of dictionaries: Allows to customize one line at a time.

Returns:

fig (matplotlib.figure) – the figure object from matplotlib See [matplotlib.pyplot.figure](https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.figure.html) for details.
ax (matplotlib.axis) – the axis object from matplotlib See [matplotlib.axes](https://matplotlib.org/api/axes_api.html) for details.

See also

pyleoclim.utils.plotting.savefig: Saving figure in Pyleoclim

Examples

soi = pyleo.utils.load_dataset('SOI')
nino = pyleo.utils.load_dataset('NINO3')
ms = soi & nino
fig, ax = ms.stackplot()

Let’s change the labels on the left

fig, ax = ms.stackplot(labels=['SOI','NINO3'])

And let’s remove them completely

fig, ax = ms.stackplot(labels=None)

Now, let’s add markers to the timeseries.

fig, ax = ms.stackplot(labels=None, plot_kwargs={'marker':'o'})

Using different marker types on each series:

fig, ax = ms.stackplot(labels=None, plot_kwargs=[{'marker':'o'},{'marker':'^'}])

By default, the y axes are kept very minimal to allow stacking many records. In some instances, however, one may want more detailed axes, with major and minor ticks. We also show how to enlarge the ylabels and adjust vertical spacing for improved readability:

fig, ax = ms.stackplot(labels=None, ylabel_fontsize = 12, v_shift_factor = 0.9,
                       yticks_minor=True, xticks_minor=True, ylims='auto')

This approach makes sense with small stacks, but quickly becomes unwieldy with large ones. Use at your own risk!

standardize()[source]

Standardize each series object in a collection

Returns:: ms – The standardized pyleoclim.MultipleSeries object
Return type:: MultipleSeries

Examples

soi = pyleo.utils.load_dataset('SOI')
nino = pyleo.utils.load_dataset('NINO3')
ms = soi & nino
ms_std = ms.standardize()

stripes(cmap='RdBu_r', sat=1.0, ref_period=None, figsize=None, savefig_settings=None, time_unit=None, labels='auto', label_color='gray', show_xaxis=False, common_time_kwargs=None, xlim=None, font_scale=0.8, x_offset=0.05)[source]

Represents a MultipleSeries object as a quilt of Ed Hawkins’ “stripes” patterns

To ensure comparability, constituent series are placed on a common time axis, using MultipleSeries.common_time(). To ensure consistent scaling, all series are Gaussianized prior to plotting.

Credit: https://showyourstripes.info/, Implementation: https://matplotlib.org/matplotblog/posts/warming-stripes/

Parameters:

cmap (str) – colormap name (https://matplotlib.org/stable/tutorials/colors/colormaps.html) Default is ‘RdBu_r’
ref_period (TYPE, optional) – dates of the reference period, in the form “(first, last)”. The default is None, which will pick the beginning and end of the common time axis.
figsize (list) – a list of two integers indicating the figure size (in inches)
sat (float > 0) – Controls the saturation of the colormap normalization by scaling the vmin, vmax in https://matplotlib.org/stable/tutorials/colors/colormapnorms.html default = 1.0
show_xaxis (bool) – flag indicating whether or not the x-axis should be shown (default = False)
savefig_settings (dictionary) –
the dictionary of arguments for plt.savefig(); some notes below:
- ’path’ must be specified; it can be any existing or non-existing path, with or without a suffix; if the suffix is not given in ‘path’, it will follow ‘format’
- ’format’ can be one of {“pdf”, ‘eps’, ‘png’, ps’} The default is None.
time_unit (str) –
the target time unit, possible inputs: {

’year’, ‘years’, ‘yr’, ‘yrs’, ‘y BP’, ‘yr BP’, ‘yrs BP’, ‘year BP’, ‘years BP’, ‘ky BP’, ‘kyr BP’, ‘kyrs BP’, ‘ka BP’, ‘ka’, ‘my BP’, ‘myr BP’, ‘myrs BP’, ‘ma BP’, ‘ma’,

} default is None, in which case the code picks the most common time unit in the collection. If no discernible winner can be found, the unit of the first series in the collection is used.
xlim (list) – The x-axis limit.
x_offset (float) – value controlling the horizontal offset between stripes and labels (default = 0.05)
labels (None, 'auto' or list) –
If None, doesn’t add labels to the subplots

If ‘auto’, uses the labels passed during the creation of pyleoclim.Series

If list, pass a list of strings for each labels. Default is ‘auto’
common_time_kwargs (dict) – Optional arguments for common_time()
font_scale (float) – The scale for the font sizes. Default is 0.8.

Returns:

fig (matplotlib.figure) – the figure object from matplotlib See [matplotlib.pyplot.figure](https://matplotlib.org/stable/api/figure_api.html) for details.
ax (matplotlib.axis) – the axis object from matplotlib See [matplotlib.axes](https://matplotlib.org/stable/api/axes_api.html) for details.

See also

pyleoclim.core.multipleseries.MultipleSeries.common_time: aligns the time axes of a MultipleSeries object
pyleoclim.utils.plotting.savefig: saving a figure in Pyleoclim
pyleoclim.core.series.Series.stripes: stripes representation in Pyleoclim
pyleoclim.utils.tsutils.gaussianize: mapping to a standard Normal distribution

Examples

co2ts = pyleo.utils.load_dataset('AACO2')
lr04 = pyleo.utils.load_dataset('LR04')
edc = pyleo.utils.load_dataset('EDC-dD')
ms = lr04.flip() & edc & co2ts # create MS object
fig, ax = ms.stripes()

The two series have different lengths, left: 2115 vs right: 1901
Metadata are different:
value_unit property -- left: ‰, right: ppm
value_name property -- left: $\delta^{18} \mathrm{O}$ x (-1), right: $CO_2$
label property -- left: LR04 benthic stack, right: EPICA Dome C CO2
archiveType property -- left: MarineSediment, right: GlacierIce
importedFrom property -- left: None, right: https://www.ncei.noaa.gov/pub/data/paleo/icecore/antarctica/antarctica2015co2composite.txt
The two series have different lengths, left: 5785 vs right: 1901
Metadata are different:
lat property -- left: -75.1011, right: None
lon property -- left: 123.3478, right: None
elevation property -- left: 3233, right: None
time_unit property -- left: y BP, right: ky BP
value_unit property -- left: ‰, right: ppm
value_name property -- left: $\delta \mathrm{D}$, right: $CO_2$
label property -- left: EPICA Dome C dD, right: EPICA Dome C CO2
sensorType property -- left: ice sheet, right: None
observationType property -- left: hydrogen isotopes, right: None
importedFrom property -- left: https://www.ncei.noaa.gov/pub/data/paleo/icecore/antarctica/epica_domec/edc3deuttemp2007.txt, right: https://www.ncei.noaa.gov/pub/data/paleo/icecore/antarctica/antarctica2015co2composite.txt
control_archiveType property -- left: False, right: None

The default style has rather thick bands, intense colors, and too many stripes. The first issue can be solved by passing a figsize tuple; the second by increasing the LIM parameter; the third by passing a step of 0.5 (500y) to common_time(). Finally, the labels are too close to the edge of the plot, which can be adjusted with x_offset, like so:

co2ts = pyleo.utils.load_dataset('AACO2')
lr04 = pyleo.utils.load_dataset('LR04')
edc = pyleo.utils.load_dataset('EDC-dD')
ms = lr04.flip() & edc & co2ts # create MS object
fig, ax = ms.stripes(figsize=(8,2.5),show_xaxis=True, sat = 0.8)

The two series have different lengths, left: 2115 vs right: 1901
Metadata are different:
value_unit property -- left: ‰, right: ppm
value_name property -- left: $\delta^{18} \mathrm{O}$ x (-1), right: $CO_2$
label property -- left: LR04 benthic stack, right: EPICA Dome C CO2
archiveType property -- left: MarineSediment, right: GlacierIce
importedFrom property -- left: None, right: https://www.ncei.noaa.gov/pub/data/paleo/icecore/antarctica/antarctica2015co2composite.txt
The two series have different lengths, left: 5785 vs right: 1901
Metadata are different:
lat property -- left: -75.1011, right: None
lon property -- left: 123.3478, right: None
elevation property -- left: 3233, right: None
time_unit property -- left: y BP, right: ky BP
value_unit property -- left: ‰, right: ppm
value_name property -- left: $\delta \mathrm{D}$, right: $CO_2$
label property -- left: EPICA Dome C dD, right: EPICA Dome C CO2
sensorType property -- left: ice sheet, right: None
observationType property -- left: hydrogen isotopes, right: None
importedFrom property -- left: https://www.ncei.noaa.gov/pub/data/paleo/icecore/antarctica/epica_domec/edc3deuttemp2007.txt, right: https://www.ncei.noaa.gov/pub/data/paleo/icecore/antarctica/antarctica2015co2composite.txt
control_archiveType property -- left: False, right: None

time_coverage_plot(figsize=[10, 3], marker=None, markersize=None, alpha=0.8, linestyle=None, linewidth=10, colors=None, cmap='turbo', norm=None, xlabel=None, ylabel=None, title=None, time_unit=None, legend=True, inline_legend=True, plot_kwargs=None, lgd_kwargs=None, label_x_offset=200, label_y_offset=0, savefig_settings=None, ax=None, ypad=None, invert_xaxis=False, invert_yaxis=False)[source]

A plot of the temporal coverage of the records in a MultipleSeries object organized by ranked length.

Inspired by Dr. Mara Y. McPartland.

Parameters:

figsize (list, optional) – Size of the figure. The default is [10, 4].
marker (str, optional) – Marker type. The default is None.
markersize (float, optional) – Marker size. The default is None.
alpha (float, optional) – Alpha of the lines
linestyle (str, optional) – Line style. The default is None.
linewidth (float, optional) – The width of the line. The default is 10.
colors (a list of, or one, Python supported color code (a string of hex code or a tuple of rgba values)) – Colors for plotting. If None, the plotting will cycle the ‘viridis’ colormap; if only one color is specified, then all curves will be plotted with that single color; if a list of colors are specified, then the plotting will cycle that color list.
cmap (str) – The colormap to use when “colors” is None. Default is ‘turbo’
norm (matplotlib.colors.Normalize) – The normalization for the colormap. If None, a linear normalization will be used.
xlabel (str, optional) – x-axis label. The default is None.
ylabel (str, optional) – y-axis label. The default is None.
title (str, optional) – Title. The default is None.
time_unit (str) –
the target time unit, possible input: {

’year’, ‘years’, ‘yr’, ‘yrs’, ‘y BP’, ‘yr BP’, ‘yrs BP’, ‘year BP’, ‘years BP’, ‘ky BP’, ‘kyr BP’, ‘kyrs BP’, ‘ka BP’, ‘ka’, ‘my BP’, ‘myr BP’, ‘myrs BP’, ‘ma BP’, ‘ma’,

} default is None, in which case the code picks the most common time unit in the collection. If no unambiguous winner can be found, the unit of the first series in the collection is used.
legend (bool, optional) – Whether the show the legend. The default is True.
inline_legend (bool, optional) – Whether to use inline labels or the default pyleoclim legend. This option overrides lgd_kwargs
plot_kwargs (dict, optional) – Plot parameters. The default is None.
lgd_kwargs (dict, optional) –
Legend parameters. The default is None.

If inline_legend is True, lgd_kwargs will be passed to ax.text() (see matplotlib.axes.Axes.text documentation) If inline_legend is False, lgd_kwargs will be passed to ax.legend() (see matplotlib.axes.Axes.legend documentation)
label_x_offset (float or list, optional) – Amount to offset label by in the x direction. Only used if inline_legend is True. Default is 200. If list, should have the same number of elements as the MultipleSeries object.
label_y_offset (float or list, optional) – Amount to offset label by in the y direction. Only used if inline_legend is True. Default is 0. If list, should have the same number of elements as the MultipleSeries object.
savefig_settings (dictionary, optional) –
the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existing or non-existing path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”} The default is None.
ax (matplotlib.ax, optional) – The matplotlib axis onto which to return the figure. The default is None.
invert_xaxis (bool, optional) – if True, the x-axis of the plot will be inverted
invert_yaxis (bool, optional) – if True, the y-axis of the plot will be inverted

Returns:

fig (matplotlib.figure) – the figure object from matplotlib See [matplotlib.pyplot.figure](https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.figure.html) for details.
ax (matplotlib.axis) – the axis object from matplotlib See [matplotlib.axes](https://matplotlib.org/api/axes_api.html) for details.

See also

pyleoclim.utils.plotting.savefig: Saving figure in Pyleoclim

Examples

co2ts = pyleo.utils.load_dataset('AACO2')
lr04 = pyleo.utils.load_dataset('LR04')
edc = pyleo.utils.load_dataset('EDC-dD')
ms = lr04.flip() & edc & co2ts # create MS object
fig, ax = ms.time_coverage_plot(label_y_offset=-.08) #Fiddling with label offsets is sometimes necessary for aesthetic

The two series have different lengths, left: 2115 vs right: 1901
Metadata are different:
value_unit property -- left: ‰, right: ppm
value_name property -- left: $\delta^{18} \mathrm{O}$ x (-1), right: $CO_2$
label property -- left: LR04 benthic stack, right: EPICA Dome C CO2
archiveType property -- left: MarineSediment, right: GlacierIce
importedFrom property -- left: None, right: https://www.ncei.noaa.gov/pub/data/paleo/icecore/antarctica/antarctica2015co2composite.txt
The two series have different lengths, left: 5785 vs right: 1901
Metadata are different:
lat property -- left: -75.1011, right: None
lon property -- left: 123.3478, right: None
elevation property -- left: 3233, right: None
time_unit property -- left: y BP, right: ky BP
value_unit property -- left: ‰, right: ppm
value_name property -- left: $\delta \mathrm{D}$, right: $CO_2$
label property -- left: EPICA Dome C dD, right: EPICA Dome C CO2
sensorType property -- left: ice sheet, right: None
observationType property -- left: hydrogen isotopes, right: None
importedFrom property -- left: https://www.ncei.noaa.gov/pub/data/paleo/icecore/antarctica/epica_domec/edc3deuttemp2007.txt, right: https://www.ncei.noaa.gov/pub/data/paleo/icecore/antarctica/antarctica2015co2composite.txt
control_archiveType property -- left: False, right: None

Awkward vertical spacing can be adjusted by varying linewidth and figure size

co2ts = pyleo.utils.load_dataset('AACO2')
lr04 = pyleo.utils.load_dataset('LR04')
edc = pyleo.utils.load_dataset('EDC-dD')
ms = lr04.flip() & edc & co2ts # create MS object
fig, ax = ms.time_coverage_plot(linewidth=20,figsize=[10,2],label_y_offset=-.1)

The two series have different lengths, left: 2115 vs right: 1901
Metadata are different:
value_unit property -- left: ‰, right: ppm
value_name property -- left: $\delta^{18} \mathrm{O}$ x (-1), right: $CO_2$
label property -- left: LR04 benthic stack, right: EPICA Dome C CO2
archiveType property -- left: MarineSediment, right: GlacierIce
importedFrom property -- left: None, right: https://www.ncei.noaa.gov/pub/data/paleo/icecore/antarctica/antarctica2015co2composite.txt
The two series have different lengths, left: 5785 vs right: 1901
Metadata are different:
lat property -- left: -75.1011, right: None
lon property -- left: 123.3478, right: None
elevation property -- left: 3233, right: None
time_unit property -- left: y BP, right: ky BP
value_unit property -- left: ‰, right: ppm
value_name property -- left: $\delta \mathrm{D}$, right: $CO_2$
label property -- left: EPICA Dome C dD, right: EPICA Dome C CO2
sensorType property -- left: ice sheet, right: None
observationType property -- left: hydrogen isotopes, right: None
importedFrom property -- left: https://www.ncei.noaa.gov/pub/data/paleo/icecore/antarctica/epica_domec/edc3deuttemp2007.txt, right: https://www.ncei.noaa.gov/pub/data/paleo/icecore/antarctica/antarctica2015co2composite.txt
control_archiveType property -- left: False, right: None

to_csv(path=None, *args, use_common_time=False, **kwargs)[source]

Export MultipleSeries to CSV

Parameters:

path (str, optional) – system path to save the file. The default is None, in which case the filename defaults to the poetic ‘MultipleSeries.csv’ in the current directory.
*args – Arguments and keyword arguments to pass to common_time.
**kwargs – Arguments and keyword arguments to pass to common_time.
use_common_time – Set to True if you want to use common_time to align the Series to a common timescale. Else, times for which some Series don’t have values will be filled with NaN (default).
bool – Set to True if you want to use common_time to align the Series to a common timescale. Else, times for which some Series don’t have values will be filled with NaN (default).

Return type:

None.

Examples

This will place the NINO3 and SOI datasets into a MultipleSeries object and export it to enso.csv.

soi = pyleo.utils.load_dataset('SOI')
nino = pyleo.utils.load_dataset('NINO3')
ms = soi & nino
ms.label = 'enso'
ms.to_csv()

to_json(path=None)[source]

Export the pyleoclim.MultipleSeries object to a json file

Parameters:: path (string, optional) – The path to the file. The default is None, resulting in a file saved in the current working directory using the label for the dataset as filename if available or ‘mulitpleseries.json’ if label is not provided.
Return type:: None.

to_pandas(paleo_style=False, *args, use_common_time=False, **kwargs)[source]

Align Series and place in DataFrame.

Column names will be taken from each Series’ label.

Parameters:

paleo_style (boolean, optional) – If True, will format datetime as the common time vector and assign as index name the time_name of the first series in the object.
*args – Arguments and keyword arguments to pass to common_time.
**kwargs – Arguments and keyword arguments to pass to common_time.
use_common_time – Pass True if you want to use common_time to align the Series to have common times. Else, times for which some Series doesn’t have values will be filled with NaN (default).
bool – Pass True if you want to use common_time to align the Series to have common times. Else, times for which some Series doesn’t have values will be filled with NaN (default).

Return type:

pandas.DataFrame

view()[source]

Generates a DataFrame version of the MultipleSeries object, suitable for viewing in a Jupyter Notebook

Return type:: pd.DataFrame

Examples

soi = pyleo.utils.load_dataset('SOI')
nino = pyleo.utils.load_dataset('NINO3')
ms = soi & nino
ms.name = 'ENSO'
ms.view()

	Southern Oscillation Index	NINO3 SST
Time
1866.002894	-0.62	NaN
1866.087769	-0.12	NaN
1866.164431	-0.62	NaN
1866.249306	-0.65	NaN
1866.331443	0.04	NaN
...	...	...
2024.752349	0.35	NaN
2024.837224	0.61	NaN
2024.919362	1.09	NaN
2025.004237	0.26	NaN
2025.089112	0.21	NaN

3506 rows × 2 columns

wavelet(method='cwt', settings={}, freq=None, freq_kwargs=None, verbose=False, mute_pbar=False)[source]

Wavelet analysis

Parameters:

method (str {wwz, cwt}) –
- cwt - the continuous wavelet transform (as per Torrence and Compo [1998])
  is appropriate only for evenly-spaced series.
- wwz - the weighted wavelet Z-transform (as per Foster [1996])
  is appropriate for both evenly and unevenly-spaced series.
Default is cwt, returning an error if the Series is unevenly-spaced.
settings (dict) – Settings for the particular method. The default is {}.
freq (str or array, optional) – Information to produce the frequency vector (highly consequential for the WWZ method) This can be ‘log’,’scale’, ‘nfft’, ‘lomb_scargle’, ‘welch’ or a NumPy array. If a string, will use make_freq_vector() with the specified frequency-generating method. If an array, this will be passed directly to the spectral method. If None (default), will use the ‘log’ method
freq_kwargs (dict) – Arguments for frequency vector
settings – Arguments for the specific spectral method
verbose (bool) – If True, will print warning messages if there is any
mute_pbar (bool, optional) – Whether to mute the progress bar. The default is False.

Returns:

scals – A Multiple Scalogram object

Return type:

MultipleScalograms

See also

pyleoclim.utils.wavelet.wwz: wwz function
pyleoclim.utils.wavelet.cwt: cwt function
pyleoclim.utils.wavelet.make_freq_vector: Functions to create the frequency vector
pyleoclim.utils.tsutils.detrend: Detrending function
pyleoclim.core.series.Series.wavelet: wavelet analysis on single object
pyleoclim.core.scalograms.MultipleScalogram: Multiple Scalogram object

References

Torrence, C. and G. P. Compo, 1998: A Practical Guide to Wavelet Analysis. Bull. Amer. Meteor. Soc., 79, 61-78. Python routines available at http://paos.colorado.edu/research/wavelets/

Examples

soi = pyleo.utils.load_dataset('SOI')
nino = pyleo.utils.load_dataset('NINO3')
ms = (soi & nino)
wav = ms.wavelet(method='wwz')

MultipleGeoSeries (pyleoclim.MultipleGeoSeries)

class pyleoclim.core.multiplegeoseries.MultipleGeoSeries(series_list, time_unit=None, label=None)[source]

MultipleGeoSeries object.

This object handles a collection of the type GeoSeries and can be created from a list of such objects. MultipleGeoSeries should be used when the need to run analysis on multiple records arises, such as running principal component analysis. Some of the methods automatically transform the time axis prior to analysis to ensure consistency.

Parameters:

series_list (list) – a list of pyleoclim.Series objects
time_unit (str) – The target time unit for every series in the list. If None, then no conversion will be applied; Otherwise, the time unit of every series in the list will be converted to the target.
label (str) – label of the collection of timeseries (e.g. ‘Euro 2k’)

Examples

from pylipd.utils.dataset import load_dir
lipd = load_dir(name='Pages2k')
df = lipd.get_timeseries_essentials()
dfs = df.query("archiveType in ('Wood','Documents','Coral','Lake sediment')")
# place in a MultipleGeoSeries object
ts_list = []
for _, row in dfs.iterrows():
    ts_list.append(pyleo.GeoSeries(time=row['time_values'],value=row['paleoData_values'],
                                   time_name=row['time_variableName'],value_name=row['paleoData_variableName'],
                                   time_unit=row['time_units'], value_unit=row['paleoData_units'],
                                   lat = row['geo_meanLat'], lon = row['geo_meanLon'],
                                   archiveType = row['archiveType'], verbose = False,
                                   label=row['dataSetName']+'_'+row['paleoData_variableName']))

Euro2k = pyleo.MultipleGeoSeries(ts_list, label='Euro2k',time_unit='years AD')
Euro2k.map()

Loading 16 LiPD files

Loaded..

(<Figure size 1600x600 with 2 Axes>,
 {'map': <GeoAxes: xlabel='lon', ylabel='lat'>, 'leg': <Axes: >})

map(marker='archiveType', hue='archiveType', size=None, cmap=None, edgecolor='k', projection='auto', proj_default=True, crit_dist=5000, colorbar=True, color_scale_type=None, background=True, borders=False, coastline=True, rivers=False, lakes=False, land=True, ocean=True, figsize=None, fig=None, scatter_kwargs=None, gridspec_kwargs=None, legend=True, gridspec_slot=None, lgd_kwargs=None, savefig_settings=None, **kwargs)[source]

Mapping of the collection of GeoSeries objects.

Parameters:

huestring, optional

Grouping variable that will produce points with different colors. Can be either categorical or numeric, although color mapping will behave differently in latter case. The default is ‘archiveType’.

sizestring, optional: Grouping variable that will produce points with different sizes. Expects to be numeric. Any data without a value for the size variable will be filtered out. The default is None.
markerstring, optional: Grouping variable that will produce points with different markers. Can have a numeric dtype but will always be treated as categorical. The default is ‘archiveType’.
edgecolorcolor (string) or list of rgba tuples, optional: Color of marker edge. The default is ‘w’.
projectionstring: the map projection. Available projections: ‘Robinson’ (default), ‘PlateCarree’, ‘AlbertsEqualArea’, ‘AzimuthalEquidistant’,’EquidistantConic’,’LambertConformal’, ‘LambertCylindrical’,’Mercator’,’Miller’,’Mollweide’,’Orthographic’, ‘Sinusoidal’,’Stereographic’,’TransverseMercator’,’UTM’, ‘InterruptedGoodeHomolosine’,’RotatedPole’,’OSGB’,’EuroPP’, ‘Geostationary’,’NearsidePerspective’,’EckertI’,’EckertII’, ‘EckertIII’,’EckertIV’,’EckertV’,’EckertVI’,’EqualEarth’,’Gnomonic’, ‘LambertAzimuthalEqualArea’,’NorthPolarStereo’,’OSNI’,’SouthPolarStereo’ By default, projection == ‘auto’, so the projection will be picked based on the degree of clustering of the sites.
proj_defaultbool, optional: If True, uses the standard projection attributes. Enter new attributes in a dictionary to change them. Lists of attributes can be found in the Cartopy documentation. The default is True.
crit_distfloat, optional: critical radius for projection choice. Default: 5000 km Only active if projection == ‘auto’
backgroundbool, optional: If True, uses a shaded relief background (only one available in Cartopy) Default is on (True).
bordersbool or dict, optional: Draws the countries border. If a dictionary of formatting arguments is supplied (e.g. linewidth, alpha), will draw according to specifications. Defaults is off (False).
coastlinebool or dict, optional: Draws the coastline. If a dictionary of formatting arguments is supplied (e.g. linewidth, alpha), will draw according to specifications. Defaults is on (True).
landbool or dict, optional: Colors land masses. If a dictionary of formatting arguments is supplied (e.g. color, alpha), will draw according to specifications. Default is off (True). Overriden if background=True.
oceanbool or dict, optional: Colors oceans. If a dictionary of formatting arguments is supplied (e.g. color, alpha), will draw according to specifications. Default is on (True). Overriden if background=True.
riversbool or dict, optional: Draws major rivers. If a dictionary of formatting arguments is supplied (e.g. linewidth, alpha), will draw according to specifications. Default is off (False).
lakesbool or dict, optional: Draws major lakes. If a dictionary of formatting arguments is supplied (e.g. color, alpha), will draw according to specifications. Default is off (False).
figsizelist or tuple, optional: Size for the figure
scatter_kwargsdict, optional: Dict of arguments available in seaborn.scatterplot. Dictionary of arguments available in matplotlib.pyplot.scatter.
legendbool, optional: Whether to draw a legend on the figure. Default is True.
colorbarbool, optional: Whether to draw a colorbar on the figure if the data associated with hue are numeric. Default is True.
color_scale_typestr, optional: Setting to “discrete” will force a discrete color scale with a default bin number of max(11, n) where n=number of unique values$^{

rac{1}{2}}$

Default is None

lgd_kwargsdict, optional

Dictionary of arguments for matplotlib.pyplot.legend.

savefig_settingsdict, optional

Dictionary of arguments for matplotlib.pyplot.saveFig.

“path” must be specified; it can be any existed or non-existed path, with or without a suffix; if the suffix is not given in “path”, it will follow “format”

“format” can be one of {“pdf”, “eps”, “png”, “ps”}

extentTYPE, optional

DESCRIPTION. The default is ‘global’.

cmapstring or list, optional

Matplotlib supported colormap id or list of colors for creating a colormap. See choosing a matplotlib colormap. The default is None.

figmatplotlib.pyplot.figure, optional

See matplotlib.pyplot.figure <https://matplotlib.org/3.5.0/api/_as_gen/matplotlib.pyplot.figure.html#matplotlib-pyplot-figure>_. The default is None.

gs_slotGridspec slot, optional

If generating a map for a multi-plot, pass a gridspec slot. The default is None.

gridspec_kwargsdict, optional

Function assumes the possibility of a colorbar, map, and legend. A list of floats associated with the keyword width_ratios will assume the first (index=0) is the relative width of the colorbar, the second to last (index = -2) is the relative width of the map, and the last (index = -1) is the relative width of the area for the legend. For information about Gridspec configuration, refer to `Matplotlib documentation <https://matplotlib.org/3.5.0/api/_as_gen/matplotlib.gridspec.GridSpec.html#matplotlib.gridspec.GridSpec>_. The default is None.

kwargs: dict, optional

‘missing_val_hue’, ‘missing_val_marker’, ‘missing_val_label’ can all be used to change the way missing values are represented (‘k’, ‘?’, are default hue and marker values will be associated with the label: ‘missing’).
‘hue_mapping’ and ‘marker_mapping’ can be used to submit dictionaries mapping hue values to colors and marker values to markers. Does not replace passing a string value for hue or marker.

Returns:

fig, ax_d: Matplotlib figure, dictionary of ax objects which includes the as many as three items: ‘cb’ (colorbar ax), ‘map’ (scatter map), and ‘leg’ (legend ax)

See also

pyleoclim.utils.mapping.scatter_map: information-rich scatterplot on Cartopy map

Examples

from pylipd.utils.dataset import load_dir
lipd = load_dir(name='Pages2k')
df = lipd.get_timeseries_essentials()
dfs = df.query("archiveType in ('Wood','Documents','Coral','Lake sediment','Borehole')")

# place in a MultipleGeoSeries object
ts_list = []
for _, row in dfs.iterrows():
    ts_list.append(pyleo.GeoSeries(time=row['time_values'],value=row['paleoData_values'],
                                   time_name=row['time_variableName'],value_name=row['paleoData_variableName'],
                                   time_unit=row['time_units'], value_unit=row['paleoData_units'],
                                   lat = row['geo_meanLat'], lon = row['geo_meanLon'],
                                   elevation = row['geo_meanElev'], observationType = row['paleoData_proxy'],
                                   archiveType = row['archiveType'], verbose = False,
                                   label=row['dataSetName']+'_'+row['paleoData_variableName']))

Euro2k = pyleo.MultipleGeoSeries(ts_list, label='Euro2k',time_unit='years AD')
Euro2k.map()

Loading 16 LiPD files

Loaded..

(<Figure size 1800x600 with 2 Axes>,
 {'map': <GeoAxes: xlabel='lon', ylabel='lat'>, 'leg': <Axes: >})

By default, a projection is picked based on the degree of geographic clustering of the sites. To focus on Europe and use a more local projection, do:

eur_coord = {'central_latitude':45, 'central_longitude':20}
Euro2k.map(projection='Orthographic',proj_default=eur_coord)

(<Figure size 1800x700 with 2 Axes>,
 {'map': <GeoAxes: xlabel='lon', ylabel='lat'>, 'leg': <Axes: >})

By default, the shape and colors of symbols denote proxy archives; however, one can use either graphical device to convey other information. For instance, if elevation is available, it may be displayed by size, like so:

Euro2k.map(projection='Orthographic', size='elevation', proj_default=eur_coord)

(<Figure size 1800x700 with 2 Axes>,
 {'map': <GeoAxes: xlabel='lon', ylabel='lat'>, 'leg': <Axes: >})

Same with observationType:

Euro2k.map(projection='Orthographic', hue = 'observationType', proj_default=eur_coord)

(<Figure size 1800x700 with 2 Axes>,
 {'map': <GeoAxes: xlabel='lon', ylabel='lat'>, 'leg': <Axes: >})

All three sources of information may be combined, but the figure height will need to be enlarged manually to fit the legend:

Euro2k.map(projection='Orthographic',hue='observationType',
           size='elevation', proj_default=eur_coord, figsize=[18, 8])

(<Figure size 1800x800 with 2 Axes>,
 {'map': <GeoAxes: xlabel='lon', ylabel='lat'>, 'leg': <Axes: >})

pca(weights=None, missing='fill-em', tol_em=0.005, max_em_iter=100, **pca_kwargs)[source]

Principal Component Analysis (Empirical Orthogonal Functions)

Decomposition of MultipleGeoSeries object in terms of orthogonal basis functions. Tolerant to missing values, infilled by an EM algorithm.

Do make sure the time axes are aligned, however! (e.g. use common_time())

Algorithm from statsmodels: https://www.statsmodels.org/stable/generated/statsmodels.multivariate.pca.PCA.html

Parameters:

weights (ndarray, optional) – Series weights to use after transforming data according to standardize or demean when computing the principal components.
missing ({str, None}) –
Method for missing data. Choices are:
- ’drop-row’ - drop rows with missing values.
- ’drop-col’ - drop columns with missing values.
- ’drop-min’ - drop either rows or columns, choosing by data retention.
- ’fill-em’ - use EM algorithm to fill missing value [ default]. ncomp should be set to the number of factors required.
- None raises if data contains NaN values.
tol_em (float) – Tolerance to use when checking for convergence of the EM algorithm.
max_em_iter (int) – Maximum iterations for the EM algorithm.

Returns:

res – Resulting pyleoclim.MultivariateDecomp object

Return type:

MultivariateDecomp

See also

pyleoclim.utils.tsutils.eff_sample_size: Effective Sample Size of timeseries y
pyleoclim.core.multivardecomp.MultivariateDecomp: The multivariate decomposition object
pyleoclim.core.mulitpleseries.MulitpleSeries.common_time: align time axes

Examples

from pylipd.utils.dataset import load_dir
lipd = load_dir(name='Pages2k') # this loads a small subset of the PAGES 2k database
lipd_euro = lipd.filter_by_geo_bbox(-20,20,40,80)
df = lipd_euro.get_timeseries_essentials()
dfs = df.query("archiveType in ('Wood') & paleoData_variableName not in ('year')")

# place in a MultipleGeoSeries object
ts_list = []
for _, row in dfs.iterrows():
    ts_list.append(pyleo.GeoSeries(time=row['time_values'],value=row['paleoData_values'],
                                   time_name=row['time_variableName'],value_name=row['paleoData_variableName'],
                                   time_unit=row['time_units'], value_unit=row['paleoData_units'],
                                   lat = row['geo_meanLat'], lon = row['geo_meanLon'],
                                   elevation = row['geo_meanElev'], observationType = row['paleoData_proxy'],
                                   archiveType = row['archiveType'], verbose = False,
                                   label=row['dataSetName']+'_'+row['paleoData_variableName']))

Euro2k = pyleo.MultipleGeoSeries(ts_list, label='Euro2k',time_unit='years AD')

res = Euro2k.common_time().pca() # carry out PCA
type(res) # the result is a MultivariateDecomp object

Loading 16 LiPD files

Loaded..

pyleoclim.core.multivardecomp.MultivariateDecomp

To plot the eigenvalue spectrum:

res.screeplot()

The provided eigenvalue array has only one dimension. UQ defaults to NB82

(<Figure size 600x400 with 1 Axes>,
 <Axes: title={'center': 'Euro2k PCA eigenvalues'}, xlabel='Mode index $i$', ylabel='$\\lambda_i$'>)

To plot the first mode, equivalent to res.modeplot(index=0):

res.modeplot()

(<Figure size 800x800 with 4 Axes>,
 {'pc': <Axes: xlabel='Time [years AD]', ylabel='$PC_1$'>,
  'psd': <Axes: xlabel='Period [years]', ylabel='PSD'>,
  'map': {'cb': <Axes: ylabel='EOF'>,
   'map': <GeoAxes: xlabel='lon', ylabel='lat'>}})

To plot the second (note the zero-based indexing):

res.modeplot(index=1)

(<Figure size 800x800 with 4 Axes>,
 {'pc': <Axes: xlabel='Time [years AD]', ylabel='$PC_2$'>,
  'psd': <Axes: xlabel='Period [years]', ylabel='PSD'>,
  'map': {'cb': <Axes: ylabel='EOF'>,
   'map': <GeoAxes: xlabel='lon', ylabel='lat'>}})

One can use map semantics to display the observation type as well:

res.modeplot(index=1, marker='observationType', size='elevation')

(<Figure size 800x800 with 5 Axes>,
 {'pc': <Axes: xlabel='Time [years AD]', ylabel='$PC_2$'>,
  'psd': <Axes: xlabel='Period [years]', ylabel='PSD'>,
  'map': {'cb': <Axes: ylabel='EOF'>,
   'map': <GeoAxes: xlabel='lon', ylabel='lat'>,
   'leg': <Axes: >}})

There are many ways to configure the map component. As a simple example, specifying the projection:

res.modeplot(index=1, marker='observationType', size='elevation',
    map_kwargs={'projection':'Robinson'})

(<Figure size 800x800 with 5 Axes>,
 {'pc': <Axes: xlabel='Time [years AD]', ylabel='$PC_2$'>,
  'psd': <Axes: xlabel='Period [years]', ylabel='PSD'>,
  'map': {'cb': <Axes: ylabel='EOF'>,
   'map': <GeoAxes: xlabel='lon', ylabel='lat'>,
   'leg': <Axes: >}})

Or dive into the nuances of gridspec and legend configurations:

res.modeplot(index=1, marker='observationType', size='elevation',
            map_kwargs={'projection':'Robinson',
                        'gridspec_kwargs': {'width_ratios': [.5, 1,14, 4], 'wspace':-.065},
                        'lgd_kwargs':{'bbox_to_anchor':[-.015,1]}})

(<Figure size 800x800 with 5 Axes>,
 {'pc': <Axes: xlabel='Time [years AD]', ylabel='$PC_2$'>,
  'psd': <Axes: xlabel='Period [years]', ylabel='PSD'>,
  'map': {'cb': <Axes: ylabel='EOF'>,
   'map': <GeoAxes: xlabel='lon', ylabel='lat'>,
   'leg': <Axes: >}})

One can also configure how hue information gets displayed:

# Dashboard with hue as a legend category
res.modeplot(index=1, size='elevation', map_kwargs= dict(colorbar=False))

# Dashboard with discrete colorbar
res.modeplot(index=1, size='elevation', map_kwargs= dict(color_scale_type='discrete'))

# Dashboard with custom scalar mappable
sm = pyleo.utils.mapping.make_scalar_mappable(cmap='vlag', hue_vect=res.eigvecs[:, 1], n=21,norm_kwargs={'vcenter': -.5})
res.modeplot(index=1, size='elevation', map_kwargs= dict(scalar_mappable=sm))

(<Figure size 800x800 with 4 Axes>,
 {'pc': <Axes: xlabel='Time [years AD]', ylabel='$PC_2$'>,
  'psd': <Axes: xlabel='Period [years]', ylabel='PSD'>,
  'map': {'cb': <Axes: ylabel='EOF'>,
   'map': <GeoAxes: xlabel='lon', ylabel='lat'>}})

and customize the marker variable:

# Dashboard with marker set to archiveType
res.modeplot(index=1, marker='archiveType', size='elevation', map_kwargs= dict(colorbar=True))

# Dashboard with marker set to archiveType
res.modeplot(index=1, marker='observationType', size='elevation', map_kwargs= dict(colorbar=True))

(<Figure size 800x800 with 5 Axes>,
 {'pc': <Axes: xlabel='Time [years AD]', ylabel='$PC_2$'>,
  'psd': <Axes: xlabel='Period [years]', ylabel='PSD'>,
  'map': {'cb': <Axes: ylabel='EOF'>,
   'map': <GeoAxes: xlabel='lon', ylabel='lat'>,
   'leg': <Axes: >}})

time_geo_plot(figsize=[10, 3], marker=None, markersize=None, alpha=0.8, y_criteria='lat', linestyle=None, linewidth=10, colors=None, cmap='turbo', norm=None, xlabel=None, ylabel=None, title=None, time_unit=None, legend=True, inline_legend=False, plot_kwargs=None, lgd_kwargs=None, label_x_offset=200, label_y_offset=0, savefig_settings=None, ax=None, invert_xaxis=False, invert_yaxis=False)[source]

A plot of the temporal coverage of the records in a MultipleGeoSeries object organized by latitude or longitude.

Similar in behaviour to MultipleSeries.time_coverage_plot

Inspired by Dr. Mara Y. McPartland.

Parameters:

figsize (list, optional) – Size of the figure. The default is [10, 4].
marker (str, optional) – Marker type. The default is None.
markersize (float, optional) – Marker size. The default is None.
alpha (float, optional) – Alpha of the lines
y_criteria (str, optional) – Criteria for the creation of the y_axis. Can be {‘lat’,’lon’,}
linestyle (str, optional) – Line style. The default is None.
linewidth (float, optional) – The width of the line. The default is 10.
colors (a list of, or one, Python supported color code (a string of hex code or a tuple of rgba values)) – Colors for plotting. If None, the plotting will cycle the ‘viridis’ colormap; if only one color is specified, then all curves will be plotted with that single color; if a list of colors are specified, then the plotting will cycle that color list.
cmap (str) – The colormap to use when “colors” is None. Default is ‘turbo’.
norm (matplotlib.colors.Normalize) – The normalization for the colormap. If None, a linear normalization will be used.
xlabel (str, optional) – x-axis label. The default is None.
ylabel (str, optional) – y-axis label. The default is None.
title (str, optional) – Title. The default is None.
time_unit (str) –
the target time unit, possible input: {

’year’, ‘years’, ‘yr’, ‘yrs’, ‘y BP’, ‘yr BP’, ‘yrs BP’, ‘year BP’, ‘years BP’, ‘ky BP’, ‘kyr BP’, ‘kyrs BP’, ‘ka BP’, ‘ka’, ‘my BP’, ‘myr BP’, ‘myrs BP’, ‘ma BP’, ‘ma’,

} default is None, in which case the code picks the most common time unit in the collection. If no unambiguous winner can be found, the unit of the first series in the collection is used.
legend (bool, optional) – Whether the show the legend. The default is True.
inline_legend (bool, optional) – Whether to use inline labels or the default pyleoclim legend. This option overrides lgd_kwargs
plot_kwargs (dict, optional) – Plot parameters. The default is None.
lgd_kwargs (dict, optional) –
Legend parameters. The default is None.

If inline_legend is True, lgd_kwargs will be passed to ax.text() (see matplotlib.axes.Axes.text documentation) If inline_legend is False, lgd_kwargs will be passed to ax.legend() (see matplotlib.axes.Axes.legend documentation)
label_x_offset (float or list, optional) – Amount to offset label by in the x direction. Only used if inline_legend is True. Default is 200. If list, should have the same number of elements as the MultipleSeries object.
label_y_offset (float or list, optional) – Amount to offset label by in the y direction. Only used if inline_legend is True. Default is 0. If list, should have the same number of elements as the MultipleSeries object.
savefig_settings (dictionary, optional) –
the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existing or non-existing path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”} The default is None.
ax (matplotlib.ax, optional) – The matplotlib axis onto which to return the figure. The default is None.
invert_xaxis (bool, optional) – if True, the x-axis of the plot will be inverted
invert_yaxis (bool, optional) – if True, the y-axis of the plot will be inverted

Returns:

fig (matplotlib.figure) – the figure object from matplotlib See [matplotlib.pyplot.figure](https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.figure.html) for details.
ax (matplotlib.axis) – the axis object from matplotlib See [matplotlib.axes](https://matplotlib.org/api/axes_api.html) for details.

See also

pyleoclim.multipleseries.MultipleSeries.time_coverage_plot

pyleoclim.utils.plotting.savefig: Saving figure in Pyleoclim

Examples

from pylipd.utils.dataset import load_dir
lipd = load_dir(name='Pages2k')
df = lipd.get_timeseries_essentials()
dfs = df.query("archiveType in ('Wood','Documents','Coral','Lake sediment')")

# place in a MultipleGeoSeries object
ts_list = []
for _, row in dfs.iloc[:5].iterrows():
    ts_list.append(pyleo.GeoSeries(time=row['time_values'],value=row['paleoData_values'],
                                time_name=row['time_variableName'],value_name=row['paleoData_variableName'],
                                time_unit=row['time_units'], value_unit=row['paleoData_units'],
                                lat = row['geo_meanLat'], lon = row['geo_meanLon'],
                                archiveType = row['archiveType'], verbose = False,
                                label=row['dataSetName']+'_'+row['paleoData_variableName']))

ms = pyleo.MultipleGeoSeries(ts_list, time_unit='years AD')
ms.time_geo_plot()

Loading 16 LiPD files

Loaded..

(<Figure size 1000x300 with 1 Axes>,
 <Axes: xlabel='Time [years AD]', ylabel='Latitude'>)

EnsembleSeries (pyleoclim.EnsembleSeries)

class pyleoclim.core.ensembleseries.EnsembleSeries(series_list, label=None)[source]

EnsembleSeries object

The EnsembleSeries object is a child of the MultipleSeries object, that is, a special case of MultipleSeries, aiming for ensembles of similar series. Ensembles usually arise from age modeling or Bayesian calibrations. All members of an EnsembleSeries object are assumed to share identical labels and units.

All methods available for MultipleSeries are available for EnsembleSeries. Some functions were modified for the special case of ensembles. The class enables ensemble-oriented methods for computation (e.g., quantiles) and visualization (e.g., envelope plot) that are unavailable to other classes.

correlation(target=None, timespan=None, alpha=0.05, method='ttest', statistic='pearsonr', number=1000, settings=None, fdr_kwargs=None, common_time_kwargs=None, mute_pbar=False, seed=None)[source]

Calculate the correlation between an EnsembleSeries object to a target.

If the target is not specified, then the 1st member of the ensemble will be the target Note that the FDR approach is applied by default to determine the significance of the p-values (more information in See Also below).

Parameters:

target (Series or EnsembleSeries) – A pyleoclim Series object or EnsembleSeries object. When the target is also an EnsembleSeries object, then the calculation of correlation is performed in a one-to-one sense, and the ourput list of correlation values and p-values will be the size of the series_list of the self object. That is, if the self object contains n Series, and the target contains n+m Series, then only the first n Series from the object will be used for the calculation; otherwise, if the target contains only n-m Series, then the first m Series in the target will be used twice in sequence.
timespan (tuple) – The time interval over which to perform the calculation
alpha (float) – The significance level (0.05 by default)
method (str, {'ttest','built-in','ar1sim','phaseran','CN'}) – method for significance testing. Default is ‘ttest’ to lower computational cost, but this is not always the best choice
statistic (str) – The name of the statistic used to measure the association, to be chosen from a subset of https://docs.scipy.org/doc/scipy/reference/stats.html#association-correlation-tests Currently supported: [‘pearsonr’,’spearmanr’,’pointbiserialr’,’kendalltau’,’weightedtau’] The default is ‘pearsonr’.
settings (dict) – Parameters for the correlation function (per scipy)
number (int) – the number of simulations (default: 1000)
fdr_kwargs (dict) – Parameters for the FDR function
common_time_kwargs (dict) – Parameters for the method MultipleSeries.common_time()
mute_pbar (bool; {True,False}) – If True, the progressbar will be muted. Default is False.
seed (float or int) – random seed for isopersistent and isospectral methods

Returns:

corr_ens – The resulting object, see pyleoclim.CorrEns

Return type:

CorrEns

See also

pyleoclim.core.series.Series.correlation: pairwise correlations
pyleoclim.utils.correlation.fdr: False Discovery Rate
pyleoclim.core.correns.CorrEns: The correlation ensemble object
pyleoclim.core.surrogateseries.SurrogateSeries: surrogate series

Examples

nn = 50 # number of noise realizations
nt = 100
series_list = []

time, signal = pyleo.utils.gen_ts(model='colored_noise',nt=nt,alpha=2.0)

ts = pyleo.Series(time=time, value = signal, verbose=False).standardize()
noise = np.random.randn(nt,nn)

for idx in range(nn):  # noise
    ts = pyleo.Series(time=time, value=ts.value+5*noise[:,idx], verbose=False)
    series_list.append(ts)

ts_ens = pyleo.EnsembleSeries(series_list)

# to set an arbitrary random seed to fix the result
corr_res = ts_ens.correlation(ts, seed=2333)
print(corr_res)

# to change the statistic:
corr_res = ts_ens.correlation(ts, statistic='kendalltau', method='phaseran', number=20)
print(corr_res)

Looping over 50 Series in the ensemble

  correlation  p-value    signif. w/o FDR (α: 0.05)    signif. w/ FDR (α: 0.05)
-------------  ---------  ---------------------------  --------------------------
129246  0.09       False                        False
171542  0.04       True                         True
196631  0.02       True                         True
204809  0.02       True                         True
215281  0.01       True                         True
311242  < 1e-3     True                         True
389569  < 1e-5     True                         True
426774  < 1e-5     True                         True
471482  < 1e-7     True                         True
491795  < 1e-7     True                         True
462757  < 1e-6     True                         True
500216  < 1e-7     True                         True
477301  < 1e-6     True                         True
483915  < 1e-7     True                         True
533258  < 1e-8     True                         True
597576  < 1e-10    True                         True
60216   < 1e-10    True                         True
600576  < 1e-9     True                         True
574013  < 1e-9     True                         True
600672  < 1e-9     True                         True
622633  < 1e-10    True                         True
654024  < 1e-11    True                         True
657441  < 1e-11    True                         True
675057  < 1e-12    True                         True
708323  < 1e-14    True                         True
691832  < 1e-13    True                         True
704572  < 1e-14    True                         True
732774  < 1e-15    True                         True
747013  < 1e-16    True                         True
759623  < 1e-17    True                         True
775057  < 1e-18    True                         True
79076   < 1e-19    True                         True
795707  < 1e-19    True                         True
81963   < 1e-22    True                         True
826044  < 1e-23    True                         True
859012  < 1e-27    True                         True
87941   < 1e-31    True                         True
878977  < 1e-31    True                         True
889743  < 1e-34    True                         True
881754  < 1e-32    True                         True
899354  < 1e-36    True                         True
915534  < 1e-39    True                         True
933996  < 1e-43    True                         True
951378  < 1e-49    True                         True
957178  < 1e-51    True                         True
957794  < 1e-51    True                         True
965517  < 1e-52    True                         True
980508  < 1e-61    True                         True
989752  < 1e-72    True                         True
       < 1e-6     True                         True
Ensemble size: 50
Looping over 50 Series in the ensemble

  correlation  p-value    signif. w/o FDR (α: 0.05)    signif. w/ FDR (α: 0.05)
-------------  ---------  ---------------------------  --------------------------
0852525  0.35       False                        False
109495   0.10       False                        False
115556   0.10       False                        False
117172   0.05       True                         False
141414   0.05       True                         False
212121   < 1e-6     True                         True
254545   < 1e-6     True                         True
269091   < 1e-6     True                         True
286465   < 1e-6     True                         True
32202    < 1e-6     True                         True
299798   < 1e-6     True                         True
314343   < 1e-6     True                         True
283232   < 1e-6     True                         True
286465   < 1e-6     True                         True
33697    < 1e-6     True                         True
390303   < 1e-6     True                         True
395556   < 1e-6     True                         True
385051   < 1e-6     True                         True
355556   < 1e-6     True                         True
390303   < 1e-6     True                         True
414949   < 1e-6     True                         True
433131   < 1e-6     True                         True
449293   < 1e-6     True                         True
462626   < 1e-6     True                         True
491313   < 1e-6     True                         True
478384   < 1e-6     True                         True
481616   < 1e-6     True                         True
520404   < 1e-6     True                         True
527273   < 1e-6     True                         True
543838   < 1e-6     True                         True
556364   < 1e-6     True                         True
577778   < 1e-6     True                         True
587879   < 1e-6     True                         True
607273   < 1e-6     True                         True
622626   < 1e-6     True                         True
656566   < 1e-6     True                         True
681212   < 1e-6     True                         True
675152   < 1e-6     True                         True
69697    < 1e-6     True                         True
668687   < 1e-6     True                         True
702222   < 1e-6     True                         True
720404   < 1e-6     True                         True
747879   < 1e-6     True                         True
79596    < 1e-6     True                         True
810101   < 1e-6     True                         True
809697   < 1e-6     True                         True
827475   < 1e-6     True                         True
874343   < 1e-6     True                         True
92       < 1e-6     True                         True
        < 1e-6     True                         True
Ensemble size: 50

The print function tabulates the output, and conveys the p-value according to the correlation test applied (“ttest”, by default). To plot the result:

corr_res.plot()

(<Figure size 400x400 with 1 Axes>, <Axes: xlabel='$r$', ylabel='Count'>)

classmethod from_AgeEnsembleArray(series, age_array, value_depth=None, age_depth=None, extrapolate=True, verbose=True)[source]

Function to create an EnsembleSeries object using an age ensemble

Function assumes that the input series and the age array share the same units. If depth vectors are passed, these are also assumed to share the same units

Parameters:

series (pyleoclim.core.series.Series) – A Series object with the values to be mapped
age_array (np.array) – An array of ages to map the values to
value_depth (vector) – An array of depths corresponding to the series values
age_depth (vector) – An array of depths corresponding to the age array
extrapolate (bool) – Whether to extrapolate the age array to the value depth. Default is True
verbose (bool) – Whether to print warnings. Default is True

Returns:

EnsembleSeries – The ensemble created using the time axes from age_array and the values from series.

Return type:

pyleoclim.core.ensembleseries.EnsembleSeries

Examples

#Create an ensemble of 100 series with random time axes of length 1000
length = 1000
age_array = np.array([pyleo.utils.tsmodel.random_time_axis(length) for i in range(100)]).T

#Create a random series
value = np.random.randn(length)
time = pyleo.utils.tsmodel.random_time_axis(length)
series = pyleo.Series(time=time,value=value, verbose=False)

#Create an ensemble using these objects
#Note that the time axis of the series object and the number of rows in the age array must match when depth is not passed
ens = pyleo.EnsembleSeries.from_AgeEnsembleArray(series = series,age_array=age_array,verbose=False)

#If we have depth vectors for our series and age array, we can pass them to the function
age_length = 1000
age_array = np.array([pyleo.utils.tsmodel.random_time_axis(age_length) for i in range(100)]).T
age_depth = np.arange(age_length)

value_length = 800
value = np.random.randn(value_length)
time = pyleo.utils.tsmodel.random_time_axis(value_length)
series = pyleo.Series(time=time, value=value, verbose=False)
value_depth = np.arange(value_length)

#Note that the length of the depth vectors must match the length of the corresponding object (number of values or number of rows in age array)
ens = pyleo.EnsembleSeries.from_AgeEnsembleArray(series = series,age_array=age_array, value_depth=value_depth, age_depth=age_depth,verbose=False)

classmethod from_PaleoEnsembleArray(series, paleo_array, paleo_depth=None, age_depth=None, extrapolate=True, verbose=True)[source]

Function to create an EnsembleSeries object using a paleo ensemble

Function assumes that the input series and the age array share the same units. If depth vectors are passed, these are also assumed to share the same units

Parameters:

series (pyleoclim.core.series.Series) – A Series object with the values to be mapped
paleo_array (np.array) – An array of paleo values to map the values to
paleo_depth (vector) – An array of depths corresponding to the paleo values
age_depth (vector) – An array of depths corresponding to the time axis of the series object
extrapolate (bool) – Whether to extrapolate the age array to the value depth. Default is True
verbose (bool) – Whether to print warnings. Default is True

Returns:

EnsembleSeries – The ensemble created using the time axes from age_array and the values from series.

Return type:

pyleoclim.core.ensembleseries.EnsembleSeries

Examples

#Create an ensemble of 100 series with random time axes of length 1000

length = 1000
paleo_array = np.array([np.random.normal(0, 1, 1000) for _ in range(100)]).T

#Create a random series
value = np.random.normal(0,1,length)
time = pyleo.utils.tsmodel.random_time_axis(length)
series = pyleo.Series(time=time, value=value, verbose=False)

#Create an ensemble using these objects
#Note that the time axis of the series object and the number of rows in the age array must match when depth is not passed
ens = pyleo.EnsembleSeries.from_PaleoEnsembleArray(series = series,paleo_array=paleo_array,verbose=False)

#In this example we'll explicitly pass the depth vectors

paleo_length = 1000
paleo_array = np.array([np.random.normal(0, 1, 1000) for _ in range(100)]).T
paleo_depth = np.arange(paleo_length)

age_length = 800
value = np.random.normal(0,1,age_length)
time = pyleo.utils.tsmodel.random_time_axis(age_length)
series = pyleo.Series(time=time, value=value, verbose=False)
age_depth = np.arange(age_length)

ens = pyleo.EnsembleSeries.from_PaleoEnsembleArray(series = series,paleo_array=paleo_array,paleo_depth=paleo_depth,age_depth=age_depth,verbose=False)

histplot(figsize=[10, 4], title=None, savefig_settings=None, ax=None, ylabel='KDE', vertical=False, edgecolor='w', **plot_kwargs)[source]

Plots the distribution of the timeseries across ensembles

Reuses the seaborn [histplot](https://seaborn.pydata.org/generated/seaborn.histplot.html) function.

Parameters:

figsize (list, optional) – The size of the figure. The default is [10, 4].
title (str, optional) – Title for the figure. The default is None.
savefig_settings (dict, optional) –
the dictionary of arguments for plt.savefig(); some notes below:
- ”path” must be specified; it can be any existed or non-existed path, with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”}.
The default is None.
ax (matplotlib.axis, optional) – A matplotlib axis. The default is None.
ylabel (str, optional) – Label for the count axis. The default is ‘KDE’.
vertical (bool; {True,False}, optional) – Whether to flip the plot vertically. The default is False.
edgecolor (matplotlib.color, optional) – The color of the edges of the bar. The default is ‘w’.
plot_kwargs (dict) – Plotting arguments for seaborn histplot: https://seaborn.pydata.org/generated/seaborn.histplot.html.

See also

pyleoclim.utils.plotting.savefig: Saving figure in Pyleoclim

Examples

nn = 30 # number of noise realizations
nt = 500
series_list = []

time, signal = pyleo.utils.gen_ts(model='colored_noise',nt=nt,alpha=1.0)

ts = pyleo.Series(time=time, value = signal, verbose=False).standardize()
noise = np.random.randn(nt,nn)

for idx in range(nn):  # noise
    ts = pyleo.Series(time=time, value=signal+noise[:,idx], verbose=False)
    series_list.append(ts)

ts_ens = pyleo.EnsembleSeries(series_list)

fig, ax = ts_ens.histplot()

make_labels()[source]

Initialization of labels

Returns:

time_header (str) – Label for the time axis
value_header (str) – Label for the value axis

plot_envelope(figsize=[10, 4], qs=[0.025, 0.25, 0.5, 0.75, 0.975], xlabel=None, ylabel=None, title=None, xlim=None, ylim=None, savefig_settings=None, ax=None, plot_legend=True, curve_clr='#d9544d', curve_lw=2, shade_clr='#d9544d', shade_alpha=0.2, inner_shade_label='IQR', outer_shade_label='95% CI', lgd_kwargs=None)[source]

Plot EnsembleSeries as an envelope.

Parameters:

figsize (list, optional) – The figure size. The default is [10, 4].
qs (list, optional) – The significance levels to consider. The default is [0.025, 0.25, 0.5, 0.75, 0.975] (median, interquartile range, and central 95% region)
xlabel (str, optional) – x-axis label. The default is None.
ylabel (str, optional) – y-axis label. The default is None.
title (str, optional) – Plot title. The default is None.
xlim (list, optional) – x-axis limits. The default is None.
ylim (list, optional) – y-axis limits. The default is None.
savefig_settings (dict, optional) –
the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existed or non-existed path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”} The default is None.
ax (matplotlib.ax, optional) – Matplotlib axis on which to return the plot. The default is None.
plot_legend (bool; {True,False}, optional) – Wether to plot the legend. The default is True.
curve_clr (str, optional) – Color of the main line (median). The default is sns.xkcd_rgb[‘pale red’].
curve_lw (str, optional) – Width of the main line (median). The default is 2.
shade_clr (str, optional) – Color of the shaded envelope. The default is sns.xkcd_rgb[‘pale red’].
shade_alpha (float, optional) – Transparency on the envelope. The default is 0.2.
inner_shade_label (str, optional) – Label for the envelope. The default is ‘IQR’.
outer_shade_label (str, optional) – Label for the envelope. The default is ‘95% CI’.
lgd_kwargs (dict, optional) – Parameters for the legend. The default is None.

Returns:

fig (matplotlib.figure) – the figure object from matplotlib See [matplotlib.pyplot.figure](https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.figure.html) for details.
ax (matplotlib.axis) – the axis object from matplotlib See [matplotlib.axes](https://matplotlib.org/api/axes_api.html) for details.

See also

pyleoclim.utils.plotting.savefig: Saving figure in Pyleoclim

Examples

nn = 30 # number of noise realizations
nt = 500
series_list = []

t,v = pyleo.utils.gen_ts(model='colored_noise',nt=nt,alpha=1.0)
signal = pyleo.Series(time=t,value=v, verbose=False)

for idx in range(nn):  # noise
    noise = np.random.randn(nt,nn)*100
    ts = pyleo.Series(time=signal.time, value=signal.value+noise[:,idx], verbose=False)
    series_list.append(ts)

ts_ens = pyleo.EnsembleSeries(series_list)

fig, ax = ts_ens.plot_envelope(curve_lw=1.5)

plot_traces(figsize=[10, 4], xlabel=None, ylabel=None, title=None, num_traces=10, seed=None, xlim=None, ylim=None, linestyle='-', savefig_settings=None, ax=None, legend=True, color='#d9544d', lw=0.5, alpha=0.3, lgd_kwargs=None)[source]

Plot EnsembleSeries as a subset of traces.

Parameters:

figsize (list, optional) – The figure size. The default is [10, 4].
xlabel (str, optional) – x-axis label. The default is None.
ylabel (str, optional) – y-axis label. The default is None.
title (str, optional) – Plot title. The default is None.
xlim (list, optional) – x-axis limits. The default is None.
ylim (list, optional) – y-axis limits. The default is None.
color (str, optional) – Color of the traces. The default is sns.xkcd_rgb[‘pale red’].
alpha (float, optional) – Transparency of the lines representing the multiple members. The default is 0.3.
linestyle ({'-', '--', '-.', ':', '', (offset, on-off-seq), ...}) – Set the linestyle of the line
lw (float, optional) – Width of the lines representing the multiple members. The default is 0.5.
num_traces (int, optional) – Number of traces to plot. The default is 10.
savefig_settings (dict, optional) –
the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existed or non-existed path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”} The default is None.
ax (matplotlib.ax, optional) – Matplotlib axis on which to return the plot. The default is None.
legend (bool; {True,False}, optional) – Whether to plot the legend. The default is True.
lgd_kwargs (dict, optional) – Parameters for the legend. The default is None.
seed (int, optional) – Set the seed for the random number generator. Useful for reproducibility. The default is None.

Returns:

fig (matplotlib.figure) – the figure object from matplotlib See [matplotlib.pyplot.figure](https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.figure.html) for details.
ax (matplotlib.axis) – the axis object from matplotlib See [matplotlib.axes](https://matplotlib.org/api/axes_api.html) for details.

See also

pyleoclim.utils.plotting.savefig: Saving figure in Pyleoclim

Examples

nn = 30 # number of noise realizations
nt = 500
series_list = []

t,v = pyleo.utils.gen_ts(model='colored_noise',nt=nt,alpha=1.0)
signal = pyleo.Series(time=t,value=v, verbose=False)

for idx in range(nn):  # noise
    noise = np.random.randn(nt,nn)*100
    ts = pyleo.Series(time=signal.time, value=signal.value+noise[:,idx], verbose=False)
    series_list.append(ts)

ts_ens = pyleo.EnsembleSeries(series_list)

fig, ax = ts_ens.plot_traces(alpha=0.2,num_traces=8)

quantiles(qs=[0.05, 0.5, 0.95], axis='value')[source]

Calculate quantiles of an EnsembleSeries object. If axis is ‘value’, the calculation requires for the time axis to be the same. You can use the common_time method to do so. In essence, it transforms the time uncertainty into a y-axis uncertainty. If axis is ‘time’, the values should be the same for all members of the emsemble.

Reuses [scipy.stats.mstats.mquantiles](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.mstats.mquantiles.html) function.

Parameters:

qs (list, optional) – List of quantiles to consider for the calculation. The default is [0.05, 0.5, 0.95].
axis (['time', 'value']) – Whether to calculate the quantiles over the values or time. Default is ‘value’.

Returns:

ens_qs – EnsembleSeries object containing empirical quantiles of original

Return type:

EnsembleSeries

See also

pyleoclim.core.multipleseries.MultipleSeries.common_time: A method to align axes

Examples

nn = 30 # number of noise realizations
nt = 500
series_list = []

t,v = pyleo.utils.gen_ts(model='colored_noise',nt=nt,alpha=1.0)
signal = pyleo.Series(t,v)

for idx in range(nn):  # noise
    noise = np.random.randn(nt,nn)*100
    ts = pyleo.Series(time=signal.time, value=signal.value+noise[:,idx], verbose=False)
    series_list.append(ts)

ts_ens = pyleo.EnsembleSeries(series_list)

ens_qs = ts_ens.quantiles()

Time axis values sorted in ascending order

To calculate in the time dimension:

nn = 30 #number of age models
time = np.arange(1,20000,100) #create a time vector
std_dev = 20 # Noise to be considered

t,v = pyleo.utils.gen_ts(model='colored_noise',nt=len(time),alpha=1.0)

series_list = []

for i in range(nn):
    noise = np.random.normal(0,std_dev,len(time))
    ts=pyleo.Series(time=np.sort(time+noise),value=v,verbose=False)
    series_list.append(ts)

time_ens = pyleo.EnsembleSeries(series_list)

ens_qs = time_ens.quantiles(axis='time')

slice(timespan)[source]

Selects a limited time span from the object

Parameters:: timespan (tuple or list) – The list of time points for slicing, whose length must be even. When there are n time points, the output Series includes n/2 segments. For example, if timespan = [a, b], then the sliced output includes one segment [a, b]; if timespan = [a, b, c, d], then the sliced output includes segment [a, b] and segment [c, d].
Returns:: new – The sliced EnsembleSeries object.
Return type:: EnsembleSeries

Examples

Select part of an object

nn = 20 # number of noise realizations
nt = 200
series_list = []

time, signal = pyleo.utils.gen_ts(model='colored_noise',nt=nt,alpha=2.0)

ts = pyleo.Series(time=time, value = signal, verbose=False).standardize()
noise = np.random.randn(nt,nn)

for idx in range(nn):  # noise
    ts = pyleo.Series(time=time, value=ts.value+5*noise[:,idx], verbose=False)
    series_list.append(ts)

ts_ens = pyleo.EnsembleSeries(series_list)

fig, ax = ts_ens.plot_envelope(curve_lw=1.5)
fig, ax = ts_ens.slice([100, 199]).plot_envelope(curve_lw=1.5)

stackplot(figsize=[5, 15], savefig_settings=None, xlim=None, fill_between_alpha=0.2, colors=None, cmap='tab10', norm=None, spine_lw=1.5, grid_lw=0.5, font_scale=0.8, label_x_loc=-0.15, v_shift_factor=0.75, linewidth=1.5)[source]

Stack plot of multiple series

Note that the plotting style is uniquely designed for this one and cannot be properly reset with pyleoclim.set_style().

Parameters:

figsize (list) – Size of the figure.
colors (list) – Colors for plotting. If None, the plotting will cycle the ‘tab10’ colormap; if only one color is specified, then all curves will be plotted with that single color; if a list of colors are specified, then the plotting will cycle that color list.
cmap (str) – The colormap to use when “colors” is None.
norm (matplotlib.colors.Normalize like) – The nomorlization for the colormap. If None, a linear normalization will be used.
savefig_settings (dictionary) –
the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existed or non-existed path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”} The default is None.
xlim (list) – The x-axis limit.
fill_between_alpha (float) – The transparency for the fill_between shades.
spine_lw (float) – The linewidth for the spines of the axes.
grid_lw (float) – The linewidth for the gridlines.
linewidth (float) – The linewidth for the curves.
font_scale (float) – The scale for the font sizes. Default is 0.8.
label_x_loc (float) – The x location for the label of each curve.
v_shift_factor (float) – The factor for the vertical shift of each axis. The default value 3/4 means the top of the next axis will be located at 3/4 of the height of the previous one.

Returns:

fig (matplotlib.figure) – the figure object from matplotlib See [matplotlib.pyplot.figure](https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.figure.html) for details.
ax (matplotlib.axis) – the axis object from matplotlib See [matplotlib.axes](https://matplotlib.org/api/axes_api.html) for details.

See also

pyleoclim.utils.plotting.savefig: Saving figure in Pyleoclim

Examples

nn = 10 # number of noise realizations
nt = 200
series_list = []

t, v = pyleo.utils.gen_ts(model='colored_noise',nt=nt,alpha=1.0)
signal, _, _ = pyleo.utils.standardize(v)
noise = np.random.randn(nt,nn)

for idx in range(nn):  # noise
    ts = pyleo.Series(time=t, value=signal+noise[:,idx], label='trace #'+str(idx+1), verbose=False)
    series_list.append(ts)

ts_ens = pyleo.EnsembleSeries(series_list)

fig, ax = ts_ens.stackplot()

to_array(axis='value', labels=True)[source]

Returns an ensemble as a numpy array with an optional list for labels. Each column in the array corresponds to an ensemble member.

Parameters:

axis (str, ['time', 'value'], optional) – Whether the return the ensemble from value or time. The default is ‘value’.
labels (bool, [True,False], optional) – Whether to retrun a separate list with the timseries labels. The default is True.

Raises:

ValueError – Axis should be either ‘time’ or ‘value’

Returns:

vals (numpy.array) – An array where each column corresponds to an ensemble member
headers (list) – A list of corresponding labels for each columm

Examples

nn = 30 #number of age models
time = np.arange(1,20000,100) #create a time vector
std_dev = 20 # Noise to be considered

t,v = pyleo.utils.gen_ts(model='colored_noise',nt=len(time),alpha=1.0)

series_list = []

for i in range(nn):
    noise = np.random.normal(0,std_dev,len(time))
    ts=pyleo.Series(time=np.sort(time+noise),value=v,verbose=False)
    series_list.append(ts)

time_ens = pyleo.EnsembleSeries(series_list)
ens_qs = time_ens.quantiles(axis='time')

vals,headers=ens_qs.to_array(axis='time')

to_dataframe(axis='value')[source]

Export the ensemble as a Pandas DataFrame, with members of the ensemble as columns. The columns are labeled according to the label in the individual series or numbered if ‘label’ is None.

Parameters:

axis (str, ['time', 'value']) – Whether the return the ensemble from value or time. each The default is ‘value’.

Raises:

ValueError – Axis should be either ‘time’ or ‘value’

Returns:

df (pandas.DataFrame) – A Pandas DataFrame containing members of the ensemble as columns.
.. jupyter-execute:: – nn = 30 #number of age models time = np.arange(1,20000,100) #create a time vector std_dev = 20 # Noise to be considered

t,v = pyleo.utils.gen_ts(model=’colored_noise’,nt=len(time),alpha=1.0)

series_list = []

for i in range(nn):
noise = np.random.normal(0,std_dev,len(time)) ts=pyleo.Series(time=np.sort(time+noise),value=v,verbose=False) series_list.append(ts)

time_ens = pyleo.EnsembleSeries(series_list) ens_qs = time_ens.quantiles(axis=’time’)

df=ens_qs.to_dataframe(axis=’time’)

EnsembleGeoSeries (pyleoclim.EnsembleGeoSeries)

class pyleoclim.core.ensemblegeoseries.EnsembleGeoSeries(series_list, label=None, lat=None, lon=None, elevation=None, archiveType=None, control_archiveType=False, sensorType=None, observationType=None, depth=None, depth_name=None, depth_unit=None)[source]

EnsembleSeries object

The EnsembleSeries object is a child of the MultipleSeries object, that is, a special case of MultipleSeries, aiming for ensembles of similar series. Ensembles usually arise from age modeling and/or Bayesian calibrations. All members of an EnsembleSeries object are assumed to share identical labels and units.

All methods available for MultipleSeries are available for EnsembleSeries. Some functions were modified for the special case of ensembles. The class enables ensemble-oriented methods for computation (e.g., quantiles) and visualization (e.g., envelope plot) that are unavailable to other classes.

Parameters:

series_list (list) – List of GeoSeries objects
lat (float) – latitude N in decimal degrees. Must be in the range [-90;+90]
lon (float) – longitude East in decimal degrees. Must be in the range [-180;+360] No conversion is applied as mapping utilities convert to [-180,+180] internally
elevation (float) – elevation of the sample, in meters above sea level. Negative numbers indicate depth below global mean sea level, therefore.
label (string) – Name of the time series (e.g., ‘Nino 3.4’) Default is None
log (dict) – Dictionary of tuples documentating the various transformations applied to the object
keep_log (bool) – Whether to keep a log of applied transformations. False by default
importedFrom (string) – source of the dataset. If it came from a LiPD file, this could be the datasetID property
archiveType (string) – climate archive, one of ‘Borehole’, ‘Coral’, ‘FluvialSediment’, ‘GlacierIce’, ‘GroundIce’, ‘LakeSediment’, ‘MarineSediment’, ‘Midden’, ‘MolluskShell’, ‘Peat’, ‘Sclerosponge’, ‘Shoreline’, ‘Speleothem’, ‘TerrestrialSediment’, ‘Wood’ Reference: https://lipdverse.org/vocabulary/archivetype/
control_archiveType ([True, False]) – Whether to standardize the name of the archiveType agains the vocabulary from: https://lipdverse.org/vocabulary/paleodata_proxy/. If set to True, will only allow for these terms and automatically convert known synonyms to the standardized name. Only standardized variable names will be automatically assigned a color scheme. Default is False.
sensorType (string) – sensor, e.g. a paleoclimate proxy sensor. This property can be used to differentiate between species of foraminifera
observationType (string) – observation type, e.g. a proxy observation. See https://lipdverse.org/vocabulary/paleodata_proxy/. Note: this is preferred terminology but not enforced
depth (array) – depth at which the values were collected
depth_name (string) – name of the field, e.g. ‘mid-depth’, ‘top-depth’, etc
depth_unit (string) – units of the depth axis, e.g. ‘cm’

dashboard(figsize=[11, 8], gs=None, plt_kwargs=None, histplt_kwargs=None, spectral_kwargs=None, spectralfig_kwargs=None, map_kwargs=None, hue='archiveType', marker='archiveType', size=None, scatter_kwargs=None, gridspec_kwargs=None, savefig_settings=None)[source]

Dashboard that plots the trace, histogram, map, and power spectrum of the ensemble.

Parameters:

figsize (list or tuple, optional) – Figure size. The default is [11,8].
gs (matplotlib.gridspec object, optional) – Requires at least two rows and 4 columns. - top row, left: timeseries - top row, right: histogram - bottom left: map - bottom right: PSD See [matplotlib.gridspec.GridSpec](https://matplotlib.org/stable/tutorials/intermediate/gridspec.html) for details.

plt_kwargsdict, optional

Optional arguments for the timeseries plot. See Series.plot() or EnsembleSeries.plot_envelope(). The default is None.

histplt_kwargsdict, optional

Optional arguments for the distribution plot. See Series.histplot() or EnsembleSeries.plot_distplot(). The default is None.

spectral_kwargsdict, optional

Optional arguments for the spectral method. Default is to use Lomb-Scargle method. See Series.spectral() or EnsembleSeries.spectral(). The default is None.

spectralfig_kwargsdict, optional

Optional arguments for the power spectrum figure. See PSD.plot() or MultiplePSD.plot_envelope(). The default is None.

map_kwargsdict, optional

Optional arguments for map configuration - projection: str; Optional value for map projection. Default ‘auto’. - proj_default: bool - lakes, land, ocean, rivers, borders, coastline, background: bool or dict; - lgd_kwargs: dict; Optional values for how the map legend is configured - gridspec_kwargs: dict; Optional values for adjusting the arrangement of the colorbar, map and legend in the map subplot - legend: bool; Whether to draw a legend on the figure. Default is True - colorbar: bool; Whether to draw a colorbar on the figure if the data associated with hue are numeric. Default is True The default is None.

huestr, optional

Variable associated with color coding for points plotted on map. May correspond to a continuous or categorical variable. The default is ‘archiveType’.

sizestr, optional

Variable associated with size. Must correspond to a continuous numeric variable. The default is None.

markerstring, optional

Grouping variable that will produce points with different markers. Can have a numeric dtype but will always be treated as categorical. The default is ‘archiveType’.

scatter_kwargsdict, optional

Optional arguments configuring how data are plotted on a map. See description of scatter_kwargs in pyleoclim.utils.mapping.scatter_map

gridspec_kwargsdict, optional

Optional dictionary for configuring dashboard layout using gridspec For information about Gridspec configuration, refer to `Matplotlib documentation <https://matplotlib.org/3.5.0/api/_as_gen/matplotlib.gridspec.GridSpec.html#matplotlib.gridspec.GridSpec>_. The default is None.

savefig_settingsdict, optional

the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existed or non-existed path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”

“format” can be one of {“pdf”, “eps”, “png”, “ps”}.

The default is None.

Returns:

fig (matplotlib.figure) – The figure
ax (dict) – dictionary of matplotlib ax

See also

pyleoclim.core.series.Series.plot: plot a timeseries
pyleoclim.core.ensembleseries.EnsembleSeries.plot_envelope: Envelope plots for an ensemble
pyleoclim.core.series.Series.histplot: plot a distribution of the timeseries
pyleoclim.core.ensembleseries.EnsembleSeries.histplot: plot a distribution of the timeseries across ensembles
pyleoclim.core.series.Series.spectral: spectral analysis method.
pyleoclim.core.multipleseries.MultipleSeries.spectral: spectral analysis method for multiple series.
pyleoclim.core.psds.PSD.signif_test: significance test for timeseries analysis
pyleoclim.core.psds.PSD.plot: plot power spectrum
pyleoclim.core.psds.MulitplePSD.plot: plot envelope of power spectrum
pyleoclim.core.geoseries.GeoSeries.map: map location of dataset
pyleoclim.utils.mapping.scatter_map: Underlying mapping function for Pyleoclim

Examples

ts = pyleo.utils.datasets.load_dataset('EDC-dD')
ts_interp = ts.convert_time_unit('kyr BP').interp(step=.5) # interpolate for a faster result
fig, ax = ts_interp.dashboard()

classmethod from_AgeEnsembleArray(geo_series, age_array, value_depth=None, age_depth=None, extrapolate=True, verbose=True)[source]

Function to create an EnsembleGeoSeries object

Function assumes that the input series and the age array share the same units. If depth vectors are passed, these are also assumed to share the same units. If age_depth is passed and value_depth is not, the function will search for a depth array in the series object.

Parameters:

geo_series (pyleoclim.core.geoseries.GeoSeries) – A Series object with the values to be mapped
age_array (np.array) – An array of ages to map the values to
value_depth (vector) – An array of depths corresponding to the series values
age_depth (vector) – An array of depths corresponding to the age array
extrapolate (bool) – Whether to extrapolate the age array to the value depth. Default is True
verbose (bool) – Whether to print warnings. Default is True

Returns:

EnsembleSeries – The ensemble created using the time axes from age_array and the values from series.

Return type:

pyleoclim.core.ensembleseries.EnsembleSeries

Examples

#Create an ensemble of 100 series with random time axes of length 1000

length = 1000
age_array = np.array([pyleo.utils.tsmodel.random_time_axis(length) for i in range(100)]).T

#Create a random geoseries
value = np.random.randn(length)
time = pyleo.utils.tsmodel.random_time_axis(length)
lat = np.random.randint(-90,90)
lon = np.random.randint(-180,180)
geo_series = pyleo.GeoSeries(time=time, value=value, lat=lat, lon=lon, verbose=False)

#Create an ensemble using these objects
#Note that the time axis of the series object and the number of rows in the age array must match when depth is not passed
ens = pyleo.EnsembleGeoSeries.from_AgeEnsembleArray(geo_series = geo_series,age_array=age_array,verbose=False)

#In this example we'll explicitly pass the depth vectors

age_length = 1000
age_array = np.array([pyleo.utils.tsmodel.random_time_axis(length) for i in range(100)]).T
age_depth = np.arange(age_length)

value_length = 800
value = np.random.normal(0,1,value_length)
time = pyleo.utils.tsmodel.random_time_axis(value_length)
lat = np.random.randint(-90,90)
lon = np.random.randint(-180,180)
geo_series = pyleo.GeoSeries(time=time, value=value, lat=lat, lon=lon, verbose=False)
value_depth = np.arange(value_length)

ens = pyleo.EnsembleGeoSeries.from_AgeEnsembleArray(geo_series = geo_series,age_array=age_array,value_depth=value_depth,age_depth=age_depth,verbose=False)

#We can also pass the depth vector through the GeoSeries object itself

age_length = 1000
age_array = np.array([pyleo.utils.tsmodel.random_time_axis(age_length) for i in range(100)]).T
age_depth = np.arange(age_length)

value_length = 800
value = np.random.randn(value_length)
time = pyleo.utils.tsmodel.random_time_axis(value_length)
lat = np.random.randint(-90,90)
lon = np.random.randint(-180,180)
geo_series = pyleo.GeoSeries(time=time, value=value,depth=np.arange(value_length), lat=lat, lon=lon, verbose=False)

ens = pyleo.EnsembleGeoSeries.from_AgeEnsembleArray(geo_series = geo_series,age_array=age_array, age_depth=age_depth,verbose=False)

classmethod from_PaleoEnsembleArray(geo_series, paleo_array, paleo_depth=None, age_depth=None, extrapolate=True, verbose=True)[source]

Function to create an EnsembleGeoSeries object using a paleo ensemble

Function assumes that the input series and the age array share the same units. If depth vectors are passed, these are also assumed to share the same units

Parameters:

geo_series (pyleoclim.core.series.Series) – A Series object with the values to be mapped
paleo_array (np.array) – An array of paleo values to map the values to
paleo_depth (vector) – An array of depths corresponding to the paleo values
age_depth (vector) – An array of depths corresponding to the time axis of the series object
extrapolate (bool) – Whether to extrapolate the age array to the value depth. Default is True
verbose (bool) – Whether to print warnings. Default is True

Returns:

EnsembleSeries – The ensemble created using the time axes from age_array and the values from series.

Return type:

pyleoclim.core.ensembleseries.EnsembleSeries

Examples

#Create an ensemble of 100 series with random time axes of length 1000

length = 1000
paleo_array = np.array([np.random.normal(0, 1, 1000) for _ in range(100)]).T

#Create a random geoseries
value = np.random.normal(0,1,length)
time = pyleo.utils.tsmodel.random_time_axis(length)
lat = np.random.randint(-90,90)
lon = np.random.randint(-180,180)
geo_series = pyleo.GeoSeries(time=time, value=value, lat=lat, lon=lon, verbose=False)

#Create an ensemble using these objects
#Note that the time axis of the series object and the number of rows in the age array must match when depth is not passed
ens = pyleo.EnsembleGeoSeries.from_PaleoEnsembleArray(geo_series = geo_series,paleo_array=paleo_array,verbose=False)

#In this example we'll explicitly pass the depth vectors

paleo_length = 1000
paleo_array = np.array([np.random.normal(0, 1, 1000) for _ in range(100)]).T
paleo_depth = np.arange(paleo_length)

age_length = 800
value = np.random.normal(0,1,age_length)
time = pyleo.utils.tsmodel.random_time_axis(age_length)
lat = np.random.randint(-90,90)
lon = np.random.randint(-180,180)
geo_series = pyleo.GeoSeries(time=time, value=value, lat=lat, lon=lon, verbose=False)
age_depth = np.arange(age_length)

ens = pyleo.EnsembleGeoSeries.from_PaleoEnsembleArray(geo_series = geo_series,paleo_array=paleo_array,paleo_depth=paleo_depth,age_depth=age_depth,verbose=False)

#Depth can also be passed via the GeoSeries object

paleo_length = 1000
paleo_array = np.array([np.random.normal(0, 1, 1000) for _ in range(100)]).T
paleo_depth = np.arange(paleo_length)

age_length = 800
value = np.random.normal(0,1,age_length)
time = pyleo.utils.tsmodel.random_time_axis(age_length)
lat = np.random.randint(-90,90)
lon = np.random.randint(-180,180)
depth = np.arange(age_length)
geo_series = pyleo.GeoSeries(time=time, value=value, depth=depth,lat=lat, lon=lon, verbose=False)

ens = pyleo.EnsembleGeoSeries.from_PaleoEnsembleArray(geo_series = geo_series,paleo_array=paleo_array,paleo_depth=paleo_depth,verbose=False)

make_labels()[source]

Initialization of labels

Returns:

time_header (str) – Label for the time axis
value_header (str) – Label for the value axis

MulEnsGeoSeries (pyleoclim.MulEnsGeoSeries)

class pyleoclim.core.mulensgeoseries.MulEnsGeoSeries(ensemble_series_list, label=None)[source]

mcpca(nsim=1000, seed=None, common_time_kwargs=None, pca_kwargs=None, align_method='dot', delta=0, phase_max=0.5)[source]

Function to conduct Monte-Carlo PCA using ensembles included in this object

Parameters:

nsim (int) – Number of simulations to carry out. Default is 1000
seed (int) – Seed to use for random calculations
common_time_kwargs (dict) – Key word arguments for MultipleSeries.common_time()
pca_kwargs (dict) – Key word arguments for MultipleGeoSeries.pca()
align_method (str; {'correlation','dot','phase','cosine'}) – How to align pcs. Correlation computes the correlation between individual principal components and the first set of PCs, flipping those that have a negative correlation. Dot computes the dot product between the eigenvectors, and flips those that have a negative dot product. Phase computes the phase difference between the eigenvectors, and flips those that have a phase difference greater than phase_max. Cosine computes the cosine similarity between the eigenvectors, and flips those that have a negative cosine similarity.
delta (float) – A number between -1 and 1, that serves as the threshold for the correlation, dot, and cosine methods. Default is 0.
phase_max (float) – A float between 0 and 1, that serves as the threshold for the phase method. Default is 0.5.

Returns:

EnsembleMvD – Ensemble Multivariate Decomposition object

Return type:

pyleo.EnsMultivarDecomp

Examples

n = 3 # number of ensembles
nn = 30 # number of noise realizations
nt = 500
ens_list = []

t,v = pyleo.utils.gen_ts(model='colored_noise',nt=nt,alpha=1.0)
signal = pyleo.Series(t,v)

for _ in range(n):
    series_list = []
    lat = np.random.randint(-90,90)
    lon = np.random.randint(-180,180)
    for idx in range(nn):  # noise
        noise = np.random.randn(nt,nn)*100
        ts = pyleo.GeoSeries(time=signal.time, value=signal.value+noise[:,idx], lat=lat, lon=lon, verbose=False)
        series_list.append(ts)

    ts_ens = pyleo.EnsembleGeoSeries(series_list)
    ens_list.append(ts_ens)

mul_ens = pyleo.MulEnsGeoSeries(ens_list)
mul_ens.mcpca(nsim=10,seed=42)

Time axis values sorted in ascending order

<pyleoclim.core.ensmultivardecomp.EnsMultivarDecomp at 0x715015b5e3f0>

stackplot(figsize=None, savefig_settings=None, time_unit=None, xlim=None, colors=None, cmap='tab10', plot_style='envelope', norm=None, labels='auto', ylabel_fontsize=8, spine_lw=1.5, grid_lw=0.5, label_x_loc=-0.15, v_shift_factor=0.75, yticks_minor=False, xticks_minor=False, ylims='auto', plot_kwargs=None, common_time_kwargs=None)[source]

Stack plot of multiple ensemble series

Time units are harmonized prior to plotting. Functionally, this method is very similar to the stackplot method of MultipleSeries, see the documentation there for more details on customization. Note that the plotting plot_style is uniquely designed for this one and cannot be properly reset with pyleoclim.set_plot_style().

Parameters:

figsize (list) – Size of the figure.
savefig_settings (dictionary) –
the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existing or non-existing path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”} The default is None.
time_unit (str) –
the target time unit, possible inputs: {

’year’, ‘years’, ‘yr’, ‘yrs’, ‘y BP’, ‘yr BP’, ‘yrs BP’, ‘year BP’, ‘years BP’, ‘ky BP’, ‘kyr BP’, ‘kyrs BP’, ‘ka BP’, ‘ka’, ‘my BP’, ‘myr BP’, ‘myrs BP’, ‘ma BP’, ‘ma’,

} default is None, in which case the code picks the most common time unit in the collection. If no discernible winner can be found, the unit of the first series in the collection is used.
xlim (list) – The x-axis limit.
colors (a list of, or one, Python supported color code (a string of hex code or a tuple of rgba values)) – Colors for plotting. If None, the plotting will cycle the ‘tab10’ colormap; if only one color is specified, then all curves will be plotted with that single color; if a list of colors are specified, then the plotting will cycle that color list.
cmap (str) – The colormap to use when “colors” is None. Note that the function will try to detect continuous or discrete colormaps, and set the norm accordingly.
plot_style (str; {'envelope', 'traces'}) – The ensemble plotting style to use. Default is ‘envelope’.
norm (matplotlib.colors.Normalize like) – The normalization for the colormap. If None, a linear normalization will be used.
labels (None, 'auto' or list) – If None, doesn’t add labels to the subplots If ‘auto’, uses the labels passed during the creation of pyleoclim.Series If list, pass a list of strings for each labels. Default is ‘auto’
spine_lw (float) – The linewidth for the spines of the axes.
grid_lw (float) – The linewidth for the gridlines.
label_x_loc (float) – The x location for the label of each curve.
v_shift_factor (float) – The factor for the vertical shift of each axis. The default value 3/4 means the top of the next axis will be located at 3/4 of the height of the previous one.
ylabel_fontsize (int) – Size for ylabel font. Default is 8, to avoid crowding.
yticks_minor (bool) – Whether the y axes should contain minor ticks (use sparingly!). Default: False
xticks_minor (bool) – Whether the x axis should contain minor ticks. Default: False
ylims (str {'spacious', 'auto'}) – Method for determining the limits of the y axes. Default is ‘spacious’, which is mean +/- 4 x std ‘auto’ activates the Matplotlib default
plot_kwargs (dict or list of dict) –
Arguments to further customize the plot from EnsembleSeries.plot_envelope or EnsembleSeries.plot_traces, depending on the chosen style.
- Dictionary: Arguments will be applied to all lines in the stackplots
- List of dictionaries: Allows to customize one line at a time.
common_time_kwargs (dict) – Arguments to pass to the common_time method of the ensemble series. Common time is called to calculate the median of the ensemble series for tick purposes, and is also used if plot_style is set to ‘envelope’.

Returns:

fig (matplotlib.figure) – the figure object from matplotlib See [matplotlib.pyplot.figure](https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.figure.html) for details.
ax (matplotlib.axis) – the axis object from matplotlib See [matplotlib.axes](https://matplotlib.org/api/axes_api.html) for details.

See also

pyleoclim.core.multipleseries.MultipleSeries.stackplot: Stack plot of multiple series
pyleoclim.core.ensembleseries.EnsembleSeries.plot_envelope: Plotting the envelope of an ensemble of series
pyleoclim.core.ensembleseries.EnsembleSeries.plot_traces: Plotting the traces of an ensemble of series
pyleoclim.utils.plotting.savefig: Saving figure in Pyleoclim

Examples

n = 3 # number of ensembles
nn = 30 # number of noise realizations
nt = 500
ens_list = []

t,v = pyleo.utils.gen_ts(model='colored_noise',nt=nt,alpha=1.0)
signal = pyleo.Series(t,v)

for _ in range(n):
    series_list = []
    lat = np.random.randint(-90,90)
    lon = np.random.randint(-180,180)
    for idx in range(nn):  # noise
        noise = np.random.randn(nt,nn)*100
        ts = pyleo.GeoSeries(time=signal.time, value=signal.value+noise[:,idx], lat=lat, lon=lon, verbose=False)
        series_list.append(ts)

    ts_ens = pyleo.EnsembleGeoSeries(series_list)
    ens_list.append(ts_ens)

mul_ens = pyleo.MulEnsGeoSeries(ens_list)
mul_ens.stackplot()

Time axis values sorted in ascending order

(<Figure size 640x480 with 4 Axes>,
 {0: <Axes: ylabel='value'>,
  1: <Axes: ylabel='value'>,
  2: <Axes: ylabel='value'>,
  'x_axis': <Axes: xlabel='Time [years CE]'>})

If you’d like to adjust the plot parameters, you can pass them via the plot_kwargs argument.

n = 3 # number of ensembles
nn = 30 # number of noise realizations
nt = 500
ens_list = []

t,v = pyleo.utils.gen_ts(model='colored_noise',nt=nt,alpha=1.0)
signal = pyleo.Series(t,v)

for _ in range(n):
    series_list = []
    lat = np.random.randint(-90,90)
    lon = np.random.randint(-180,180)
    for idx in range(nn):  # noise
        noise = np.random.randn(nt,nn)*100
        ts = pyleo.GeoSeries(time=signal.time, value=signal.value+noise[:,idx], lat=lat, lon=lon, verbose=False)
        series_list.append(ts)

    ts_ens = pyleo.EnsembleGeoSeries(series_list)
    ens_list.append(ts_ens)

mul_ens = pyleo.MulEnsGeoSeries(ens_list)
mul_ens.stackplot(plot_kwargs={'shade_alpha':0.5})

Time axis values sorted in ascending order

(<Figure size 640x480 with 4 Axes>,
 {0: <Axes: ylabel='value'>,
  1: <Axes: ylabel='value'>,
  2: <Axes: ylabel='value'>,
  'x_axis': <Axes: xlabel='Time [years CE]'>})

If you’d like to plot traces instead, you can use the modify the plot_style argument

n = 3 # number of ensembles
nn = 30 # number of noise realizations
nt = 500
ens_list = []

t,v = pyleo.utils.gen_ts(model='colored_noise',nt=nt,alpha=1.0)
signal = pyleo.Series(t,v)

for _ in range(n):
    series_list = []
    lat = np.random.randint(-90,90)
    lon = np.random.randint(-180,180)
    for idx in range(nn):  # noise
        noise = np.random.randn(nt,nn)*100
        ts = pyleo.GeoSeries(time=signal.time, value=signal.value+noise[:,idx], lat=lat, lon=lon, verbose=False)
        series_list.append(ts)

    ts_ens = pyleo.EnsembleGeoSeries(series_list)
    ens_list.append(ts_ens)

mul_ens = pyleo.MulEnsGeoSeries(ens_list)
mul_ens.stackplot(plot_style='traces')

Time axis values sorted in ascending order

(<Figure size 640x480 with 4 Axes>,
 {0: <Axes: ylabel='value'>,
  1: <Axes: ylabel='value'>,
  2: <Axes: ylabel='value'>,
  'x_axis': <Axes: xlabel='Time [years CE]'>})

SurrogateSeries (pyleoclim.SurrogateSeries)

class pyleoclim.core.surrogateseries.SurrogateSeries(series_list=[], number=1, method=None, label=None, param=None, seed=None)[source]

Object containing surrogate timeseries, obtained either by emulating an existing series (see from_series()) or from a parametric model (see from_param())

Surrogate Series is a child of EnsembleSeries. All methods available for EnsembleSeries are available for surrogate series.

from_param(param, length=50, time_pattern='even', settings=None)[source]

Simulate the SurrogateSeries from a parametric model

Parameters:

param (list) – model parameters (e.g. [tau, sigma0] for an AR(1) model, [beta] for colored noise)
length (int) – Length of the series. Default: 50
time_pattern (str {even, random, specified}) –
The pattern used to generate the surrogate time axes * ‘even’ uses an evenly-spaced time with spacing delta_t specified in settings (if not specified, defaults to 1.0) * ‘random’ uses random_time_axis() with specified distribution and parameters

Relevant settings are delta_t_dist and param. (default: ‘exponential’ with parameter = 1.0)
- ’specified’: uses time axis time from settings.
seed (int) – Control random seed option for reproducibility
settings (dict) – Parameters for surrogate generator. See individual methods for details.

Returns:

surr

Return type:

SurrogateSeries

See also

pyleoclim.utils.tsmodel.ar1_sim: AR(1) simulator
pyleoclim.utils.tsmodel.colored_noise: simulating from a power law spectrum, $S(f) propto f^{-eta}$
pyleoclim.utils.tsmodel.random_time_axis: Generate time increment vector according to a specific probability model

Examples

ar1 = pyleo.SurrogateSeries(method='ar1sim', number=10)
ar1.from_param(length=100, param = [2,2])
ar1.plot_envelope(title=rf'AR(1) synthetic series ($        au={2},\sigma^2={2}$)')

(<Figure size 1000x400 with 1 Axes>,
 <Axes: title={'center': 'AR(1) synthetic series ($        au=2,\\sigma^2=2$)'}, xlabel='Time [years CE]', ylabel='value'>)

from_series(target_series)[source]

Fashion the SurrogateSeries object after a target series

Parameters:

target_series (Series object) – target Series used to infer surrogate properties
seed (int) – Control random seed option for reproducibility

Returns:

surr

Return type:

SurrogateSeries

See also

pyleoclim.utils.tsmodel.ar1_sim: AR(1) simulator
pyleoclim.utils.tsmodel.uar1_sim: maximum likelihood AR(1) simulator
pyleoclim.utils.tsutils.phaseran: phase randomization

Examples

SOI = pyleo.utils.load_dataset(‘SOI’).interp() SOI_surr = pyleo.SurrogateSeries(method=’phaseran’, number=4) SOI_surr.from_series(SOI) fig, ax = SOI_surr.plot_traces() SOI.plot(ax=ax,color=’black’,linewidth=1, ylim=[-6,3]) ax.legend(loc=’lower center’, bbox_to_anchor=(0.5, 0), ncol=2) ax.set_title(SOI_surr.label)

PSD (pyleoclim.PSD)

class pyleoclim.core.psds.PSD(frequency, amplitude, label=None, timeseries=None, plot_kwargs=None, spec_method=None, spec_args=None, signif_qs=None, signif_method=None, period_unit=None, beta_est_res=None)[source]

The PSD (Power spectral density) class is intended for conveniently manipulating the result of spectral methods, including performing significance tests, estimating scaling coefficients, and plotting.

See examples in pyleoclim.core.series.Series.spectral to see how to create and manipulate these objects

Parameters:

frequency (numpy.array, list, or float) – One or more frequencies in power spectrum
amplitude (numpy.array, list, or float) – The amplitude at each (frequency, time) point; note the dimension is assumed to be (frequency, time)
label (str, optional) – Descriptor of the PSD. Default is None
timeseries (pyleoclim.Series, optional) – Default is None
plot_kwargs (dict, optional) – Plotting arguments for seaborn histplot: https://seaborn.pydata.org/generated/seaborn.histplot.html. Default is None
spec_method (str, optional) – The name of the spectral method to be applied on the timeseries Default is None
spec_args (dict, optional) – Arguments for wavelet analysis (‘freq’, ‘scale’, ‘mother’, ‘param’) Default is None
signif_qs (pyleoclim.MultipleScalogram, optional) – Pyleoclim MultipleScalogram object containing the quantiles qs of the surrogate scalogram distribution. Default is None
signif_method (str, optional) – The method used to obtain the significance level. Default is None
period_unit (str, optional) – Unit of time. Default is None
beta_est_res (list or numpy.array, optional) – Results of the beta estimation calculation. Default is None.

See also

pyleoclim.core.series.Series.spectral: Spectral analysis
pyleoclim.core.scalograms.Scalogram: Scalogram object
pyleoclim.core.scalograms.MultipleScalogram: Object storing multiple scalogram objects
pyleoclim.core.psds.MultiplePSD: Object storing several PSDs from different Series or ensemble members in an age model

anti_alias(avgs=2)[source]

Apply the anti-aliasing filter

Parameters:: avgs (int) – flag for whether spectrum is derived from instantaneous point measurements (avgs<>1) OR from measurements averaged over each sampling interval (avgs==1)
Returns:: new – New PSD object with the spectral aliasing effect alleviated.
Return type:: pyleoclim.core.psds.PSD

Examples

Generate colored noise with scaling exponent equals to unity, and test the impact of anti-aliasing filter

t, v = pyleo.utils.tsmodel.gen_ts('colored_noise', alpha=1, m=1e5) # m=1e5 leads to aliasing
ts = pyleo.Series(time=t, value=v, label='colored noise', verbose=False)

# without the anti-aliasing filter
fig, ax = ts.spectral(method='mtm').beta_est().plot()

# with the anti-aliasing filter
fig, ax = ts.spectral(method='mtm').anti_alias().beta_est().plot()

References

Kirchner, J. W. Aliasing in 1/f(alpha) noise spectra: origins, consequences, and remedies. Phys Rev E Stat Nonlin Soft Matter Phys 71, 66110 (2005).

See also

pyleoclim.utils.wavelet.AliasFilter.alias_filter: anti-aliasing filter

beta_est(fmin=None, fmax=None, logf_binning_step='max', verbose=False)[source]

Estimate the scaling exponent (beta) of the PSD

For a power law S(f) ~ f^beta in log-log space, beta is simply the slope.

Parameters:

fmin (float, optional) – the minimum frequency edge for beta estimation; the default is the minimum of the frequency vector of the PSD obj
fmax (float, optional) – the maximum frequency edge for beta estimation; the default is the maximum of the frequency vector of the PSD obj
logf_binning_step (str, {'max', 'first'}) – if ‘max’, then the maximum spacing of log(f) will be used as the binning step if ‘first’, then the 1st spacing of log(f) will be used as the binning step
verbose (bool; {True, False}) – If True, will print warning messages if there is any

Returns:

new – New PSD object with the estimated scaling slope information, which is stored as a dictionary that includes: - beta: the scaling factor - std_err: the one standard deviation error of the scaling factor - f_binned: the binned frequency series, used as X for linear regression - psd_binned: the binned PSD series, used as Y for linear regression - Y_reg: the predicted Y from linear regression, used with f_binned for the slope curve plotting

Return type:

pyleoclim.core.psds.PSD

Examples

Generate fractal noise and verify that its scaling exponent is close to unity

t, v = pyleo.utils.tsmodel.gen_ts(model='colored_noise')
ts = pyleo.Series(time=t, value= v, label = 'fractal noise, unit slope', verbose=False)
psd = ts.detrend().spectral(method='cwt')

# estimate the scaling slope
psd_beta = psd.beta_est(fmin=1/50, fmax=1/2)

fig, ax = psd_beta.plot(color='tab:blue',beta_kwargs={'color':'tab:red','linewidth':2})

See also

pyleoclim.core.series.Series.spectral: spectral analysis
pyleoclim.utils.spectral.beta_estimation: Estimate the scaling exponent of a power spectral density
pyleoclim.core.psds.PSD.plot: plotting method for PSD objects

copy()[source]: Copy object

plot(in_loglog=True, in_period=True, label=None, xlabel=None, ylabel='PSD', title=None, marker=None, markersize=None, color=None, linestyle=None, linewidth=None, transpose=False, xlim=None, ylim=None, figsize=[10, 4], savefig_settings=None, ax=None, legend=True, lgd_kwargs=None, xticks=None, yticks=None, alpha=None, zorder=None, plot_kwargs=None, signif_clr='red', signif_linestyles=['--', '-.', ':'], signif_linewidth=1, plot_beta=True, beta_kwargs=None)[source]

Plots the PSD estimates and signif level if included

Parameters:

in_loglog (bool; {True, False}, optional) – Plot on loglog axis. The default is True.
in_period (bool; {True, False}, optional) – Plot the x-axis as periodicity rather than frequency. The default is True.
label (str, optional) – label for the series. The default is None.
xlabel (str, optional) – Label for the x-axis. The default is None. Will guess based on Series
ylabel (str, optional) – Label for the y-axis. The default is ‘PSD’.
title (str, optional) – Plot title. The default is None.
marker (str, optional) – marker to use. The default is None.
markersize (int, optional) – size of the marker. The default is None.
color (str, optional) – Line color. The default is None.
linestyle (str, optional) – linestyle. The default is None.
linewidth (float, optional) – Width of the line. The default is None.
transpose (bool; {True, False}, optional) – Plot periodicity on y-. The default is False.
xlim (list, optional) – x-axis limits. The default is None.
ylim (list, optional) – y-axis limits. The default is None.
figsize (list, optional) – Figure size. The default is [10, 4].
savefig_settings (dict, optional) –
save settings options. The default is None. the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existed or non-existed path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”}
ax (ax, optional) – The matplotlib.Axes object onto which to return the plot. The default is None.
legend (bool; {True, False}, optional) – whether to plot the legend. The default is True.
lgd_kwargs (dict, optional) – Arguments for the legend. The default is None.
xticks (list, optional) – xticks to use. The default is None.
yticks (list, optional) – yticks to use. The default is None.
alpha (float, optional) – Transparency setting. The default is None.
zorder (int, optional) – Order for the plot. The default is None.
plot_kwargs (dict, optional) – Other plotting argument. The default is None.
signif_clr (str, optional) – Color for the significance line. The default is ‘red’.
signif_linestyles (list of str, optional) – Linestyles for significance. The default is [’–’, ‘-.’, ‘:’].
signif_linewidth (float, optional) – width of the significance line. The default is 1.
plot_beta (bool; {True, False}, optional) – If True and self.beta_est_res is not None, then the scaling slope line will be plotted
beta_kwargs (dict, optional) – The visualization keyword arguments for the scaling slope

Return type:

fig, ax

Examples

Generate fractal noise, assess significance against an AR(1) benchmark, and plot:

import matplotlib.pyplot as plt

t, v = pyleo.utils.tsmodel.gen_ts(model='colored_noise')
ts = pyleo.Series(time = t, value = v, label = 'fractal noise', verbose=False)
tsn = ts.standardize()

psd_sim = tsn.spectral(method='mtm').signif_test(number=20)
psd_sim.plot()

(<Figure size 1000x400 with 1 Axes>,
 <Axes: xlabel='Period [years]', ylabel='PSD'>)

If you add the estimate of the scaling exponent, the line of best fit will be added to the plot, and the estimated exponent to its legend. For instance:

psd_beta = psd_sim.beta_est(fmin=1/100, fmax=1/2)
fig, ax = psd_beta.plot()

See also

pyleoclim.core.series.Series.spectral: spectral analysis
pyleoclim.core.psds.PSD.signif_test: significance testing for PSD objects
pyleoclim.core.psds.PSD.beta_est: scaling exponent estimation for PSD objects

signif_test(method='ar1sim', number=None, seed=None, qs=[0.95], settings=None, scalogram=None)[source]

Parameters:

number (int, optional) – Number of surrogate series to generate for significance testing. The default is None.
method (str; {'ar1asym','ar1sim','uar1'}) – Method to generate surrogates. AR1sim uses simulated timeseries with similar persistence. AR1asymp represents the closed form solution. The default is AR1sim
seed (int, optional) – Option to set the seed for reproducibility. The default is None.
qs (list, optional) – Significance levels to return. The default is [0.95].
settings (dict, optional) – Parameters for the specific significance test. The default is None. Note that the default value for the asymptotic solution is time-average
scalogram (pyleoclim.Scalogram object, optional) – Scalogram containing signif_scals exported during significance testing of scalogram. If number is None and signif_scals are present, will use length of scalogram list as number of significance tests

Returns:

new – New PSD object with appropriate significance test

Return type:

pyleoclim.core.psds.PSD

Examples

Compute the spectrum of the Southern Oscillation Index and assess significance against an AR(1) benchmark:

nino3 = pyleo.utils.load_dataset('NINO3')
psd = nino3.standardize().spectral('mtm',settings={'NW':2})
psd_sim = psd.signif_test(number=20)
fig, ax = psd_sim.plot()

By default, this method uses 200 Monte Carlo simulations of an AR(1) process. For a smoother benchmark, up the number of simulations. Also, you may obtain and visualize several quantiles at once, e.g. 90% and 95%:

psd_1000 = psd.signif_test(number=100, qs=[0.90, 0.95])
fig, ax = psd_1000.plot()

Another option is to use a closed-form, asymptotic solution for the AR(1) spectrum:

psd_asym = psd.signif_test(method='ar1asym',qs=[0.90, 0.95])
fig, ax = psd_asym.plot()

If significance tests from a comparable scalogram have been saved, they can be passed here to speed up the generation of noise realizations for significance testing. Setting export_scal to True saves the noise realizations generated during significance testing for future use:

scalogram = nino3.standardize().wavelet().signif_test(number=20, export_scal=True)

The psd can be calculated by using the previously generated scalogram

psd_scal = nino3.standardize().spectral(scalogram=scalogram)

The same scalogram can then be passed to do significance testing. Pyleoclim will dig through the scalogram object to find the saved noise realizations and reuse them flexibly.

fig, ax = psd.signif_test(scalogram=scalogram).plot()

See also

pyleoclim.utils.wavelet.tc_wave_signif: asymptotic significance calculation
pyleoclim.core.psds.MultiplePSD: Object storing several PSDs from different Series or ensemble members in an age model
pyleoclim.core.scalograms.Scalogram: Scalogram object
pyleoclim.core.series.Series.surrogates: Generate surrogates with increasing time axis
pyleoclim.core.series.Series.spectral: Performs spectral analysis on Pyleoclim Series
pyleoclim.core.series.Series.wavelet: Performs wavelet analysis on Pyleoclim Series

MultiplePSD (pyleoclim.MultiplePSD)

class pyleoclim.core.psds.MultiplePSD(psd_list, beta_est_res=None)[source]
MultiplePSD objects store several PSDs from different Series or ensemble members from a posterior distribution (e.g. age model, Bayesian climate reconstruction, etc). This is used extensively for Monte Carlo significance tests.

anti_alias(avgs=2, mute_pbar=False)[source]

Apply the anti-aliasing filter

Parameters:

avgs (int) – flag for whether spectrum is derived from instantaneous point measurements (avgs<>1) OR from measurements averaged over each sampling interval (avgs==1)

mute_pbar (bool; {True,False}) – If True, the progressbar will be muted. Default is False.

Returns:

new – New MultiplePSD object with the spectral aliasing effect alleviated.

Return type:

pyleoclim.core.psds.MultiplePSD

References

Kirchner, J. W. Aliasing in 1/f(alpha) noise spectra: origins, consequences, and remedies. Phys Rev E Stat Nonlin Soft Matter Phys 71, 66110 (2005).

See also

pyleoclim.utils.wavelet.AliasFilter.alias_filter
anti-aliasing filter

beta_est(fmin=None, fmax=None, logf_binning_step='max', verbose=False)[source]

Estimate the scaling exponent of each constituent PSD

This function calculates the scaling exponent (beta) for each of the PSDs stored in the object. The scaling exponent represents the slope of the spectrum in log-log space.

Parameters:

fmin (float) – the minimum frequency edge for beta estimation; the default is the minimum of the frequency vector of the PSD object

fmax (float) – the maximum frequency edge for beta estimation; the default is the maximum of the frequency vector of the PSD object

logf_binning_step (str; {'max', 'first'}) – if ‘max’, then the maximum spacing of log(f) will be used as the binning step. if ‘first’, then the 1st spacing of log(f) will be used as the binning step.

verbose (bool) – If True, will print warning messages if there is any

Returns:

new – New MultiplePSD object with the estimated scaling slope information, which is stored as a dictionary that includes: - beta: the scaling factor - std_err: the one standard deviation error of the scaling factor - f_binned: the binned frequency series, used as X for linear regression - psd_binned: the binned PSD series, used as Y for linear regression - Y_reg: the predicted Y from linear regression, used with f_binned for the slope curve plotting

Return type:

pyleoclim.core.psds.MultiplePSD

See also

pyleoclim.core.psds.PSD.beta_est
scaling exponent estimation for a single PSD object

copy()[source]

Copy object

plot(figsize=[10, 4], in_loglog=True, in_period=True, xlabel=None, ylabel='PSD', title=None, xlim=None, ylim=None, savefig_settings=None, ax=None, xticks=None, yticks=None, legend=True, colors=None, cmap=None, norm=None, plot_kwargs=None, lgd_kwargs=None)[source]

Plot multiple PSDs on the same plot

Parameters:

figsize (list, optional) – Figure size. The default is [10, 4].

in_loglog (bool, optional) – Whether to plot in loglog. The default is True.

in_period (bool, {True, False} optional) – Plots against periods instead of frequencies. The default is True.

xlabel (str, optional) – x-axis label. The default is None.

ylabel (str, optional) – y-axis label. The default is ‘Amplitude’.

title (str, optional) – Title for the figure. The default is None.

xlim (list, optional) – Limits for the x-axis. The default is None.

ylim (list, optional) – limits for the y-axis. The default is None.

colors (a list of, or one, Python supported color code (a string of hex code or a tuple of rgba values)) – Colors for plotting. If None, the plotting will cycle the ‘tab10’ colormap; if only one color is specified, then all curves will be plotted with that single color; if a list of colors are specified, then the plotting will cycle that color list.

cmap (str) – The colormap to use when “colors” is None.

norm (matplotlib.colors.Normalize like) – The nomorlization for the colormap. If None, a linear normalization will be used.

savefig_settings (dict, optional) –
the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existing or non-existing path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”

”format” can be one of {“pdf”, “eps”, “png”, “ps”}

ax (matplotlib axis, optional) – The matplotlib axis object on which to retrun the figure. The default is None.

xticks (list, optional) – x-ticks label. The default is None.

yticks (list, optional) – y-ticks label. The default is None.

legend (bool, optional) – Whether to plot the legend. The default is True.

plot_kwargs (dictionary, optional) – Parameters for plot function. The default is None.

lgd_kwargs (dictionary, optional) – Parameters for legend. The default is None.

Returns:

fig (matplotlib.pyplot.figure)

ax (matplotlib.pyplot.axis)

See also

pyleoclim.core.psds.PSD.plot
plotting method for PSD objects
plot_envelope(figsize=[10, 4], qs=[0.025, 0.5, 0.975], in_loglog=True, in_period=True, xlabel=None, ylabel='PSD', title=None, xlim=None, ylim=None, savefig_settings=None, ax=None, xticks=None, yticks=None, plot_legend=True, curve_clr='#d9544d', curve_lw=3, shade_clr='#d9544d', shade_alpha=0.3, shade_label=None, lgd_kwargs=None, members_plot_num=10, members_alpha=0.3, members_lw=1, seed=None)[source]
Plot an envelope statistics for mulitple PSD

This function plots an envelope statistics from multiple PSD. This is especially useful when the PSD are coming from an ensemble of possible solutions (e.g., age ensembles)

Parameters:

figsize (list, optional) – The figure size. The default is [10, 4].

qs (list, optional) – The significance levels to consider. The default is [0.025, 0.5, 0.975].

in_loglog (bool, optional) – Plot in log space. The default is True.

in_period (bool, optional) – Whether to plot periodicity instead of frequency. The default is True.

xlabel (str, optional) – x-axis label. The default is None.

ylabel (str, optional) – y-axis label. The default is ‘Amplitude’.

title (str, optional) – Plot title. The default is None.

xlim (list, optional) – x-axis limits. The default is None.

ylim (list, optional) – y-axis limits. The default is None.

savefig_settings (dict, optional) –
the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existing or non-existing path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”

”format” can be one of {“pdf”, “eps”, “png”, “ps”} The default is None.

ax (matplotlib.ax, optional) – Matplotlib axis on which to return the plot. The default is None.

xticks (list, optional) – xticks label. The default is None.

yticks (list, optional) – yticks label. The default is None.

plot_legend (bool, optional) – Whether to plot the legend. The default is True.

curve_clr (str, optional) – Color of the main PSD. The default is sns.xkcd_rgb[‘pale red’].

curve_lw (str, optional) – Width of the main PSD line. The default is 3.

shade_clr (str, optional) – Color of the shaded envelope. The default is sns.xkcd_rgb[‘pale red’].

shade_alpha (float, optional) – Transparency on the envelope. The default is 0.3.

shade_label (str, optional) – Label for the envelope. The default is None.

lgd_kwargs (dict, optional) – Parameters for the legend. The default is None.

members_plot_num (int, optional) – Number of individual members to plot. The default is 10.

members_alpha (float, optional) – Transparency of the lines representing the multiple members. The default is 0.3.

members_lw (float, optional) – With of the lines representing the multiple members. The default is 1.

seed (int, optional) – Set the seed for random number generator. Useful for reproducibility. The default is None.

Returns:

fig (matplotlib.pyplot.figure)

ax (matplotlib.pyplot.axis)

See also

pyleoclim.core.psds.PSD.plot
plotting method for PSD objects

pyleoclim.core.ensembleseries.EnsembleSeries.plot_envelope
envelope plot for ensembles

Examples
nn = 30 # number of noise realizations
nt = 500 # timeseries length
psds = []

time, signal = pyleo.utils.gen_ts(model='colored_noise',nt=nt,alpha=1.0)

ts = pyleo.Series(time=time, value = signal, verbose=False).standardize()
noise = np.random.randn(nt,nn)

for idx in range(nn):  # noise
    ts = pyleo.Series(time=time, value=signal+10*noise[:,idx], verbose=False)
    psd = ts.spectral()
    psds.append(psd)

mPSD = pyleo.MultiplePSD(psds)

fig, ax = mPSD.plot_envelope()
plot_traces(figsize=[10, 4], in_loglog=True, in_period=True, xlabel=None, ylabel='PSD', title=None, num_traces=10, seed=None, xticks=None, yticks=None, xlim=None, ylim=None, linestyle='-', savefig_settings=None, ax=None, legend=True, color='#d9544d', lw=0.5, alpha=0.3, lgd_kwargs=None)[source]
Plot MultiplePSD as a subset of traces.

Parameters:

figsize (list, optional) – The figure size. The default is [10, 4].

in_loglog (bool, optional) – Plot in log space. The default is True.

in_period (bool, optional) – Whether to plot periodicity instead of frequency. The default is True.

xlabel (str, optional) – x-axis label. The default is None.

ylabel (str, optional) – y-axis label. The default is None.

title (str, optional) – Plot title. The default is None.

xticks (list, optional) – x-ticks label. The default is None.

yticks (list, optional) – y-ticks label. The default is None.

xlim (list, optional) – x-axis limits. The default is None.

ylim (list, optional) – y-axis limits. The default is None.

color (str, optional) – Color of the traces. The default is sns.xkcd_rgb[‘pale red’].

alpha (float, optional) – Transparency of the lines representing the multiple members. The default is 0.3.

linestyle ({'-', '--', '-.', ':', '', (offset, on-off-seq), ...}) – Set the linestyle of the line

lw (float, optional) – Width of the lines representing the multiple members. The default is 0.5.

num_traces (int, optional) – Number of traces to plot. The default is 10.

savefig_settings (dict, optional) –
the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existed or non-existed path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”

”format” can be one of {“pdf”, “eps”, “png”, “ps”} The default is None.

ax (matplotlib.ax, optional) – Matplotlib axis on which to return the plot. The default is None.

legend (bool; {True,False}, optional) – Whether to plot the legend. The default is True.

lgd_kwargs (dict, optional) – Parameters for the legend. The default is None.

seed (int, optional) – Set the seed for the random number generator. Useful for reproducibility. The default is None.

Returns:

fig (matplotlib.figure) – the figure object from matplotlib See [matplotlib.pyplot.figure](https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.figure.html) for details.

ax (matplotlib.axis) – the axis object from matplotlib See [matplotlib.axes](https://matplotlib.org/api/axes_api.html) for details.

See also

pyleoclim.utils.plotting.savefig
Saving figure in Pyleoclim

Examples
nn = 30 # number of noise realizations
nt = 500 # timeseries length
psds = []

time, signal = pyleo.utils.gen_ts(model='colored_noise',nt=nt,alpha=1.0)

ts = pyleo.Series(time=time, value = signal, verbose=False).standardize()
noise = np.random.randn(nt,nn)

for idx in range(nn):  # noise
    ts = pyleo.Series(time=time, value=signal+10*noise[:,idx], verbose=False)
    psd = ts.spectral()
    psds.append(psd)

mPSD = pyleo.MultiplePSD(psds)

fig, ax = mPSD.plot_traces()
quantiles(qs=[0.05, 0.5, 0.95], lw=[0.5, 1.5, 0.5])[source]

Calculate the quantiles of the significance testing

Parameters:

qs (list, optional) – List of quantiles to consider for the calculation. The default is [0.05, 0.5, 0.95].

lw (list, optional) – Linewidth to use for plotting each level. Should be the same length as qs. The default is [0.5, 1.5, 0.5].

Raises:

ValueError – Frequency axis not consistent across the PSD list!

Returns:

psds

Return type:

pyleoclim.core.psds.MultiplePSD

Scalogram (pyleoclim.Scalogram)

class pyleoclim.core.scalograms.Scalogram(frequency, scale, time, amplitude, coi=None, label=None, Neff_threshold=3, wwz_Neffs=None, timeseries=None, wave_method=None, wave_args=None, signif_qs=None, signif_method=None, freq_method=None, freq_kwargs=None, scale_unit=None, time_label=None, signif_scals=None, qs=None)[source]

The Scalogram class is analogous to PSD, but for wavelet spectra (scalograms).

copy()[source]

Copy object

Returns:: scal – The copied version of the pyleoclim.Scalogram object
Return type:: pyleoclim.core.scalograms.Scalogram

Examples

series = pyleo.utils.load_dataset('NINO3')
scalogram = series.wavelet()
scalogram_copy = scalogram.copy()

plot(variable='amplitude', in_scale=True, xlabel=None, ylabel=None, title=None, ylim=None, xlim=None, yticks=None, figsize=[10, 8], signif_clr='white', signif_linestyles='-', signif_linewidths=1, contourf_style={}, cbar_style={}, plot_cb=True, savefig_settings={}, ax=None, signif_thresh=0.95)[source]

Plot the scalogram

Parameters:

in_scale (bool, optional) – Plot the in scale instead of frequency space. The default is True.
variable ({'amplitude','power'}) – Whether to plot the amplitude or power. Default is amplitude
xlabel (str, optional) – Label for the x-axis. The default is None.
ylabel (str, optional) – Label for the y-axis. The default is None.
title (str, optional) – Title for the figure. The default is ‘default’, which auto-generates a title.
ylim (list, optional) – Limits for the y-axis. The default is None.
xlim (list, optional) – Limits for the x-axis. The default is None.
yticks (list, optional) – yticks label. The default is None.
figsize (list, optional) – Figure size The default is [10, 8].
signif_clr (str, optional) – Color of the singificance line. The default is ‘white’.
signif_thresh (float in [0, 1]) – Significance threshold. Default is 0.95. If this quantile is not found in the qs field of the Coherence object, the closest quantile will be picked.
signif_linestyles (str, optional) – Linestyle of the significance line. The default is ‘-‘.
signif_linewidths (float, optional) – Width for the significance line. The default is 1.
contourf_style (dict, optional) – Arguments for the contour plot. The default is {}.
cbar_style (dict, optional) – Arguments for the colarbar. The default is {}.
savefig_settings (dict, optional) – saving options for the figure. The default is {}.
ax (ax, optional) – Matplotlib Axis on which to return the figure. The default is None.

Returns:

fig (matplotlib.figure) – the figure object from matplotlib See [matplotlib.pyplot.figure](https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.figure.html) for details.
ax (matplotlib.axis) – the axis object from matplotlib See [matplotlib.axes](https://matplotlib.org/api/axes_api.html) for details.

See also

pyleoclim.core.series.Series.wavelet: Wavelet analysis
pyleoclim.utils.plotting.savefig: Saving figure in Pyleoclim

Examples

ts = pyleo.utils.load_dataset('NINO3')
scalogram = ts.wavelet()
fig,ax = scalogram.plot()

signif_test(method='ar1sim', number=None, seed=None, qs=[0.95], settings=None, export_scal=False)[source]

Significance test for scalograms

Parameters:

method ({'ar1asym', 'ar1sim','uar1', 'CN'}) – Method to use to generate the surrogates. ar1sim uses simulated timeseries with similar persistence. ar1asym represents the theoretical, closed-form solution. The default is ar1sim
number (int) – Number of surrogates to generate for significance analysis based on simulations. The default is 200.
seed (int, optional) – Set the seed for the random number generator. Useful for reproducibility The default is None.
qs (list, optional) – Significane level to consider. The default is [0.95].
settings (dict, optional) – Parameters for the model. The default is None.
export_scal (bool; {True,False}) – Whether or not to export the scalograms used in the noise realizations. Note: For the wwz method, the scalograms used for wavelet analysis are slightly different than those used for spectral analysis (different decay constant). As such, this functionality should be used only to expedite exploratory analysis.

Raises:

ValueError – qs should be a list with at least one value.

Returns:

new – A new Scalogram object with the significance level

Return type:

pyleoclim.core.scalograms.Scalogram

See also

pyleoclim.core.series.Series.wavelet: Wavelet analysis
pyleoclim.core.scalograms.MultipleScalogram: MultipleScalogram object
pyleoclim.utils.wavelet.tc_wave_signif: Asymptotic significance calculation
pyleoclim.core.surrogateseries.SurrogateSeries: surrogate series objects

Examples

Generating scalogram, running significance tests, and saving the output for future use in generating psd objects or in summary_plot()

ts = pyleo.utils.load_dataset('NINO3')
scalogram = ts.wavelet().signif_test(number=2, export_scal=True)

By setting export_scal to True, the noise realizations used to generate the significance test will be saved. These come in handy for generating summary plots and for running significance tests on spectral objects.

MultipleScalogram (pyleoclim.MultipleScalogram)

class pyleoclim.core.scalograms.MultipleScalogram(scalogram_list)[source]

MultipleScalogram objects are used to store the results of significance testing for wavelet analysis

copy()[source]

Copy the object

See also

pyleoclim.core.scalograms.Scalogram.copy: Scalogram object copy

quantiles(qs=[0.05, 0.5, 0.95])[source]

Calculate quantiles

Parameters:

qs (list, optional) – List of quantiles to consider for the calculation. The default is [0.05, 0.5, 0.95].

Raises:

ValueError – Frequency axis not consistent across the PSD list!
Value Error – Time axis not consistent across the scalogram list!

Returns:

scals

Return type:

pyleoclim.core.scalograms.MultipleScalogram

Coherence (pyleoclim.Coherence)

class pyleoclim.core.coherences.Coherence(frequency, scale, time, wtc, xwt, phase, coi=None, wave_method=None, wave_args=None, timeseries1=None, timeseries2=None, signif_qs=None, signif_method=None, qs=None, freq_method=None, freq_kwargs=None, Neff_threshold=3, scale_unit=None, time_label=None)[source]

Coherence object, meant to receive the WTC and XWT part of Series.wavelet_coherence()

See also

pyleoclim.core.series.Series.wavelet_coherence: Wavelet coherence method

copy()[source]: Copy object

dashboard(title=None, figsize=[9, 12], overlap=True, phase_style={}, line_colors=['tab:blue', 'tab:orange'], savefig_settings={}, ts_plot_kwargs=None, wavelet_plot_kwargs=None)[source]

Cross-wavelet dashboard, including the two series, their WTC and XWT.

Note: this design balances many considerations, and is not easily customizable.

Parameters:

title (str, optional) – Title of the plot. The default is None.
figsize (list, optional) – Figure size. The default is [9, 12], as this is an information-rich figure.
overlap (boolean, optional) – whether to restrict the plot to the period of overlap between the series. Defaults to True
phase_style (dict, optional) – Arguments for the phase arrows. The default is {}. It includes: - ‘pt’: the default threshold above which phase arrows will be plotted - ‘skip_x’: the number of points to skip between phase arrows along the x-axis - ‘skip_y’: the number of points to skip between phase arrows along the y-axis - ‘scale’: number of data units per arrow length unit (see matplotlib.pyplot.quiver) - ‘width’: shaft width in arrow units (see matplotlib.pyplot.quiver) - ‘color’: arrow color (see matplotlib.pyplot.quiver)
line_colors (list, optional) – Colors for the 2 traces For nomenclature, see https://matplotlib.org/stable/gallery/color/named_colors.html
savefig_settings (dict, optional) –
The default is {}. the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existed or non-existed path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”}
ts_plot_kwargs (dict) – arguments to be passed to the timeseries subplot, see pyleoclim.core.series.Series.plot for details
wavelet_plot_kwargs (dict) – arguments to be passed to the contour subplots (XWT and WTC), [see pyleoclim.core.coherence.Coherence.plot for details]

Return type:

fig, ax

See also

pyleoclim.core.coherence.Coherence.plot: creates a coherence plot
pyleoclim.core.series.Series.wavelet_coherence: computes the coherence between two timeseries.
pyleoclim.core.series.Series.plot: plots a timeseries
matplotlib.pyplot.quiver: makes a quiver plot

Examples

Calculate the coherence of NINO3 and All India Rainfall and plot it as a dashboard:

ts_air = pyleo.utils.load_dataset('AIR')
ts_nino = pyleo.utils.load_dataset('NINO3')

coh = ts_air.wavelet_coherence(ts_nino)
coh_sig = coh.signif_test(number=10)

coh_sig.dashboard()

(<Figure size 900x1200 with 6 Axes>,
 {'ts1': <Axes: ylabel='AIR [mm/month]'>,
  'ts2': <Axes: xlabel='Time [year C.E.]', ylabel='NINO3 [$^{\\circ}$C]'>,
  'wtc': <Axes: ylabel='Scale [yrs]'>,
  'xwt': <Axes: xlabel='Time [year C.E.]', ylabel='Scale [yrs]'>})

You may customize colors like so:

coh_sig.dashboard(line_colors=['teal','gold'])

(<Figure size 900x1200 with 6 Axes>,
 {'ts1': <Axes: ylabel='AIR [mm/month]'>,
  'ts2': <Axes: xlabel='Time [year C.E.]', ylabel='NINO3 [$^{\\circ}$C]'>,
  'wtc': <Axes: ylabel='Scale [yrs]'>,
  'xwt': <Axes: xlabel='Time [year C.E.]', ylabel='Scale [yrs]'>})

To export the figure, use savefig_settings:

coh_sig.dashboard(savefig_settings={'path':'./coh_dash.png','dpi':300})

Figure saved at: "coh_dash.png"

(<Figure size 900x1200 with 6 Axes>,
 {'ts1': <Axes: ylabel='AIR [mm/month]'>,
  'ts2': <Axes: xlabel='Time [year C.E.]', ylabel='NINO3 [$^{\\circ}$C]'>,
  'wtc': <Axes: ylabel='Scale [yrs]'>,
  'xwt': <Axes: xlabel='Time [year C.E.]', ylabel='Scale [yrs]'>})

phase_stats(scales, number=1000, level=0.05)[source]

Estimate phase angle statistics of a Coherence object

As per [1], the strength (consistency) of a phase relationship is assessed using:

sigma, the circular standard deviation
kappa, an estimate of the Von Mises distribution’s concentration parameter.
It is a reciprocal measure of dispersion, so 1/kappa is analogous to the variance) [3].

Because of inherent persistence of geophysical signals and of the reproducing kernel of the continuous wavelet transform [3], phase statistics are assessed relative to an AR(1) model fit to the angle deviations observed at the requested scale(s).

Specifically, if number is specified, the method simulates number Monte Carlo realizations of an AR(1) process fit to fluctuations around the mean angle. This ensemble is used to obtain the confidence limits: sigma_lo (level quantile) and kappa_hi (1-level quantile). These correspond to 1-tailed tests of the strength of the relationship.

Parameters:

scales (float) – scale at which to evaluate the phase angle
number (int, optional) – number of AR(1) series to create for significance testing. The default is 1000.
level (float, optional) – significance level against which to gauge sigma and kappa. default: 0.05

Returns:

result – contains angle_mean (the mean angle for those scales), sigma (the circular standard deviation), kappa, sigma_lo (alpha-level quantile for sigma) and kappa_hi, the (1-alpha)-level quantile for kappa.

Return type:

dict

See also

pyleoclim.core.series.Series.wavelet_coherence: Wavelet coherence
pyleoclim.core.scalograms.Scalogram: Scalogram object
pyleoclim.core.scalograms.MultipleScalogram: Multiple Scalogram object
pyleoclim.core.coherence.Coherence.plot: plotting method for Coherence objects
pyleoclime.utils.wavelet.angle_sig: significance of phase angle statistics
pyleoclim.utils.wavelet.angle_stats: phase angle statistics

References

[1] Grinsted, A., J. C. Moore, and S. Jevrejeva (2004), Application of the cross wavelet transform and wavelet coherence to geophysical time series, Nonlinear Processes in Geophysics, 11, 561–566.

[2] Huber, R., Dutra, L. V., & da Costa Freitas, C. (2001). SAR interferogram phase filtering based on the Von Mises distribution. In IGARSS 2001. Scanning the Present and Resolving the Future. Proceedings. IEEE 2001 International Geoscience and Remote Sensing Symposium (Cat. No. 01CH37217) (Vol. 6, pp. 2816-2818). IEEE.

[3] Farge, M. and Schneider, K. (2006): Wavelets: application to turbulence Encyclopedia of Mathematical Physics (Eds. J.-P. Françoise, G. Naber and T.S. Tsun) pp 408-420.

Examples

Calculate the phase angle between NINO3 and All India Rainfall at 5y scales:

ts_air = pyleo.utils.load_dataset('AIR')
ts_nino = pyleo.utils.load_dataset('NINO3')
coh = ts_air.wavelet_coherence(ts_nino)
coh.phase_stats(scales=5)

Results(mean_angle=np.float64(-2.6915968149088703), kappa=np.float64(3.479951408674022), sigma=np.float64(0.5877766646462549), kappa_hi=np.float64(0.6849797644327951), sigma_lo=np.float64(1.5013474577057688))

One may also obtain phase angle statistics over an interval, like the 2-8y ENSO band:

phase = coh.phase_stats(scales=[2,8])
print("The mean angle is {:4.2f}°".format(phase.mean_angle/np.pi*180))
print(phase)

The mean angle is -154.29°
Results(mean_angle=np.float64(-2.6927825175802447), kappa=np.float64(3.360315074427689), sigma=np.float64(0.6012824554612396), kappa_hi=np.float64(0.4630633245129424), sigma_lo=np.float64(1.7258092470407183))

From this example, one diagnoses a strong anti-phased relationship in the ENSO band, with high von Mises concentration (kappa ~ 3.35 >> kappa_hi) and low circular dispersion (sigma ~ 0.6 << sigma_lo). This would be strong evidence of a consistent anti-phasing between NINO3 and AIR at those scales.

plot(var='wtc', xlabel=None, ylabel=None, title='auto', figsize=[10, 8], ylim=None, xlim=None, in_scale=True, yticks=None, contourf_style={}, phase_style={}, cbar_style={}, savefig_settings={}, ax=None, signif_clr='white', signif_linestyles='-', signif_linewidths=1, signif_thresh=0.95, under_clr='ivory', over_clr='black', bad_clr='dimgray')[source]

Plot the cross-wavelet results

Parameters:

var (str {'wtc', 'xwt'}) – variable to be plotted as color field. Default: ‘wtc’, the wavelet transform coherency. ‘xwt’ plots the cross-wavelet transform instead.
xlabel (str, optional) – x-axis label. The default is None.
ylabel (str, optional) – y-axis label. The default is None.
title (str, optional) – Title of the plot. The default is ‘auto’, where it is made from object metadata. To mute, pass title = None.
figsize (list, optional) – Figure size. The default is [10, 8].
ylim (list, optional) – y-axis limits. The default is None.
xlim (list, optional) – x-axis limits. The default is None.
in_scale (bool, optional) – Plots scales instead of frequencies The default is True.
yticks (list, optional) – y-ticks label. The default is None.
contourf_style (dict, optional) – Arguments for the contour plot. The default is {}.
phase_style (dict, optional) – Arguments for the phase arrows. The default is {}. It includes: - ‘pt’: the default threshold above which phase arrows will be plotted - ‘skip_x’: the number of points to skip between phase arrows along the x-axis - ‘skip_y’: the number of points to skip between phase arrows along the y-axis - ‘scale’: number of data units per arrow length unit (see matplotlib.pyplot.quiver) - ‘width’: shaft width in arrow units (see matplotlib.pyplot.quiver) - ‘color’: arrow color (see matplotlib.pyplot.quiver)
cbar_style (dict, optional) – Arguments for the color bar. The default is {}.
savefig_settings (dict, optional) –
The default is {}. the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existed or non-existed path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”}
ax (ax, optional) – Matplotlib axis on which to return the figure. The default is None.
signif_thresh (float in [0, 1]) – Significance threshold. Default is 0.95. If this quantile is not found in the qs field of the Coherence object, the closest quantile will be picked.
signif_clr (str, optional) – Color of the significance line. The default is ‘white’.
signif_linestyles (str, optional) – Style of the significance line. The default is ‘-‘.
signif_linewidths (float, optional) – Width of the significance line. The default is 1.
under_clr (str, optional) – Color for under 0. The default is ‘ivory’.
over_clr (str, optional) – Color for over 1. The default is ‘black’.
bad_clr (str, optional) – Color for missing values. The default is ‘dimgray’.

Return type:

fig, ax

See also

pyleoclim.core.coherence.Coherence.dashboard: plots a a dashboard showing the coherence and the cross-wavelet transform.
pyleoclim.core.series.Series.wavelet_coherence: computes the coherence from two timeseries.
matplotlib.pyplot.quiver: quiver plot

Examples

Calculate the wavelet coherence of NINO3 and All India Rainfall and plot it: .. jupyter-execute:

ts_air = pyleo.utils.load_dataset('AIR')
ts_nino = pyleo.utils.load_dataset('NINO3')

coh = ts_air.wavelet_coherence(ts_nino)
coh.plot()

Establish significance against an AR(1) benchmark:

coh_sig = coh.signif_test(number=20, qs=[.9,.95,.99])
coh_sig.plot()

(<Figure size 1000x800 with 2 Axes>,
 <Axes: title={'center': 'Lines:95% threshold'}, xlabel='Time [year C.E.]', ylabel='Scale [yrs]'>)

Note that specifiying 3 significance thresholds does not take any more time as the quantiles are simply estimated from the same ensemble. By default, the plot function looks for the closest quantile to 0.95, but this is easy to adjust, e.g. for the 99th percentile:

coh_sig.plot(signif_thresh = 0.99)

(<Figure size 1000x800 with 2 Axes>,
 <Axes: title={'center': 'Lines:99% threshold'}, xlabel='Time [year C.E.]', ylabel='Scale [yrs]'>)

By default, the function plots the wavelet transform coherency (WTC), which quantifies where two timeseries exhibit similar behavior in time-frequency space, regardless of whether this corresponds to regions of high common power. To visualize the latter, you want to plot the cross-wavelet transform (XWT) instead, like so:

coh_sig.plot(var='xwt')

(<Figure size 1000x800 with 2 Axes>,
 <Axes: title={'center': 'Lines:95% threshold'}, xlabel='Time [year C.E.]', ylabel='Scale [yrs]'>)

signif_test(number=200, method='ar1sim', seed=None, qs=[0.95], settings=None, mute_pbar=False)[source]

Significance testing for Coherence objects

The method obtains quantiles qs of the distribution of coherence between number pairs of Monte Carlo simulations of a process that resembles the original series. Currently, only AR(1) surrogates are supported.

Parameters:

number (int, optional) – Number of surrogate series to create for significance testing. The default is 200.
method ({'ar1sim','phaseran','CN'}, optional) – Method through which to generate the surrogate series. The default is ‘phaseran’.
seed (int, optional) – Fixes the seed for NumPy’s random number generator. Useful for reproducibility. The default is None, so fresh, unpredictable entropy will be pulled from the operating system.
qs (list, optional) – Significance levels to return. The default is [0.95].
settings (dict, optional) – Parameters for surrogate model. The default is None.
mute_pbar (bool, optional) – Mute the progress bar. The default is False.

Returns:

new – original Coherence object augmented with significance levels signif_qs, a list with the following MultipleScalogram objects: * 0: MultipleScalogram for the wavelet transform coherency (WTC) * 1: MultipleScalogram for the cross-wavelet transform (XWT)

Each object contains as many Scalogram objects as qs contains values

Return type:

pyleoclim.core.coherence.Coherence

See also

pyleoclim.core.series.Series.wavelet_coherence: Wavelet coherence
pyleoclim.core.scalograms.Scalogram: Scalogram object
pyleoclim.core.scalograms.MultipleScalogram: Multiple Scalogram object
pyleoclim.core.coherence.Coherence.plot: plotting method for Coherence objects

Examples

Calculate the coherence of NINO3 and All India Rainfall and assess significance:

ts_air = pyleo.utils.load_dataset('AIR')
ts_nino = pyleo.utils.load_dataset('NINO3')

coh = ts_air.wavelet_coherence(ts_nino)
coh_sig = coh.signif_test(number=20)
coh_sig.plot()

(<Figure size 1000x800 with 2 Axes>,
 <Axes: title={'center': 'Lines:95% threshold'}, xlabel='Time [year C.E.]', ylabel='Scale [yrs]'>)

By default, significance is assessed against a 95% benchmark derived from an AR(1) process fit to the data, using 200 Monte Carlo simulations. To customize, one can increase the number of simulations (more reliable, but slower), and the quantile levels.

coh_sig2 = coh.signif_test(number=100, qs=[.9,.95,.99])
coh_sig2.plot()

(<Figure size 1000x800 with 2 Axes>,
 <Axes: title={'center': 'Lines:95% threshold'}, xlabel='Time [year C.E.]', ylabel='Scale [yrs]'>)

The plot() function will represent the 95% level as contours by default. If you need to show 99%, say, use the signif_thresh argument:

coh_sig2.plot(signif_thresh=0.99)

(<Figure size 1000x800 with 2 Axes>,
 <Axes: title={'center': 'Lines:99% threshold'}, xlabel='Time [year C.E.]', ylabel='Scale [yrs]'>)

Note that if the 99% quantile is not present, the plot method will look for the closest match, but lines are always labeled appropriately. For reproducibility purposes, it may be good to specify the (pseudo)random number generator’s seed, like so:

coh_sig27 = coh.signif_test(number=20, seed=27)

This will generate exactly the same set of draws from the (pseudo)random number at every execution, which may be important for marginal features in small ensembles. In general, however, we recommend increasing the number of draws to check that features are robust.

One can also specifiy a different method to obtain surrogates, e.g. phase randomization:

coh.signif_test(method='phaseran').plot()

(<Figure size 1000x800 with 2 Axes>,
 <Axes: title={'center': 'Lines:95% threshold'}, xlabel='Time [year C.E.]', ylabel='Scale [yrs]'>)

GlobalCoherence (pyleoclim.GlobalCoherence)

class pyleoclim.core.coherences.GlobalCoherence(global_coh, coh, signif_qs=None, signif_method=None, qs=None, label='Coherence')[source]

Class to store the results of cross spectral analysis

Parameters:

global_coh (numpy array) – coherence values
scale (numpy array) – scale values
frequency (numpy array) – frequency values
coi (numpy array) – cone of influence values
coh (Coherence) – Original coherence object

See also

pyleoclim.core.series.Series.global_coherence: method to compute the spectral coherence

copy()[source]: Copy object

plot(figsize=(8, 8), xlim=None, xlabel=None, label=None, psd_y_label='PSD', coh_y_label='Coherence', coh_line_color='grey', ax=None, coh_ylim=(0.4, 1), fill_alpha=0.3, fill_color='grey', coh_plot_kwargs=None, savefig_settings=None, spectral_kwargs=None, legend=True, legend_kwargs=None, spec1_plot_kwargs=None, spec2_plot_kwargs=None)[source]

Plot the coherence as a function of scale or frequency, alongside the spectrum of the two timeseries (using the same method used for the coherence).

Parameters:

figsize (tuple) – size of the figure. Default is (8,8). Only used if ax is None
xlim (tuple) – x limits for the plot. Default is None
label (str) – label of the plot
xlabel (str) – x label of the plot
psd_y_label (str) – y label of the power spectral density plot (left hand side)
coh_y_label (str) – y label of the coherence plot (right hand side)
coh_line_color (str) – color of the coherence line
coh_ylim (tuple) – y limits for the coherence plot. Default is (.4,1)
fill_alpha (float) – alpha value for the fill_between plot. Default is .3
fill_color (str) – color of the fill_between plot
coh_plot_kwargs (dict) – additional arguments to pass to the pyleoclim.utils.plotting.plot_xy
savefig_settings (dict) – settings to pass to the pyleoclim.utils.plotting.savefig function
spectral_kwargs (dict) – additional arguments to pass to the pyleo.Series.spectral method
spec1_plot_kwargs (dict) – additional arguments to pass to the pyleo.Series.spectral method
spec2_plot_kwargs (dict) – additional arguments to pass to the pyleo.Series.spectral method
legend (bool) – whether to include a legend or not
legend_kwargs (dict) – additional arguments to pass to ax.legend
ax (matplotlib axis) – axis to plot on

Returns:

ax – axis with the plot

Return type:

matplotlib axis

Examples

soi = pyleo.utils.load_dataset('SOI')
nino3 = pyleo.utils.load_dataset('NINO3')

gcoh = soi.global_coherence(nino3)
gcoh.plot()

(<Figure size 800x800 with 2 Axes>,
 <Axes: xlabel='Period [year]', ylabel='PSD'>)

signif_test(method='ar1sim', number=200, qs=[0.95])[source]

Perform a significance test on the coherence values

Parameters:

method (str; {'ar1sim','CN','phaseran'}) – method to use for the surrogate test. Default is ‘ar1sim’.
number (int) – number of surrogates to generate. Default is 200
qs (list) – list of quantiles to compute. Default is [.95]

Returns:

global_coh – Global coherence with significance field filled in

Return type:

pyleoclim.core.globalcoherence.GlobalCoherence

Examples

soi = pyleo.utils.load_dataset('SOI')
nino3 = pyleo.utils.load_dataset('NINO3')

gcoh = soi.global_coherence(nino3)
gcoh_sig = gcoh.signif_test(number=10)
gcoh_sig.plot()

(<Figure size 800x800 with 2 Axes>,
 <Axes: xlabel='Period [year]', ylabel='PSD'>)

Corr (pyleoclim.Corr)

class pyleoclim.core.corr.Corr(r, p, r_crit, signif, alpha, p_fmt_td=0.01, p_fmt_style='exp')[source]

The object for correlation results in order to format the print message

Parameters:

r (float) – the metric of association (see https://docs.scipy.org/doc/scipy/reference/stats.html#association-correlation-tests)
p (float) – the p-value
r_crit (float) – critical value for significance
p_fmt_td (float) – the threshold for p-value formatting (0.01 by default, i.e., if p<0.01, will print “< 0.01” instead of “0”)
p_fmt_style (str) – the style for p-value formatting (exponential notation by default)
signif (bool) – the significance
alpha (float) – The significance level (0.05 by default)

See also

pyleoclim.utils.correlation.fdr: FDR function
pyleoclim.utils.correlation.association: workhorse function to compute various metrics of association

copy()[source]: Copy object

CorrEns (pyleoclim.CorrEns)

class pyleoclim.core.correns.CorrEns(r, p, signif, signif_fdr, alpha, p_fmt_td=0.01, p_fmt_style='exp')[source]

CorrEns objects store the result of an ensemble correlation calculation between timeseries and/or ensemble of timeseries. The class enables a print and plot function to easily visualize the result.

Parameters:

r (list) – the list of correlation coefficients
p (list) – the list of p-values
p_fmt_td (float) – the threshold for p-value formating (0.01 by default, i.e., if p<0.01, will print “< 0.01” instead of “0”)
p_fmt_style (str) – the style for p-value formating (exponential notation by default)
signif (list) – the list of significance without FDR
signif_fdr (list) – the list of significance with FDR
signif_fdr – the list of significance with FDR
alpha (float) – The significance level

See also

pyleoclim.utils.correlation.association: workhorse function to compute various metrics of association
pyleoclim.utils.correlation.fdr: FDR (False Discovery Rate) function

copy()[source]: Copy object

plot(figsize=[4, 4], title=None, ax=None, savefig_settings=None, hist_kwargs=None, title_kwargs=None, xlim=None, alpha=0.8, multiple='layer', clr_insignif='silver', clr_signif='#029386', clr_signif_fdr='darkorange', clr_percentile='#ff796c')[source]

Plot the distribution of correlation values as a histogram

Uses seaborn’s histplot

Color-coding is used to indicate significance, with or without applying the False Discovery Rate (FDR) method.

Parameters:

figsize (list, optional) – The figure size. The default is [4, 4].
title (str, optional) – Plot title. The default is None.
multiple (str, optional) – Approach to organizing the 3 different histrograms on the plot. possible values: “layer”[default], “dodge”, “stack”, “fill”
alpha (float in [0, 1]) – transparency parameter for histrogram bars. Default: 0.8
savefig_settings (dict) –
the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existing or new path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”}
hist_kwargs (dict) – additional keyword arguments for sns.histplot() [experimental]
title_kwargs (dict) – the keyword arguments for ax.set_title()
ax (matplotlib.axis, optional) – the axis object from matplotlib See [matplotlib.axes](https://matplotlib.org/api/axes_api.html) for details.
xlim (list, optional) – x-axis limits. The default is None.

See also

pyleoclim.core.series.Series.correlation: correlation with significance
pyleoclim.utils.correlation.fdr: False Discovery Rate
seaborn.histplot: https://seaborn.pydata.org/generated/seaborn.histplot.html
pyleoclim.utils.plotting.savefig: save figures in Pyleoclim

MultivarDecomp (pyleoclim.MultivariateDecomp)

class pyleoclim.core.multivardecomp.MultivariateDecomp(name, eigvals, eigvecs, pctvar, pcs, neff, orig)[source]

Class to hold the results of multivariate decompositions applies to : pca(), mcpca(), mssa()

Parameters:

time (float) – the common time axis
name (str) – name of the dataset/analysis to use in plots
eigvals (1d array) – vector of eigenvalues from the decomposition
eigvecs (2d array) – array of eigenvectors from the decomposition (e.g. EOFs)
pcs (1d array) – array containing the temporal expansion coefficients (e.g. “principal components” in the climate lore)
pctvar (float) – array of pct variance accounted for by each mode
orig (MultipleSeries, or MultipleGeoSeries object) – original data, on a common time axis
neff (float) – scalar representing the effective sample size of the leading mode

modeplot(index=0, figsize=[8, 8], fig=None, savefig_settings=None, gs=None, title=None, title_kwargs=None, spec_method='mtm', cmap=None, hue='EOF', marker=None, size=None, scatter_kwargs=None, flip=False, map_kwargs=None, gridspec_kwargs=None)[source]

Dashboard visualizing the properties of a given mode.

Includes: The temporal coefficient (PC or similar), its spectrum, and the loadings (EOF or similar), possibly geolocated. If the object does not have geolocation information, a spaghetti plot of the standardized series is displayed.

To see how to use the different parameters, consult the following tutorial: http://linked.earth/PyleoTutorials/notebooks/L2_principal_component_analysis.html

Parameters:

indexint

the (0-based) index of the mode to visualize. Default is 0, corresponding to the first mode.

figsizelist, optional

The figure size. The default is [8, 8].

savefig_settingsdict

the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existed or non-existed path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”

“format” can be one of {“pdf”, “eps”, “png”, “ps”}

titlestr, optional

text for figure title

title_kwargsdict

the keyword arguments for ax.set_title()

gsmatplotlib.gridspec object, optional

Requires at least two rows and two columns. - top row, left: timeseries of principle component - top row, right: PSD - bottom row: spaghetti plot or map See [matplotlib.gridspec.GridSpec](https://matplotlib.org/stable/tutorials/intermediate/gridspec.html) for details.

gridspec_kwargsdict, optional

Dictionary with custom gridspec values. - wspace changes space between columns (default: wspace=0.05) - hspace changes space between rows (default: hspace=0.03) - width_ratios: relative width of each column (default: width_ratios=[5,1,3] where middle column serves as a spacer) - height_ratios: relative height of each row (default: height_ratios=[2,1,5] where middle row serves as a spacer)

spec_method: str, optional

The name of the spectral method to be applied on the PC. Default: MTM Note that the data are evenly-spaced, so any spectral method that assumes even spacing is applicable here: ‘mtm’, ‘welch’, ‘periodogram’ ‘wwz’ is relevant if scaling exponents need to be estimated, but ill-advised otherwise, as it is very slow.

cmap: str, optional

if ‘hue’ is specified, will be used for map scatter plot values. colormap name for the loadings (https://matplotlib.org/stable/tutorials/colors/colormaps.html)

map_kwargsdict, optional

Optional arguments for map configuration - projection: str; Optional value for map projection. Default ‘auto’. - proj_default: bool - lakes, land, ocean, rivers, borders, coastline, background: bool or dict; - lgd_kwargs: dict; Optional values for how the map legend is configured - gridspec_kwargs: dict; Optional values for adjusting the arrangement of the colorbar, map and legend in the map subplot - legend: bool; Whether to draw a legend on the figure. Default is True - colorbar: bool; Whether to draw a colorbar on the figure if the data associated with hue are numeric. Default is True - color_scale_type : str; Setting to “discrete” will force a discrete color scale with a default bin number of max(11, n) where n=number of unique values$^{

rac{1}{2}}$. Default is None

scalar_mappable: matplotlib.cm.ScalarMappable; can be used to pass a matplotlib scalar mappable. See pyleoclim.utils.plotting.make_scalar_mappable for documentation on using the Pyleoclim utility, or the Matplotlib tutorial on customizing colorbars.

scatter_kwargsdict, optional: Optional arguments configuring how data are plotted on a map. This allows you to configure information about the style of markers (e.g., triangle, squares, circles), their size (through markersize), edgecolor, etc…) See description of scatter_kwargs in pyleoclim.utils.mapping.scatter_map
huestr, optional: (only applicable if using scatter map) Variable associated with color coding for points plotted on map. May correspond to a continuous or categorical variable. The default is ‘EOF’.
sizestr, optional: (only applicable if using scatter map) Variable associated with size. Must correspond to a continuous numeric variable. The default is None.
markerstring, optional: (only applicable if using scatter map) Grouping variable that will produce points with different markers. Can have a numeric dtype but will always be treated as categorical. The default is None, which will produce circle markers. Alternatively, pass the name of a categorical variable, e.g. ‘archiveType’. If ‘archiveType’ is specified, will attempt to use pyleoclim archiveType markers mapping, defaulting to ‘?’ where values are unavailable. Caution: this is different from the marker argument in matplotlib. If you want to specify a marker style (e.g., cirle, square, triangle, use scatter_kwargs)

Returns:

figmatplotlib.figure: The figure
axdict: dictionary of matplotlib ax

See also

pyleoclim.core.MultipleSeries.pca: Principal Component Analysis
pyleoclim.core.MultipleGeoSeries.pca: Principal Component Analysis
pyleoclim.utils.tsutils.eff_sample_size: Effective sample size
pyleoclim.utils.mapping.scatter_map: mapping
pyleoclim.utils.plotting.make_scalar_mappable: Custom scalar mappable

screeplot(figsize=[6, 4], uq='N82', title=None, ax=None, savefig_settings=None, title_kwargs=None, xlim=[0, 10], clr_eig='C0')[source]

Plot the eigenvalue spectrum with uncertainties

Parameters:

figsize (list, optional) – The figure size. The default is [6, 4].
title (str, optional) – Plot title. The default is ‘scree plot’.
savefig_settings (dict) –
the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existed or non-existed path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”}
title_kwargs (dict, optional) – the keyword arguments for ax.set_title()
ax (matplotlib.axis, optional) – the axis object from matplotlib See [matplotlib.axes](https://matplotlib.org/api/axes_api.html) for details.
xlim (list, optional) – x-axis limits. The default is [0, 10] (first 10 eigenvalues)
uq (str, optional) – Method used for uncertainty quantification of the eigenvalues. ‘N82’ uses the North et al “rule of thumb” [1] with effective sample size computed as in [2]. ‘MC’ uses Monte-Carlo simulations (e.g. MC-EOF). Returns an error if no ensemble is found.
clr_eig (str, optional) – color to be used for plotting eigenvalues

See also

pyleoclim.core.MultipleSeries.pca: Principal Component Analysis

References

[1]_ North, G. R., T. L. Bell, R. F. Cahalan, and F. J. Moeng (1982), Sampling errors in the estimation of empirical orthogonal functions, Mon. Weather Rev., 110, 699–706.

[2]_ Hannachi, A., I. T. Jolliffe, and D. B. Stephenson (2007), Empirical orthogonal functions and related techniques in atmospheric science: A review, International Journal of Climatology, 27(9), 1119–1152, doi:10.1002/joc.1499.

EnsMultivarDecomp (pyleoclim.EnsMultivarDecomp)

class pyleoclim.core.ensmultivardecomp.EnsMultivarDecomp(pca_list, label=None)[source]

modeplot(index=0, flip=False, plot_envelope_kwargs=None, psd_envelope_kwargs=None, figsize=[8, 8], fig=None, savefig_settings=None, gs=None, title=None, title_kwargs=None, spec_method='mtm', cmap='coolwarm', marker=None, scatter_kwargs=None, map_kwargs=None, gridspec_kwargs=None, quantiles=[0.25, 0.5, 0.75])[source]

Plot relevant information about the specific mode

Parameters:

index (int) – The 1-based index of the mode to visualize. Default is 1, corresponding to the first mode.
flip (bool) – Whether or not to flip the PC
plot_envelope_kwargs (dict) – Dictionary of key word arguments for plot envelope
plot_envelope_kwargs – Dictionary of key word arguments for psd envelope
figsize (list, optional) – The figure size. The default is [8, 8].
savefig_settings (dict) –
the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existed or non-existed path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”}
title (str, optional) – text for figure title
title_kwargs (dict) – the keyword arguments for ax.set_title()
gs (matplotlib.gridspec object, optional) – Requires at least two rows and two columns. - top row, left: timeseries of principle component - top row, right: PSD - bottom row: spaghetti plot or map See [matplotlib.gridspec.GridSpec](https://matplotlib.org/stable/tutorials/intermediate/gridspec.html) for details.
gridspec_kwargs (dict, optional) – Dictionary with custom gridspec values. - wspace changes space between columns (default: wspace=0.05) - hspace changes space between rows (default: hspace=0.03) - width_ratios: relative width of each column (default: width_ratios=[5,1,3] where middle column serves as a spacer) - height_ratios: relative height of each row (default: height_ratios=[2,1,5] where middle row serves as a spacer)
spec_method (str, optional) – The name of the spectral method to be applied on the PC. Default: MTM Note that the data are evenly-spaced, so any spectral method that assumes even spacing is applicable here: ‘mtm’, ‘welch’, ‘periodogram’ ‘wwz’ is relevant if scaling exponents need to be estimated, but ill-advised otherwise, as it is very slow.
cmap (str, optional) – if ‘hue’ is specified, will be used for map scatter plot values. colormap name for the loadings (https://matplotlib.org/stable/tutorials/colors/colormaps.html)
map_kwargs (dict, optional) – Optional arguments for map configuration - projection: str; Optional value for map projection. Default ‘auto’. - proj_default: bool - lakes, land, ocean, rivers, borders, coastline, background: bool or dict; - lgd_kwargs: dict; Optional values for how the map legend is configured - gridspec_kwargs: dict; Optional values for adjusting the arrangement of the colorbar, map and legend in the map subplot - legend: bool; Whether to draw a legend on the figure. Default is True - colorbar: bool; Whether to draw a colorbar on the figure if the data associated with hue are numeric. Default is True The default is None.
scatter_kwargs (dict, optional) – Optional arguments configuring how data are plotted on a map. See description of scatter_kwargs in pyleoclim.utils.mapping.scatter_map
hue (str, optional) – (only applicable if using scatter map) Variable associated with color coding for points plotted on map. May correspond to a continuous or categorical variable. The default is ‘EOF’.
marker (string, optional) – (only applicable if using scatter map) Grouping variable that will produce points with different markers. Can have a numeric dtype but will always be treated as categorical. The default is ‘archiveType’.
quantiles (list,array) – Quantiles to use for plotting EOFs

Returns:

fig (matplotlib.figure) – The figure
ax (dict) – dictionary of matplotlib ax

See also

pyleoclim.core.MultipleSeries.pca: Principal Component Analysis
pyleoclim.core.MultipleGeoSeries.pca: Principal Component Analysis
pyleoclim.utils.tsutils.eff_sample_size: Effective sample size
pyleoclim.utils.mapping.scatter_map: mapping

Examples

n = 100 # number of series

soi = pyleo.utils.load_dataset('SOI')
soi_time_axes = [pyleo.utils.random_time_axis(n=len(soi.time)) for _ in range(n)]
soi_ens = pyleo.EnsembleGeoSeries([pyleo.GeoSeries(time=time, value=soi.value,lat=-5,lon=-85,auto_time_params=True,verbose=False) for time in soi_time_axes])

nino3 = pyleo.utils.load_dataset('NINO3')
nino3_time_axes = [pyleo.utils.random_time_axis(n=len(nino3.time)) for _ in range(n)]
nino3_ens = pyleo.EnsembleGeoSeries([pyleo.GeoSeries(time=time, value=nino3.value,lat=-5,lon=-85,auto_time_params=True,verbose=False) for time in nino3_time_axes])

mul_ens = pyleo.MulEnsGeoSeries([nino3_ens,soi_ens])
mcpca = mul_ens.mcpca(nsim=10,seed=42)
mcpca.modeplot()

(<Figure size 800x800 with 4 Axes>,
 {'pc': <Axes: xlabel='Time [years CE]', ylabel='$PC_1$'>,
  'psd': <Axes: xlabel='Period [years CE]', ylabel='PSD'>,
  'map': {'map': <GeoAxes: xlabel='lon', ylabel='lat'>, 'leg': <Axes: >}})

screeplot(clr_eig='C0', linewidth=0.3, title='Screeplot', violin_kwargs=None, figsize=[8, 8], savefig_settings=None, ax=None)[source]

Function to plot the scree plot of the PCA

Parameters:

quantiles (list,array) – Quantiles to use for plotting range of pctvar for each mode.
clr_eig (str) – Color to use for the eigenvalues. Default is ‘C0’.
linewidth (float) – Linewidth to use for the violin plot. Default is 0.3.
violin_kwargs (dict) – Dictionary of key word arguments for violin plot. If exposed plot arguments {‘color’,’linewidth’,’title’} are included in kwargs, kwargs will overwrite the exposed argument values.
figsize (list, optional) – The figure size. The default is [8, 8].
savefig_settings (dict) –
the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existed or non-existed path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”}
ax (matplotlib.axis, optional) – the axis object from matplotlib See [matplotlib.axes](https://matplotlib.org/api/axes_api.html) for details.

Returns:

fig (matplotlib.figure) – The figure
ax (dict) – dictionary of matplotlib ax

Examples

n = 100 # number of series

soi = pyleo.utils.load_dataset('SOI')
soi_time_axes = [pyleo.utils.random_time_axis(n=len(soi.time)) for _ in range(n)]
soi_ens = pyleo.EnsembleGeoSeries([pyleo.GeoSeries(time=time, value=soi.value,lat=0,lon=0,auto_time_params=True,verbose=False) for time in soi_time_axes])

nino3 = pyleo.utils.load_dataset('NINO3')
nino3_time_axes = [pyleo.utils.random_time_axis(n=len(nino3.time)) for _ in range(n)]
nino3_ens = pyleo.EnsembleGeoSeries([pyleo.GeoSeries(time=time, value=nino3.value,lat=0,lon=0,auto_time_params=True,verbose=False) for time in nino3_time_axes])

mul_ens = pyleo.MulEnsGeoSeries([nino3_ens,soi_ens])
mcpca = mul_ens.mcpca(nsim=10,seed=42)
mcpca.screeplot()

(<Figure size 800x800 with 1 Axes>,
 <Axes: title={'center': 'Screeplot'}, xlabel='Mode index $i$', ylabel='$\\lambda_i$'>)

SsaRes (pyleoclim.SsaRes)

class pyleoclim.core.ssares.SsaRes(orig, label, eigvals, eigvecs, pctvar, PC, RCmat, RCseries, mode_idx, eigvals_q=None)[source]

This class is meant to hold the output of the Singular Spectrum Analysis (SSA) method, which applies to Series objets. Two functions are enabled by this class:

screeplot, which plots the eigenvalue spectrum to help determine what modes to keep
modeplot, which plots the individual mode temporal EOF and temporal PC

Parameters:

orig (Series) – timeseries on which SSA was performed
eigvals (float (M, 1)) – a vector of real eigenvalues derived from the signal
pctvar (float (M, 1)) – same vector, expressed in % variance accounted for by each mode.
eigvals_q (float (M, 2)) – array containing the 5% and 95% quantiles of the Monte-Carlo eigenvalue spectrum [ assigned NaNs if unused ]
eigvecs (float (M, M)) – a matrix of the temporal eigenvectors (T-EOFs), i.e. the temporal patterns that explain most of the variations in the original series.
PC (float (N - M + 1, M)) – array of principal components, i.e. the loadings that, convolved with the T-EOFs, produce the reconstructed components, or RCs
RCmat (float (N, M)) – array of reconstructed components, One can think of each RC as the contribution of each mode to the timeseries, weighted by their eigenvalue (loosely speaking, their “amplitude”). Summing over all columns of RC recovers the original series. (synthesis, the reciprocal operation of analysis).
mode_idx (list) – index of retained modes
RCseries (float (N, 1)) – reconstructed series based on the RCs of mode_idx (scaled to original series; mean must be added after the fact)

See also

pyleoclim.utils.decomposition.ssa: Singular Spectrum Analysis

copy()[source]

Make a copy of the SsaRes object

Returns:: SsaRes – A copy of the SsaRes object
Return type:: pyleoclim.SsaRes

modeplot(index=0, figsize=[10, 5], savefig_settings=None, title_kwargs=None, spec_method='mtm', plot_original=True)[source]

Dashboard visualizing the properties of a given SSA mode, including:

the analyzing function (T-EOF)
the reconstructed component (RC)
its spectrum

Parameters:

index (int) – the (0-based) index of the mode to visualize. Default is 0, corresponding to the first mode.
figsize (list, optional) – The figure size. The default is [10, 5].
savefig_settings (dict) –
the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existed or non-existed path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”}
title_kwargs (dict) – the keyword arguments for ax.set_title()
spec_method (str, optional) – The name of the spectral method to be applied on the PC. Default: MTM Note that the data are evenly-spaced, so any spectral method that assumes even spacing is applicable here: ‘mtm’, ‘welch’, ‘periodogram’. ‘wwz’ is relevant too if scaling exponents need to be estimated.

See also

pyleoclim.core.series.Series.ssa: Singular Spectrum Analysis for timeseries objects
pyleoclim.utils.decomposition.ssa: Singular Spectrum Analysis utility
pyleoclim.core.ssares.SsaRes.screeplot: plot SSA eigenvalue spectrum

Examples

Plot the first SSA mode of the NINO3 index:

ts  = pyleo.utils.load_dataset('NINO3')
ssa = ts.ssa()

fig, ax = ssa.modeplot()

Plot the second mode (note 0-based indexing):

fig, ax = ssa.modeplot(index=1)

screeplot(figsize=[6, 4], title='SSA scree plot', ax=None, savefig_settings=None, title_kwargs=None, xlim=None, clr_mcssa='#e50000', clr_signif='#029386', clr_eig='black')[source]

Scree plot for SSA, visualizing the eigenvalue spectrum and indicating which modes were retained.

Parameters:

figsize (list, optional) – The figure size. The default is [6, 4].
title (str, optional) – Plot title. The default is ‘SSA scree plot’.
savefig_settings (dict) –
the dictionary of arguments for plt.savefig(); some notes below:
- ”path” must be specified; it can be any existed or non-existed path, with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”}
title_kwargs (dict) – the keyword arguments for ax.set_title()
ax (matplotlib.axis, optional) – the axis object from matplotlib See matplotlib.axes. for details.
xlim (list, optional) – x-axis limits. The default is None.
clr_mcssa (str, optional) – color of the Monte Carlo SSA AR(1) shading (if data are provided) default: red
clr_eig (str, optional) – color of the eigenvalues, default: black
clr_signif (str, optional) – color of the highlights for significant eigenvalue. (default: teal)

See also

pyleoclim.core.series.Series.ssa: Singular Spectrum Analysis for timeseries objects
pyleoclim.core.utils.decomposition.ssa: Singular Spectrum Analysis utility
pyleoclim.core.ssares.SsaRes.modeplot: plot SSA modes

Examples

Plot the SSA eigenvalue spectrum of the NINO3 index:

ts = pyleo.utils.load_dataset('NINO3')
ssa = ts.ssa()
fig, ax = ssa.screeplot()

Resolution (pyleoclim.Resolution)

class pyleoclim.core.resolutions.Resolution(resolution, time=None, resolution_unit=None, label=None, timeseries=None)[source]

Resolution object

Resolution objects store time axis resolution information derived from Series objects. They are generated via the resolution method applied to a Series object and contain methods relevant to the analysis of resolution information.

Parameters:

resolution (list or numpy.array) – An array containing information on the time step between subsequent values in the original time series.
time (list or numpy.array) – An array containing the time axis from the original series. This should be same length as resolution, so removing one value from the original axis is required.
resolution_unit (str) – Unit of resolution. This should be the same as the original time axis.
label (str) – Label for the resolution object.
timeseries (pyleoclim.Series) – Original pyleoclim timeseries object.

dashboard(figsize=[11, 8], title=None, plot_kwargs=None, histplot_kwargs=None, savefig_settings=None)[source]

Resolution plot dashboard

Parameters:

figsize (list or tuple, optional) – Figure size. The default is [11,8].
title (str) – Figure title
plot_kwargs (dict) – the dictionary of keyword arguments for ax.plot()
histplot_kwargs (dict, optional) – The dictionary of keyword arguments for ax.histplot()
savefig_settings (dict, optional) –
the dictionary of arguments for plt.savefig(); some notes below:
- ”path” must be specified; it can be any existed or non-existed path, with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”}.
The default is None.

Returns:

fig (matplotlib.figure) – The figure
ax (matplolib.axis) – The axis

Examples

ts = pyleo.utils.load_dataset('EDC-dD')
resolution = ts.resolution()
resolution.dashboard()

(<Figure size 1100x800 with 2 Axes>,
 {'res': <Axes: xlabel='Age [y BP]', ylabel='resolution [y BP]'>,
  'res_hist': <Axes: xlabel='Counts'>})

describe()[source]

Describe the stats of the time series

Returns:: stats – Dictionary of relevant stats produced by scipy.stats.describe
Return type:: dict

Examples

ts = pyleo.utils.load_dataset('EDC-dD')
resolution = ts.resolution()
resolution.describe()

{'nobs': np.int64(5784),
 'minmax': (np.float64(8.244210000000002), np.float64(1364.0)),
 'mean': np.float64(138.5932963710235),
 'variance': np.float64(29806.73648249974),
 'skewness': np.float64(2.661861461835658),
 'kurtosis': np.float64(8.705801510819656),
 'median': np.float64(58.132250000006024)}

histplot(figsize=[10, 4], title=None, savefig_settings=None, ax=None, ylabel='KDE', vertical=False, edgecolor='w', **plot_kwargs)[source]

Plot the distribution of the resolution values

Parameters:

figsize (list) – a list of two integers indicating the figure size
title (str) – the title for the figure
savefig_settings (dict) –
the dictionary of arguments for plt.savefig(); some notes below:
- ”path” must be specified; it can be any existed or non-existed path,
  with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”}
ax (matplotlib.axis, optional) – A matplotlib axis
ylabel (str) – Label for the count axis
vertical ({True,False}) – Whether to flip the plot vertically
edgecolor (matplotlib.color) – The color of the edges of the bar
plot_kwargs (dict) – Plotting arguments for seaborn histplot: https://seaborn.pydata.org/generated/seaborn.histplot.html

See also

pyleoclim.utils.plotting.savefig: saving figure in Pyleoclim

Examples

ts = pyleo.utils.load_dataset('EDC-dD')
res = ts.resolution()
res.histplot()

(<Figure size 1000x400 with 1 Axes>,
 <Axes: xlabel='resolution [y BP]', ylabel='KDE'>)

make_labels()[source]

Initialization of plot labels based on Series metadata

Returns:

time_header (str) – Label for the time axis
value_header (str) – Label for the value axis

plot(figsize=[10, 4], marker=None, markersize=None, color=None, linestyle=None, linewidth=None, xlim=None, ylim=None, label=None, xlabel=None, ylabel=None, title=None, zorder=None, legend=True, plot_kwargs=None, lgd_kwargs=None, alpha=None, savefig_settings=None, ax=None, invert_xaxis=False, invert_yaxis=False)[source]

Plot the timeseries

Parameters:

figsize (list) – a list of two integers indicating the figure size
marker (str) – e.g., ‘o’ for dots See [matplotlib.markers](https://matplotlib.org/stable/api/markers_api.html) for details
markersize (float) – the size of the marker
color (str, list) – the color for the line plot e.g., ‘r’ for red See [matplotlib colors](https://matplotlib.org/stable/gallery/color/color_demo.html) for details
linestyle (str) – e.g., ‘–’ for dashed line See [matplotlib.linestyles](https://matplotlib.org/stable/gallery/lines_bars_and_markers/linestyles.html) for details
linewidth (float) – the width of the line
label (str) – the label for the line
xlabel (str) – the label for the x-axis
ylabel (str) – the label for the y-axis
title (str) – the title for the figure
zorder (int) – The default drawing order for all lines on the plot
legend ({True, False}) – plot legend or not
invert_xaxis (bool, optional) – if True, the x-axis of the plot will be inverted
invert_yaxis (bool, optional) – same for the y-axis
plot_kwargs (dict) – the dictionary of keyword arguments for ax.plot() See [matplotlib.pyplot.plot](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.plot.html) for details
lgd_kwargs (dict) – the dictionary of keyword arguments for ax.legend() See [matplotlib.pyplot.legend](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.legend.html) for details
alpha (float) – Transparency setting
savefig_settings (dict) –
the dictionary of arguments for plt.savefig(); some notes below: - “path” must be specified; it can be any existed or non-existed path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”}
ax (matplotlib.axis, optional) – the axis object from matplotlib See [matplotlib.axes](https://matplotlib.org/api/axes_api.html) for details.

Returns:

fig (matplotlib.figure) – the figure object from matplotlib See [matplotlib.pyplot.figure](https://matplotlib.org/stable/api/figure_api.html) for details.
ax (matplotlib.axis) – the axis object from matplotlib See [matplotlib.axes](https://matplotlib.org/stable/api/axes_api.html) for details.

Notes

When ax is passed, the return will be ax only; otherwise, both fig and ax will be returned.

See also

pyleoclim.utils.plotting.savefig: saving a figure in Pyleoclim

Examples

ts = pyleo.utils.load_dataset('EDC-dD')
resolution = ts.resolution()
resolution.plot()

(<Figure size 1000x400 with 1 Axes>,
 <Axes: xlabel='Age [y BP]', ylabel='resolution [y BP]'>)

MultipleResolution (pyleoclim.MultipleResolution)

class pyleoclim.core.resolutions.MultipleResolution(resolution_list, time_unit)[source]

MultipleResolution object

MultipleResolution objects store time axis resolution information derived from MultipleSeries objects. They are generated via the resolution method applied to a MultipleSeries object and contain methods relevant to the analysis of resolution information.

Parameters:

resolution_list (list) – List of resolution objects.
time_unit (str) – The unit of time for the resolution.

See also

Resolution: The base class from which MultipleResolution is derived.

describe()[source]

Describe the stats of the resolution list

Returns:

resolution_dict – Dictionary of relevant stats produced by scipy.stats.describe

Return type:

dict

Examples

..jupyter-execute:

co2ts = pyleo.utils.load_dataset('AACO2')
edc = pyleo.utils.load_dataset('EDC-dD')
ms = edc & co2ts # create MS object
ms_resolution = ms.resolution(statistic=None)
ms_resolution.describe()

summary_plot(figsize=(10, 8), xlabel=None, ylabel=None, legend=False, ax=None, log_scale=None, boxplot_whis=[0, 100], boxplot_width=0.6, boxplot_dodge=False, boxplot_palette='viridis', stripplot_size=2, stripplot_color='.3', stripplot_alpha=0.8, stripplot_dodge=False, boxplot_kwargs=None, stripplot_kwargs=None, savefig_settings=None)[source]

Summary plot showing distribution of resolutions from each resolution object as box plots with overlaid strip plots.

Parameters:

figsize (tuple, list) – Size of the figure.
xlabel (str) – Label for the x axis. “Resolution [{time_unit}]” will be used by default.
ylabel (str) – Label for the y axis. Left empty by default.
legend (bool) – Whether or not to plot the legend. Default is False.
ax (matplotlib.ax) – The matplotlib axis onto which to return the figure. The default is None.
log_scale (bool) – Whether to plot the y-axis on a log scale. Default is None.
boxplot_whis (float or pair of floats) – If scalar, whiskers are drawn to the farthest datapoint within whis * IQR from the nearest hinge. If a tuple, it is interpreted as percentiles that whiskers represent. Default is [0,100]
boxplot_width (float) – Width allotted to each element on the orient axis.
boxplot_dodge (bool) – Whether boxplot elements should be narrowed and shifted along the orient axis to eliminate overlap. Default is False.
boxplot_palette (palette name, list, or dict) – Colors to use for the different levels of the hue variable. Should be something that can be interpreted by color_palette(), or a dictionary mapping hue levels to matplotlib colors. Gets passed to seaborn.boxplot for details.
stripplot_size (float) – Size of stripplot markers.
stripplot_color (str) – Color for the stripplot markers.
stripplot_alpha (float) – Alpha for the stripplot markers.
stripplot_dodge (bool) – Whether stripplot elements should be narrowed and shifted along the orient axis to eliminate overlap. Default is False.
boxplot_kwargs (dict) – Dictionary of arguments for seaborn.boxplot. Arguments that are passed here will overwrite explicit arguments (e.g. whis, width, etc.).
stripplot_kwargs (dict) – Dictionary of argument for seaborn.stripplot. Arguments that are passed here will overwrite explicit arguments (e.g. size, color, etc.).
savefig_settings (dictionary, optional) –
the dictionary of arguments for pyleo.utils.plotting.savefig(); some notes below: - “path” must be specified; it can be any existing or non-existing path,

with or without a suffix; if the suffix is not given in “path”, it will follow “format”
- ”format” can be one of {“pdf”, “eps”, “png”, “ps”} The default is None.

Returns:

fig (matplotlib.figure) – the figure object from matplotlib See [matplotlib.pyplot.figure](https://matplotlib.org/stable/api/figure_api.html) for details.
ax (matplotlib.axis) – the axis object from matplotlib See [matplotlib.axes](https://matplotlib.org/stable/api/axes_api.html) for details.

Examples

co2ts = pyleo.utils.load_dataset('AACO2')
edc = pyleo.utils.load_dataset('EDC-dD')
ms = edc & co2ts # create MS object
ms_resolution = ms.resolution(statistic=None)
ms_resolution.summary_plot()

Time unit not found, attempting conversion.
Converted to ky BP

(<Figure size 1000x800 with 1 Axes>, <Axes: xlabel='Resolution [ky BP]'>)