diff --git a/cfr/proxy.py b/cfr/proxy.py index 8edda375..27c6e7ec 100644 --- a/cfr/proxy.py +++ b/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): diff --git a/cfr/utils.py b/cfr/utils.py index 39c32183..f868bfcf 100644 --- a/cfr/utils.py +++ b/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 diff --git a/tests/test_proxy.py b/tests/test_proxy.py index 450c1715..caf2b086 100644 --- a/tests/test_proxy.py +++ b/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() diff --git a/tests/test_utils.py b/tests/test_utils.py index 2239b953..7994cea5 100644 --- a/tests/test_utils.py +++ b/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