Merge pull request #662 from LinkedEarth/annualize

Annualize

Author: jordanplanders

condaenv.6d4f03jk.requirements.txt

@@ -126,7 +126,7 @@ class Series:
          set to True to remove the NaNs and make time axis strictly prograde with duplicated timestamps reduced by averaging the values
          Default is None (marked for deprecation)
 
-    auto_time_params : bool,
+    auto_time_params : bool
         If True, uses tsbase.disambiguate_time_metadata to ensure that time_name and time_unit are usable by Pyleoclim. This may override the provided metadata.
         If False, the provided time_name and time_unit are used. This may break some functionalities (e.g. common_time and convert_time_unit), so use at your own risk.
         If not provided, code will set to True for internal consistency.
@@ -2399,6 +2399,250 @@ def sort(self, verbose=False, ascending = True, keep_log = False):
                 new.log=()
             new.log += ({len(new.log):'sort', 'ascending': ascending},)
         return new
+    
+    def annualize(self, months=list(range(1, 13)), min_res=0.25, frac_req_months=2/3):
+        '''
+        Annualize subannual data by averaging values within specified months for each year.
+        
+        This method converts subannual time series data to annual resolution by computing
+        weighted averages of values from specified months using the custom_year_averages
+        utility function.
+        
+        Parameters
+        ----------
+        months : list of int, optional
+            List of months to include in the annual aggregation. Must be consecutive
+            months either within a calendar year or spanning across years (e.g., 
+            [10, 11, 12, 1, 2, 3] for Oct-Mar). Default is a calendar average, [1 ... 12].
+        min_res : float, optional
+            Minimum resolution threshold (in years) to determine if data is subannual.
+            Default is 0.25 (quarterly or finer resolution required).
+        frac_req_months : float, optional
+            Minimum fraction of requested months that must be present for a year to be included.
+            Default is 2/3 (67%). For example, if requesting 12 months, at least 8 months must be 
+            present. If requesting 3 months (e.g., DJF), at least 2 months must be present.
+            Years that do not meet this threshold are dropped from the resulting series.
+            Set to 1.0 to require all months, or lower values (e.g., 1/3) for more lenient inclusion.
+        
+        Returns
+        -------
+        Series
+            A new Series object with annualized data.
+        
+        Examples
+        --------
+        1) Annual average
+        .. jupyter-execute::
+ 
+            soi = pyleo.utils.load_dataset('SOI')
+            dt = soi.resolution().describe()['median']
+            print(f"The series' resolution is {dt*12:.0f} month")
+
+            soi_a = soi.annualize()
+            dta = soi_a.resolution().describe()['median']
+            print(f"The series' resolution is now {dta:.0f} year")
+            fig, ax = soi.plot(title='Jan-Dec averaging')
+            soi_a.plot(marker='o',ax=ax)
+
+        2) JJA average : straightforward
+        .. jupyter-execute::
+            
+            soi_jja = soi.annualize(months=[6, 7 , 8])
+            fig, ax = soi.plot(title='JJA averaging')
+            soi_jja.plot(marker='o',ax=ax, label='JJA average')
+
+        3) DJF average : straddles a year; handles it gracefully
+        .. jupyter-execute::
+            
+            soi_djf = soi.annualize(months=[12, 1 , 2])
+            fig, ax = soi.plot(title='DJF averaging')
+            soi_djf.plot(marker='o',ax=ax, label='DJF average')
+
+        4) Varying the fraction of required months
+        .. jupyter-execute::
+            
+            AprMar = [4,5,6,7,8,9,10,11,12,1,2,3]
+            soi_am_default = soi.annualize(months=AprMar)
+            soi_am_stringent = soi.annualize(months=AprMar,frac_req_months=0.9)
+
+            fig, ax = soi.plot(title='Apr-Mar averaging')
+            soi_am_default.plot(marker='o',ax=ax, label='Apr-Mar, $f=2/3$')
+            soi_am_stringent.plot(marker='o',ax=ax, label='Apr-Mar, $f=0.9$')
+        
+        We see that insisting on a very high fraction of available months will result in dropped years 
+        
+        Notes
+        -----
+        - This method requires subannual data (median resolution <= min_res)
+        - Months must be consecutive to define a clear averaging period
+        - Uses weighted averaging to account for potentially uneven temporal spacing
+        '''
+        
+        # Check if data is subannual using Resolution class
+        if not hasattr(self, 'time_unit') or 'year' not in self.time_unit.lower():
+            raise ValueError("time_unit must contain 'Year' to attempt annualization")
+        
+        if not 0 <= min_res <= 0.5: 
+            raise ValueError(f"min_res must be in (0, 0.5); got {min_res}")
+        
+        # Check for suitable time axis properties
+        dt = self.resolution().describe()['median']
+        if dt >= min_res:
+            raise ValueError(f"Data appears to be too coarse (median resolution: {dt:.1f} years >= {min_res}) for meaningful averages.")
+        if dt < 0:
+            warnings.warn("Series was retrograde; it has been sorted in ascending order so annualization runs properly", RuntimeWarning)
+            self = self.sort()
+        
+        # Set default months if not provided
+        if months is None:
+            months = list(range(1, 13))
+        
+        # Validate months parameter
+        if not isinstance(months, (list, tuple)):
+            raise ValueError("months must be a list or tuple of integers")
+        
+        if not all(isinstance(m, int) and 1 <= m <= 12 for m in months):
+            raise ValueError("months must contain integers between 1 and 12")
+        
+        if len(months) == 0:
+            raise ValueError("months cannot be empty")
+        
+        # Check if months are consecutive (allowing for year-end wrapping)
+        def is_consecutive_with_wrap(months_list):
+            """Check if months are consecutive, allowing wrapping around year boundary"""
+            if len(months_list) <= 1:
+                return True
+            
+            # Check normal consecutive case
+            is_normal_consecutive = True
+            for i in range(1, len(months_list)):
+                if months_list[i] != months_list[i-1] + 1:
+                    is_normal_consecutive = False
+                    break
+            
+            if is_normal_consecutive:
+                return True
+            
+            # Check wrap-around case (e.g., [11, 12, 1, 2] or [12, 1, 2])
+            # Find where the wrap occurs (where month decreases)
+            wrap_point = None
+            for i in range(1, len(months_list)):
+                if months_list[i] < months_list[i-1]:
+                    if wrap_point is None:
+                        wrap_point = i
+                    else:
+                        return False  # Multiple wraps not allowed
+            
+            if wrap_point is None:
+                return False  # No wrap found but not consecutive
+            
+            # Check consecutive before wrap
+            for i in range(1, wrap_point):
+                if months_list[i] != months_list[i-1] + 1:
+                    return False
+            
+            # Check consecutive after wrap
+            for i in range(wrap_point + 1, len(months_list)):
+                if months_list[i] != months_list[i-1] + 1:
+                    return False
+            
+            # Check that the wrap makes sense (December -> January)
+            if months_list[wrap_point-1] == 12 and months_list[wrap_point] == 1:
+                return True
+            
+            return False
+        
+        if not is_consecutive_with_wrap(months):
+            raise ValueError("months must be consecutive (e.g., [1,2,3] or [11,12,1,2])")
+        
+        # Determine start and end months (use original order, not sorted)
+        start_month = months[0]
+        end_month = months[-1]
+        
+        # Convert to pandas Series using existing method
+        try:
+            series = self.to_pandas()
+        except Exception as e:
+            raise ValueError(f"Could not convert Series to pandas: {e}")
+        
+        # Ensure the index is datetime
+        if not hasattr(series.index, 'month'):
+            try:
+                series.index = pd.to_datetime(series.index)
+            except Exception as e:
+                raise ValueError(f"Could not convert time axis to datetime: {e}")
+        
+        # Use the new custom_year_averages function
+        try:
+            annual_data = tsutils.custom_year_averages(
+                data=series,
+                start_month=start_month,
+                end_month=end_month,
+                years=None  # Use all available years
+            )
+        except Exception as e:
+            raise ValueError(f"Failed to compute annual averages: {e}")
+        
+        # Apply frac_req_months filtering
+        # Filter out years that don't meet the minimum fraction requirement
+        sufficient_years = []
+        
+        for year in annual_data.index:
+            # Determine the date range for this year's period
+            if start_month <= end_month:
+                # Same calendar year
+                period_start = pd.Timestamp(year=year, month=start_month, day=1)
+                period_end = pd.Timestamp(year=year, month=end_month, day=28)  # Use 28 to be safe
+            else:
+                # Spans calendar years
+                period_start = pd.Timestamp(year=year-1, month=start_month, day=1)
+                period_end = pd.Timestamp(year=year, month=end_month, day=28)  # Use 28 to be safe
+            
+            # Get data for this period
+            period_mask = (series.index >= period_start) & (series.index <= period_end)
+            period_data = series[period_mask]
+            
+            if len(period_data) == 0:
+                continue
+            
+            # Count available months (works for any temporal resolution, not just monthly)
+            available_months = set(period_data.index.month)
+            required_months = set(months)
+            available_required_months = available_months.intersection(required_months)
+            
+            # Calculate fraction of required months that are present
+            fraction_present = len(available_required_months) / len(required_months)
+            
+            if fraction_present >= frac_req_months:
+                sufficient_years.append(year)
+        
+        # Filter annual_data to only include years with sufficient data
+        if sufficient_years:
+            annual_data = annual_data[annual_data.index.isin(sufficient_years)]
+        else:
+            annual_data = pd.Series([], dtype=float, name=annual_data.name)
+        
+        # Report dropped years
+        total_years_before = len(tsutils.custom_year_averages(series, start_month, end_month, years=None).index)
+        dropped_years = total_years_before - len(annual_data)
+        
+        if dropped_years > 0:
+            warnings.warn(f"Dropped {dropped_years} of {total_years_before} years due to insufficient data coverage (< {frac_req_months:.1%} of requested months)", RuntimeWarning)
+        
+        if annual_data.empty:
+            raise ValueError(f"No years found with sufficient data coverage (>= {frac_req_months:.1%} of requested months). Consider lowering frac_req_months.")
+        
+        # Create new Series by copying self and replacing time/value
+        annual_series = self.copy()
+        annual_series.time = annual_data.index.values.astype(float)
+        annual_series.value = annual_data.values
+        
+        # Update label to indicate annualization
+        if self.label:
+            annual_series.label = f"{self.label} (annualized)"
+        
+        return annual_series
+    
 
     def gaussianize(self, keep_log = False):
         ''' Gaussianizes the timeseries (i.e. maps its values to a standard normal)