annualize: support years at or before 0; always warn about dropped records

cfr/proxy.py

@@ -1729,7 +1729,10 @@ def annualize(self, months=list(range(1, 13)), force=False, verbose=False):
             new += spobj
 
         new = new.filter(by='tag', keys=['annualized'])
-        if verbose and not force: p_warning(f'{self.nrec - new.nrec} records have been dropped as they cannot be properly annualized with the given months.')
+        dropped = [pid for pid in self.records if pid not in new.records]
+        if len(dropped) > 0:
+            # warn even when not verbose, since the dropped records are otherwise easy to miss
+            p_warning(f'{len(dropped)} records have been dropped as they cannot be properly annualized with the given months: {", ".join(dropped)}')
         return new
 
     def slice(self, timespan):

cfr/utils.py

@@ -219,15 +219,43 @@ def gcd(lat1, lon1, lat2, lon2):
 	km = 6367 * c
 	return km
 
+def _year_float2year_month(year_float):
+    ''' Split floats in unit of year into integer years and months, as year_float2datetime does,
+    but without datetime objects, so that years <= 0 are supported
+    '''
+    year_float = np.asarray(year_float, dtype=float)
+    year = np.floor(year_float).astype(int)
+    is_leap = (year % 4 == 0) & ((year % 100 != 0) | (year % 400 == 0))
+    year_length = np.where(is_leap, 366, 365)
+    # days since 1 January; one minute is added to fix the numerical error, as in year_float2datetime
+    doy = np.floor((year_float - year)*year_length + 1/1440).astype(int)
+    rollover = doy >= year_length
+    year = np.where(rollover, year+1, year)
+    doy = np.where(rollover, 0, doy)
+    month_starts = np.cumsum([0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30])
+    month = np.searchsorted(month_starts, doy, side='right')
+    month_starts_leap = month_starts + (month_starts >= 59)
+    month_leap = np.searchsorted(month_starts_leap, doy, side='right')
+    month = np.where(is_leap & ~rollover, month_leap, month)
+    return year, month
+
 def annualize(time=None, value=None, da=None, months=None):
     # convert negative months to support expressions like [-12, 1, 2]
     months = list(range(1, 13)) if months is None else np.abs(months)
 
     if da is None:
-        dates = year_float2datetime(time)
-        da_wk = xr.DataArray(
-            value, dims=['time'], coords={'time': dates}
-        )
+        # bin on numeric years rather than calendar dates, so that years <= 0 are supported;
+        # a month after the last month of the season belongs to the next year's season (e.g., Dec for [12, 1, 2])
+        year, month = _year_float2year_month(time)
+        mask = np.isin(month, months)
+        if not np.any(mask):
+            raise ValueError(f'No data in months {list(months)}.')
+
+        season_year = year[mask] + (month[mask] > months[-1])
+        sr = pd.Series(np.asarray(value, dtype=float)[mask]).groupby(season_year).mean()
+        time_ann = np.arange(sr.index.min(), sr.index.max()+1).astype(float)
+        value_ann = sr.reindex(time_ann.astype(int)).values
+        return time_ann, value_ann
     else:
         da_wk = da
 
@@ -236,13 +264,8 @@ def annualize(time=None, value=None, da=None, months=None):
     anchor = ['JAN', 'FEB', 'MAR', 'APR', 'MAY', 'JUN', 'JUL', 'AUG', 'SEP', 'OCT', 'NOV', 'DEC']
     idx = months[-1]-1
     sda_ann = sda.resample(time=f'A-{anchor[idx]}').mean()
-    time_ann = np.floor(datetime2year_float(sda_ann.time.values))
-    value_ann = sda_ann.values
 
-    if da is None:
-        return time_ann, value_ann
-    else:
-        return sda_ann
+    return sda_ann
 
 def annualize_var(year_float, var, resolution='month', weights=None):
     ''' Annualize a variable array

tests/test_proxy.py

@@ -411,6 +411,14 @@ def test_annualize_returns_database(self):
         assert isinstance(ann, ProxyDatabase)
         assert ann.nrec <= pdb.nrec  # might drop records if can't annualize
 
+    def test_annualize_keeps_records_before_year_zero(self):
+        pdb = ProxyDatabase()
+        for i, start in enumerate([1.5, 0.5, -4.5]):
+            pdb += ProxyRecord(pid=f'early_{i}', lat=0.0, lon=0.0, ptype='tree.TRW',
+                               time=start + np.arange(50), value=np.arange(50, dtype=float))
+        ann = pdb.annualize()
+        assert ann.nrec == 3
+
 
 # ---------------------------------------------------------------------------
 # ProxyRecord.get_clim() / del_clim() / del_pseudo()

tests/test_utils.py

@@ -187,6 +187,35 @@ def test_jja_season(self):
         # JJA mean should be 1.0 since we only picked JJA months
         np.testing.assert_allclose(val_ann, 1.0, atol=1e-10)
 
+    @pytest.mark.parametrize('start', [1.5, 0.5, -4.5, -1161.0])
+    def test_years_at_or_before_zero(self, start):
+        # annual records starting at or before year 0 used to raise 'year 0 is out of range'
+        time = start + np.arange(50)
+        value = np.arange(50, dtype=float)
+        time_ann, val_ann = utils.annualize(time, value)
+        np.testing.assert_array_equal(time_ann, np.floor(start) + np.arange(50))
+        np.testing.assert_allclose(val_ann, value)
+
+    def test_subannual_across_year_zero(self):
+        time = -2 + np.arange(48)/12
+        value = np.repeat(np.arange(4, dtype=float), 12)
+        time_ann, val_ann = utils.annualize(time, value)
+        np.testing.assert_array_equal(time_ann, [-2, -1, 0, 1])
+        np.testing.assert_allclose(val_ann, [0, 1, 2, 3])
+
+    def test_djf_season_across_year_zero(self):
+        # December goes to the following year's DJF
+        time = -1 + np.arange(24)/12
+        value = np.tile(np.arange(1, 13, dtype=float), 2)
+        time_ann, val_ann = utils.annualize(time, value, months=[-12, 1, 2])
+        np.testing.assert_array_equal(time_ann, [-1, 0, 1])
+        np.testing.assert_allclose(val_ann, [1.5, (12+1+2)/3, 12])
+
+    def test_missing_months_raise(self):
+        time = 1990 + np.arange(10, dtype=float)  # whole years fall in January
+        with pytest.raises(ValueError):
+            utils.annualize(time, np.ones(10), months=[6, 7, 8])
+
 
 # ---------------------------------------------------------------------------
 # colored_noise