diff --git a/environment.yml b/environment.yml index 9d4164f1..7f149585 100644 --- a/environment.yml +++ b/environment.yml @@ -1,8 +1,8 @@ -# Development environment. This environment includes all the required +# Development environment. This environment includes all the required # packages for this repo including running the test suite. name: pyleo channels: - - LinkedEarth + - LinkedEarth - conda-forge - default dependencies: @@ -18,7 +18,7 @@ dependencies: - pyyaml - requests - scikit-learn - - scipy + - scipy=1.10.1 - seaborn - statsmodels - tabulate diff --git a/pyleoclim/core/multivardecomp.py b/pyleoclim/core/multivardecomp.py index 75c61108..940db5af 100644 --- a/pyleoclim/core/multivardecomp.py +++ b/pyleoclim/core/multivardecomp.py @@ -1,14 +1,13 @@ import numpy as np -import pandas as pd +#import pandas as pd from matplotlib import pyplot as plt, gridspec -from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable +#from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable from matplotlib.ticker import MaxNLocator -import cartopy.crs as ccrs -import cartopy.feature as cfeature +#import cartopy.crs as ccrs +#import cartopy.feature as cfeature from ..core import series -from ..utils import plotting, lipdutils -from ..utils import mapping +from ..utils import plotting, mapping, tsbase class MultivariateDecomp: @@ -306,7 +305,17 @@ def modeplot(self, index=0, figsize=[8, 8], fig=None, savefig_settings=None,gs=N ax['pc'] = fig.add_subplot(gs[0, 0]) label = rf'$PC_{index + 1}$' t = self.orig.series_list[0].time - ts = series.Series(time=t, value=PC, verbose=False) # define timeseries object for the PC + # get time unit + if self.orig.time_unit is not None: + time_unit = self.orig.time_unit + else: + time_unit = self.orig.series_list[0].time_unit + + time_name, _ = tsbase.disambiguate_time_metadata(time_unit) + + ts = series.Series(time=t, value=PC, verbose=False, + time_name=time_name, + time_unit=time_unit) # define timeseries object for the PC ts.plot(ax=ax['pc']) ax['pc'].set_ylabel(label) diff --git a/pyleoclim/core/psds.py b/pyleoclim/core/psds.py index 97515a93..841bbb97 100644 --- a/pyleoclim/core/psds.py +++ b/pyleoclim/core/psds.py @@ -1,6 +1,7 @@ -from ..utils import plotting, lipdutils +from ..utils import plotting from ..utils import wavelet as waveutils from ..utils import spectral as specutils +#from ..utils import tsbase #from ..core.multiplepsd import * #import ..core.multiplepsd @@ -16,6 +17,27 @@ from matplotlib.ticker import ScalarFormatter, FormatStrFormatter +# def infer_period_unit_from_time_unit(time_unit): +# ''' infer a period unit based on the given time unit + +# ''' +# if time_unit is None: +# period_unit = None +# else: +# unit_group = lipdutils.timeUnitsCheck(time_unit) +# if unit_group != 'unknown': +# if unit_group == 'kage_units': +# period_unit = 'kyrs' +# else: +# period_unit = 'yrs' +# else: +# if time_unit[-1] == 's': +# period_unit = time_unit +# else: +# period_unit = f'{time_unit}s' + +# return period_unit + def infer_period_unit_from_time_unit(time_unit): ''' infer a period unit based on the given time unit @@ -23,18 +45,8 @@ def infer_period_unit_from_time_unit(time_unit): if time_unit is None: period_unit = None else: - unit_group = lipdutils.timeUnitsCheck(time_unit) - if unit_group != 'unknown': - if unit_group == 'kage_units': - period_unit = 'kyrs' - else: - period_unit = 'yrs' - else: - if time_unit[-1] == 's': - period_unit = time_unit - else: - period_unit = f'{time_unit}s' - + tu = time_unit.lower().replace(".","").split() + period_unit = tu[0] return period_unit class PSD: diff --git a/pyleoclim/core/series.py b/pyleoclim/core/series.py index 169427dd..7333485f 100644 --- a/pyleoclim/core/series.py +++ b/pyleoclim/core/series.py @@ -144,14 +144,19 @@ def __init__(self, time, value, time_unit=None, time_name=None, time = np.array(time) value = np.array(value) - # assign time metadata if they are not provided + # assign time metadata if they are not provided or provided incorrectly + offending = [tsbase.MATCH_CE, tsbase.MATCH_BP] + if time_unit is None: time_unit='years CE' if verbose: warnings.warn(f'No time_unit parameter provided. Assuming {time_unit}.', UserWarning) + elif time_unit.lower().replace(".","") in frozenset().union(*offending): + # fix up time name and units for offending cases + time_name, time_unit = tsbase.disambiguate_time_metadata(time_unit) else: # give a proper time name to those series that confuse that notion with time units - time_name, time_unit = tsbase.disambiguate_time_metadata(time_unit) + time_name, _ = tsbase.disambiguate_time_metadata(time_unit) if time_name is None: if verbose: @@ -1257,7 +1262,7 @@ def ssa(self, M=None, nMC=0, f=0.3, trunc = None, var_thresh=80, online=True): res = decomposition.ssa(self.value, M=M, nMC=nMC, f=f, trunc = trunc, var_thresh=var_thresh, online=online) - resc = SsaRes(label=self.label, original=self.value, time = self.time, eigvals = res['eigvals'], eigvecs = res['eigvecs'], + resc = SsaRes(label=self.label, orig=self, eigvals = res['eigvals'], eigvecs = res['eigvecs'], pctvar = res['pctvar'], PC = res['PC'], RCmat = res['RCmat'], RCseries=res['RCseries'], mode_idx=res['mode_idx']) if nMC >= 0: diff --git a/pyleoclim/core/ssares.py b/pyleoclim/core/ssares.py index 44c6c907..b9d9cf82 100644 --- a/pyleoclim/core/ssares.py +++ b/pyleoclim/core/ssares.py @@ -24,10 +24,7 @@ class SsaRes: Parameters ---------- - time : array - time axis - - original : array + orig: Series timeseries on which SSA was performed eigvals: float (M, 1) @@ -59,9 +56,8 @@ class SsaRes: pyleoclim.utils.decomposition.ssa : Singular Spectrum Analysis ''' - def __init__(self, time, original, label, eigvals, eigvecs, pctvar, PC, RCmat, RCseries,mode_idx, eigvals_q=None): - self.time = time - self.original = original + def __init__(self, orig, label, eigvals, eigvecs, pctvar, PC, RCmat, RCseries,mode_idx, eigvals_q=None): + self.orig = orig self.label = label self.eigvals = eigvals self.eigvals_q = eigvals_q @@ -251,14 +247,14 @@ def modeplot(self, index=0, figsize=[10, 5], savefig_settings=None, gs = gridspec.GridSpec(2, 2) # plot RC ax = fig.add_subplot(gs[0, :]) - - ax.plot(self.time,RC,label='mode '+str(index+1),zorder=99) + ts_rc = self.orig.copy() + ts_rc.value = RC + ts_rc.label = r'$RC_'+str(index+1)+'$' + ts_rc.plot(zorder=99,ax=ax) if plot_original: - ax.plot(self.time,self.original,color='Silver',lw=1,label='original') - ax.legend() - ax.set_xlabel('Time') - ax.set_ylabel(r'$RC_'+str(index+1)+'$') - ax.set_title('SSA reconstructed component '+str(index+1)+', '+ '{:3.2f}'.format(self.pctvar[index]) + '% variance explained',weight='bold') + self.orig.plot(ax=ax,color='Silver',linewidth=1) + + ax.set_title('SSA mode '+str(index+1)+', '+ '{:3.2f}'.format(self.pctvar[index]) + '% variance explained',weight='bold') # plot T-EOF ax = fig.add_subplot(gs[1, 0]) ax.plot(self.eigvecs[:,index]) @@ -266,10 +262,9 @@ def modeplot(self, index=0, figsize=[10, 5], savefig_settings=None, ax.set_xlabel('Time'), ax.set_ylabel('T-EOF') # plot spectrum ax = fig.add_subplot(gs[1, 1]) - ts_rc = series.Series(time=self.time, value=RC, verbose=False) # define timeseries object for the RC psd_mtm_rc = ts_rc.interp().spectral(method=spec_method) psd_mtm_rc.plot(ax=ax) - ax.set_xlabel('Period') + #ax.set_xlabel('Period') ax.set_title('RC Spectrum ('+spec_method+')') if 'path' in savefig_settings: diff --git a/pyleoclim/data/metadata.yml b/pyleoclim/data/metadata.yml index 2cc20a91..893f7d8c 100644 --- a/pyleoclim/data/metadata.yml +++ b/pyleoclim/data/metadata.yml @@ -136,7 +136,7 @@ cenogrid_d18O: pyleo_kwargs: time_unit: 'My BP' time_name: 'Tuned Time' - value_unit: '‰ (VPDB)' + value_unit: '‰ VPDB' value_name: '$\delta^{18} \mathrm{O}$' label: 'CENOGRID' archiveType: 'Marine Sediment' @@ -152,7 +152,7 @@ cenogrid_d13C: pyleo_kwargs: time_unit: 'My BP' time_name: 'Tuned Time' - value_unit: '‰ (VPDB)' + value_unit: '‰ VPDB' value_name: '$\delta^{13} \mathrm{C}$' label: 'CENOGRID' archiveType: 'Marine Sediment'