Merge pull request #449 from LinkedEarth/modeplot_PCunits

Modeplot pc units

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

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)
 

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,25 +17,36 @@
 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
 
     '''
     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:

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:

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,25 +247,24 @@ 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])
         ax.set_title('Analyzing function')
         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:

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'