Add days_in_year and decimal_year to dt accessor

doc/api.rst

@@ -527,9 +527,11 @@ Datetimelike properties
    DataArray.dt.quarter
    DataArray.dt.days_in_month
    DataArray.dt.daysinmonth
+   DataArray.dt.days_in_year
    DataArray.dt.season
    DataArray.dt.time
    DataArray.dt.date
+   DataArray.dt.decimal_year
    DataArray.dt.calendar
    DataArray.dt.is_month_start
    DataArray.dt.is_month_end

doc/whats-new.rst

@@ -23,6 +23,9 @@ v2024.05.1 (unreleased)
 New Features
 ~~~~~~~~~~~~
 
+- Add :py:property:`~core.accessor_dt.DatetimeAccessor.days_in_year` and :py:property:`~core.accessor_dt.DatetimeAccessor.decimal_year` to the Datetime accessor on DataArrays. (:pull:`9105`).
+  By `Pascal Bourgault <https://github.com/aulemahal>`_.
+
 Performance
 ~~~~~~~~~~~
 

xarray/coding/calendar_ops.py

@@ -3,6 +3,8 @@
 import numpy as np
 import pandas as pd
 
+from xarray.core.common import full_like
+from xarray.core.computation import where
 from xarray.coding.cftime_offsets import date_range_like, get_date_type
 from xarray.coding.cftimeindex import CFTimeIndex
 from xarray.coding.times import _should_cftime_be_used, convert_times
@@ -22,14 +24,35 @@
 ]
 
 
-def _days_in_year(year, calendar, use_cftime=True):
-    """Return the number of days in the input year according to the input calendar."""
-    date_type = get_date_type(calendar, use_cftime=use_cftime)
-    if year == -1 and calendar in _CALENDARS_WITHOUT_YEAR_ZERO:
-        difference = date_type(year + 2, 1, 1) - date_type(year, 1, 1)
+def _leap_gregorian(years):
+    # A year is a leap year if either (i) it is divisible by 4 but not by 100 or (ii) it is divisible by 400
+    return ((years % 4 == 0) & (years % 100 != 0)) | (years % 400 == 0)
+
+
+def _leap_julian(years):
+    # A year is a leap year if it is divisible by 4, even if it is also divisible by 100
+    return years % 4 == 0
+
+
+def _days_in_year(years, calendar):
+    """The number of days in each year for the corresponding calendar."""
+    if calendar in ['standard', 'gregorian']:
+        return 365 + where(years < 1582, _leap_julian(years), _leap_gregorian(years)) * 1
+    if calendar == 'proleptic_gregorian':
+        return 365 + _leap_gregorian(years) * 1
+    if calendar == 'julian':
+        return 365 + _leap_julian(years) * 1
+    if calendar in ['noleap', '365_day']:
+        const = 365
+    elif calendar in ['all_leap', '366_day']:
+        const = 366
+    elif calendar == '360_day':
+        const = 360
     else:
