diff --git a/pyleoclim/core/series.py b/pyleoclim/core/series.py index 21976d27..4b8e6e84 100644 --- a/pyleoclim/core/series.py +++ b/pyleoclim/core/series.py @@ -7,6 +7,8 @@ How to create and manipulate such objects is described in a short example below, while `this notebook `_ demonstrates how to apply various Pyleoclim methods to Series objects. """ +import operator + from ..utils import tsutils, plotting, tsmodel, tsbase, mapping, lipdutils from ..utils import wavelet as waveutils from ..utils import spectral as specutils @@ -163,6 +165,55 @@ def __init__(self, time, value, time_name=None, time_unit=None, value_name=None, self.mean=np.mean(self.value) else: self.mean = mean + + @property + def datetime_index(self): + datum, exponent, direction = tsutils.time_unit_to_datum_exp_dir(self.time_unit) + if direction == 'prograde': + op = operator.add + elif direction == 'retrograde': + op = operator.sub + else: + raise ValueError(f'Expected one of {"prograde", "retrograde"}, got {direction}') + + timedelta = self.time * 10**exponent + years = timedelta.astype('int') + seconds = ((timedelta - timedelta.astype('int')) * tsutils.SECONDS_PER_YEAR).astype('timedelta64[s]') + + np_times = op(op(int(datum), years).astype(str).astype('datetime64[s]'), seconds) + return pd.DatetimeIndex(np_times, name=self.time_name) + + @property + def metadata(self): + return dict( + time_unit = self.time_unit, + value_unit = self.value_unit, + label = self.label, + ) + + @classmethod + def from_pandas(cls, ser, metadata): + time = tsutils.convert_datetime_index_to_time(ser.index, metadata['time_unit']) + return cls( + time=time, + value=ser.to_numpy(), + time_name=ser.index.name, + value_name=ser.name, + **metadata, + ) + + def to_pandas(self): + ser = pd.Series(self.value, index=self.datetime_index, name=self.value_name) + # Could be a dataclass instead? + return (ser, self.metadata) + + def pandas_method(self, method): + ser, metadata = self.to_pandas() + result = method(ser) + if not isinstance(result, pd.Series): + raise ValueError('Given method does not return a pandas Series and cannot be applied') + return self.from_pandas(result, metadata) + def convert_time_unit(self, time_unit='years', keep_log=False): ''' Convert the time unit of the Series object diff --git a/pyleoclim/utils/tsutils.py b/pyleoclim/utils/tsutils.py index 9826f869..7be17214 100644 --- a/pyleoclim/utils/tsutils.py +++ b/pyleoclim/utils/tsutils.py @@ -46,6 +46,59 @@ clean_ts ) +SECONDS_PER_YEAR = 365.25 * 60 * 60 * 24 + +def time_unit_to_datum_exp_dir(time_unit): + if time_unit == 'years': + datum = '0' + exponent = 0 + direction = 'prograde' + elif time_unit in ('ky', 'kyr', 'kyrs', 'kiloyear', 'ka'): + datum = '1950' + exponent = 3 + direction = 'retrograde' + elif time_unit in ('y', 'yr', 'yrs', 'year', 'year CE', 'years CE', 'year(s) AD'): + datum = '0' + exponent = 0 + direction = 'prograde' + elif time_unit == 'BP': + datum = '1950' + exponent = 0 + direction = 'retrograde' + elif time_unit in ('yr BP', 'yrs BP', 'years BP'): + datum ='1950' + exponent = 0 + direction = 'retrograde' + elif time_unit in ('Ma', 'My'): + datum ='1950' + exponent = 6 + direction = 'retrograde' + elif time_unit in ('Ga', 'Gy'): + datum ='1950' + exponent = 3 + direction = 'retrograde' + else: + raise ValueError(f'Time unit {time_unit} not supported') + + return (datum, exponent, direction) + +def convert_datetime_index_to_time(datetime_index, time_unit): + datum, exponent, direction = time_unit_to_datum_exp_dir(time_unit) + import operator + if direction == 'prograde': + multiplier = 1 + elif direction == 'retrograde': + multiplier = -1 + else: + raise ValueError(f'Expected one of {"prograde", "retrograde"}, got {direction}') + year_diff = (datetime_index.year - int(datum)) + seconds_diff = (datetime_index.to_numpy() - datetime_index.year.astype(str).astype('datetime64[s]').to_numpy()).astype('int') + diff = year_diff + seconds_diff / SECONDS_PER_YEAR + time = multiplier * diff / 10**exponent + + return time + + def simple_stats(y, axis=None): """ Computes simple statistics