refactor(scripts): replace dayofyear with day_of_year (pandas)

CHANGELOG-v0.md

@@ -87,6 +87,8 @@ This file covers every 0.x release; version 1.0 starts `CHANGELOG-v1.md`. The fo
 - The Fortran sources use one comment format: `!--- text` for a comment line, a comment of
   several lines between two `!---~---` lines, and a plain `! text` after a statement. Results are
   bit for bit; only comments change. (#48)
+- The Python tools and the phenology example use pandas' `day_of_year`, not its deprecated
+  `dayofyear`, which pandas 4 warns about and will remove. Results are bit for bit. (#52)
 - Each tier's `file_chunk` is honoured as written. A tier with a cohort or patch variable used to be
   cut to monthly files (sub-daily: daily), so the default monthly output was a file a month; it is now
   a file a year, its cohort and patch axes growing to the file's largest live count. (#55)

examples/example02_canopy_phenology/run_phenology.py

@@ -51,11 +51,11 @@ def crossing_day(series, first_doy, last_doy):
     """Per year, the (interpolated) day of year the series first rises through 0.5."""
     out = {}
     for year, g in series.groupby(series.index.year):
-        g = g[(g.index.dayofyear >= first_doy) & (g.index.dayofyear <= last_doy)]
+        g = g[(g.index.day_of_year >= first_doy) & (g.index.day_of_year <= last_doy)]
         k = np.flatnonzero(np.diff(np.sign(g.values - 0.5)) > 0)
         if len(k):
             i = k[0]
-            x0, x1 = g.index[i].dayofyear, g.index[i + 1].dayofyear
+            x0, x1 = g.index[i].day_of_year, g.index[i + 1].day_of_year
             out[year] = x0 + (0.5 - g.values[i]) * (x1 - x0) / (g.values[i + 1] - g.values[i])
     return pd.Series(out)
 
@@ -133,7 +133,7 @@ def simulate(self, values):
         tair = self.tair if self.tair is not None else np.full(n, 298.15)
         psi = self.psi if self.psi is not None else np.zeros(n)
         days = [ph.step(temp_day=t, par_hours=h, predawn_leaf_psi=w, doy=d)
-                for t, h, w, d in zip(tair, hours, psi, self.dates.dayofyear)]
+                for t, h, w, d in zip(tair, hours, psi, self.dates.day_of_year)]
         return pd.DataFrame({"cover": [x.leaf_cover for x in days],
                              "litter": [x.senescence + x.background for x in days],
                              "flush": [x.leaf_flush_tendency for x in days],
@@ -210,7 +210,7 @@ def __init__(self):
         lai = lai.dropna().rolling(3, center=True).median().dropna().loc[f"{y0 + 1}":f"{y1}"]
         rel = []
         for _, g in lai.groupby(lai.index.year):
-            doy = g.index.dayofyear
+            doy = g.index.day_of_year
             winter = g[(doy <= 90) | (doy >= 330)].median()
             summer = g[(doy >= 170) & (doy <= 240)].quantile(0.9)
             rel.append(((g - winter) / (summer - winter)).clip(-0.1, 1.2))
@@ -241,7 +241,7 @@ def free(self):
     def fallen(self, sim):
         """Fraction of the year's peak canopy lost, from August on (as HF003 counts it)."""
         peak = sim.cover.groupby(sim.index.year).transform("max")
-        return (1.0 - sim.cover / peak.clip(lower=1e-6)).where(sim.index.dayofyear >= 213, 0.0)
+        return (1.0 - sim.cover / peak.clip(lower=1e-6)).where(sim.index.day_of_year >= 213, 0.0)
 
     def basket_fraction(self, sim):
         """Cumulative fraction of each leaf year's model litter (from 1 April) at the basket dates."""

scripts/calibrate_fast/data_rules.py

@@ -265,7 +265,7 @@ def candidate_windows(scores: pd.Series, days: int, per_year: int, leaf_on: list
         end = start + pd.Timedelta(days=days - 1)
         if end.year != start.year:
             continue
-        s0, s1 = slot_of(start.dayofyear, per_year), slot_of(end.dayofyear, per_year)
+        s0, s1 = slot_of(start.day_of_year, per_year), slot_of(end.day_of_year, per_year)
         span = pd.date_range(start, end, freq="D")
         if s0 != s1 or not set(span.month) <= set(leaf_on):
             continue

scripts/calibrate_fast/tower.py

@@ -57,7 +57,7 @@ def solar_elevation(index_utc: pd.DatetimeIndex, lat: float, lon: float, step: f
     """The sun's elevation [degrees] at the middle of each record (NOAA's approximation:
     declination and equation of time from the fractional year; well within a degree)."""
     t = index_utc + pd.Timedelta(seconds=step / 2.0)
-    doy = t.dayofyear.to_numpy()
+    doy = t.day_of_year.to_numpy()
     hour = t.hour.to_numpy() + t.minute.to_numpy() / 60.0 + t.second.to_numpy() / 3600.0
     g = 2.0 * np.pi / 365.0 * (doy - 1 + (hour - 12.0) / 24.0)
     eot = 229.18 * (0.000075 + 0.001868 * np.cos(g) - 0.032077 * np.sin(g) - 0.014615 * np.cos(2 * g)