-        difference = date_type(year + 1, 1, 1) - date_type(year, 1, 1)
-    return difference.days
+        raise ValueError(f'Calendar {calendar} not recognized.')
+    if isinstance(years, (float, int)):
+        return const
+    return full_like(years, const)
 
 
 def convert_calendar(
@@ -188,11 +211,7 @@ def convert_calendar(
         # Special case for conversion involving 360_day calendar
         if align_on == "year":
             # Instead of translating dates directly, this tries to keep the position within a year similar.
-            new_doy = time.groupby(f"{dim}.year").map(
-                _interpolate_day_of_year,
-                target_calendar=calendar,
-                use_cftime=use_cftime,
-            )
+            new_doy = _interpolate_day_of_year(time, target_calendar=calendar)
         elif align_on == "random":
             # The 5 days to remove are randomly chosen, one for each of the five 72-days periods of the year.
             new_doy = time.groupby(f"{dim}.year").map(
@@ -232,16 +251,13 @@ def convert_calendar(
     return out
 
 
-def _interpolate_day_of_year(time, target_calendar, use_cftime):
-    """Returns the nearest day in the target calendar of the corresponding
-    "decimal year" in the source calendar.
-    """
-    year = int(time.dt.year[0])
-    source_calendar = time.dt.calendar
+def _interpolate_day_of_year(times, target_calendar):
+    """Returns the nearest day in the target calendar of the corresponding "decimal year" in the source calendar."""
+    source_calendar = times.dt.calendar
     return np.round(
-        _days_in_year(year, target_calendar, use_cftime)
-        * time.dt.dayofyear
-        / _days_in_year(year, source_calendar, use_cftime)
+        _days_in_year(times.dt.year, target_calendar)
+        * times.dt.dayofyear
+        / _days_in_year(times.dt.year, source_calendar)
     ).astype(int)
 
 
@@ -250,18 +266,18 @@ def _random_day_of_year(time, target_calendar, use_cftime):
 
     Removes Feb 29th and five other days chosen randomly within five sections of 72 days.
     """
-    year = int(time.dt.year[0])
+    year = time.dt.year[0]
     source_calendar = time.dt.calendar
     new_doy = np.arange(360) + 1
     rm_idx = np.random.default_rng().integers(0, 72, 5) + 72 * np.arange(5)
     if source_calendar == "360_day":
         for idx in rm_idx:
             new_doy[idx + 1 :] = new_doy[idx + 1 :] + 1
-        if _days_in_year(year, target_calendar, use_cftime) == 366:
+        if _days_in_year(year, target_calendar) == 366:
             new_doy[new_doy >= 60] = new_doy[new_doy >= 60] + 1
     elif target_calendar == "360_day":
         new_doy = np.insert(new_doy, rm_idx - np.arange(5), -1)
-        if _days_in_year(year, source_calendar, use_cftime) == 366:
+        if _days_in_year(year, source_calendar) == 366:
             new_doy = np.insert(new_doy, 60, -1)
     return new_doy[time.dt.dayofyear - 1]
 
@@ -304,10 +320,11 @@ def _datetime_to_decimal_year(times, dim="time", calendar=None):
     """
     from xarray.core.dataarray import DataArray
 
-    calendar = calendar or times.dt.calendar
+    if calendar is None:
+        calendar = times.dt.calendar
 
     if is_np_datetime_like(times.dtype):
-        times = times.copy(data=convert_times(times.values, get_date_type("standard")))
+        times = times.copy(data=convert_times(times.values, get_date_type("proleptic_gregorian")))
 
     def _make_index(time):
         year = int(time.dt.year[0])

xarray/core/accessor_dt.py

@@ -7,6 +7,7 @@
 import pandas as pd
 
 from xarray.coding.times import infer_calendar_name
+from xarray.coding.calendar_ops import _days_in_year, _datetime_to_decimal_year
 from xarray.core import duck_array_ops
 from xarray.core.common import (
     _contains_datetime_like_objects,
@@ -533,6 +534,26 @@ def calendar(self) -> CFCalendar:
         """
         return infer_calendar_name(self._obj.data)
 
+    @property
+    def days_in_year(self) -> T_DataArray:
+        """The number of days in the year."""
+        obj_type = type(self._obj)
+        result = _days_in_year(self.year, self.calendar)
+        return obj_type(
+            result, name="days_in_year", coords=self._obj.coords,
+            dims=self._obj.dims, attrs=self._obj.attrs
+        )
+
+    @property
+    def decimal_year(self) -> T_DataArray:
+        """Convert the dates as a fractional year."""
+        obj_type = type(self._obj)
+        result = _datetime_to_decimal_year(self._obj)
+        return obj_type(
+            result, name="decimal_year", coords=self._obj.coords,
+            dims=self._obj.dims, attrs=self._obj.attrs
+        )
+
 
 class TimedeltaAccessor(TimeAccessor[T_DataArray]):
     """Access Timedelta fields for DataArrays with Timedelta-like dtypes.

xarray/tests/test_accessor_dt.py

@@ -141,6 +141,21 @@ def test_strftime(self) -> None:
             "2000-01-01 01:00:00" == self.data.time.dt.strftime("%Y-%m-%d %H:%M:%S")[1]
         )
 
+    @requires_cftime
+    @pytest.mark.parametrize("calendar,expected", [('standard', 366), ('noleap', 365), ('360_day', 360), ('all_leap', 366)])
+    def test_days_in_year(self, calendar, expected) -> None:
+        assert (
+            self.data.convert_calendar(calendar, align_on='year').time.dt.days_in_year == expected
+        ).all()
+
+    @requires_cftime
+    def test_decimal_year(self) -> None:
+        h_per_yr = 366 * 24
+        np.testing.assert_array_equal(
+            self.data.time.dt.decimal_year[0:3],
+            [2000, 2000 + 1 / h_per_yr, 2000 + 2 / h_per_yr]
+        )
+
     def test_not_datetime_type(self) -> None:
         nontime_data = self.data.copy()
         int_data = np.arange(len(self.data.time)).astype("int8")