temporal_average implement seasonal intervals

* adds a seasonal interval that iterates over the seasons
  DJF, MAM, JJA, SON. Output is on a time axis centered on
  midnight of the first day of the month, of the first month
  of each season.
* includes a test, and extensions to the command line app

.travis.yml

@@ -17,7 +17,7 @@ env:
     - BUILD_TYPE=Debug
     - TECA_DIR=/travis_teca_dir
     - TECA_PYTHON_VERSION=3
-    - TECA_DATA_REVISION=98
+    - TECA_DATA_REVISION=100
   jobs:
     - DOCKER_IMAGE=ubuntu IMAGE_VERSION=20.04 IMAGE_NAME=ubuntu_20_04 REQUIRE_NETCDF_MPI=TRUE
     - DOCKER_IMAGE=ubuntu IMAGE_VERSION=20.04 IMAGE_NAME=ubuntu_20_04 REQUIRE_NETCDF_MPI=FALSE

alg/teca_temporal_reduction.py

@@ -36,6 +36,170 @@ def __str__(self):
             return strg
 
     class interval_iterator:
+        class season_iterator:
+            """
+            An iterator over seasons (DJF, MAM, JJA, SON) between 2
+            time_point's.  A pair of time steps bracketing the current season
+            are returned at each iteration. Only full seasonal intervals are
+            processed.  If the input data doesn't start or end on a seasonal
+            boundary it is skipped.
+            """
+
+            def __init__(self, t, units, calendar):
+                """
+                t - an array of floating point time values
+                units - string units of the time values
+                calendar - string name of the calendar system
+                """
+                self.t = t
+                self.units = units
+
+                calendar = calendar.lower()
+                self.calendar = calendar
+
+                # time point's to iterate between
+                self.t0 = teca_temporal_reduction_internals.time_point(
+                              t[0], units, calendar)
+
+                self.t1 = teca_temporal_reduction_internals.time_point(
+                              t[-1], units, calendar)
+
+                # current time state
+                self.year, self.month = \
+                     self.get_first_season(self.t0.year, self.t0.month)
+
+            def get_season_name(self, month):
+                """
+                returns one of DJF,MAM,JJA,SON based on the month passed in
+                """
+                if (month == 12) or ((month >= 1) and (month <= 2)):
+                    return 'DJF'
+                elif (month >= 3) and (month <= 5):
+                    return 'MAM'
+                elif (month >= 6) and (month <= 8):
+                    return 'JJA'
+                elif (month >= 9) and (month <= 11):
+                    return 'SON'
+
+                raise RuntimeError('Invalid month %d' % (month))
+
+            def get_first_season(self, y, m):
+                """
+                given a year and month, checks that the values fall on
+                a seasonal boundary. if not, returns the year and month
+                of the start of the next season.
+                """
+                if (m == 12) or (m == 3) or (m == 6) or (m == 9):
+                    return y, m
+                else:
+                    return self.get_next_season(y, m)
+
+            def get_season_end(self, year, month):
+                """
+                Given a year and month returns the year month and day
+                of the end of the season. the input month need not be on
+                a seasonal boundary.
+                """
+                if (month == 12):
+                    y = year + 1
+                    m = 2
+                elif (month >= 1) and (month <= 2):
+                    y = year
+                    m = 2
+                elif (month >= 3) and (month <= 5):
+                    y = year
+                    m = 5
+                elif (month >= 6) and (month <= 8):
+                    y = year
+                    m = 8
+                elif (month >= 9) and (month <= 11):
+                    y = year
+                    m = 11
+                else:
+                    raise RuntimeError('Invalid month %d' % (month))
+
+                d = self.last_day_of_month(y, m)
+
+                return y, m, d
+
+            def get_next_season(self, year, month):
+                """
+                Given a year and month returns the year and month
+                of the next season. the input momnth doesn't need to be
+                on a seasonal boundary.
+                """
+                if (month == 12):
+                    y = year + 1
+                    m = 3
+                elif (month >= 1) and (month <= 2):
+                    y = year
+                    m = 3
+                elif (month >= 3) and (month <= 5):
+                    y = year
+                    m = 6
+                elif (month >= 6) and (month <= 8):
+                    y = year
+                    m = 9
+                elif (month >= 9) and (month <= 11):
+                    y = year
+                    m = 12
+                else:
+                    raise RuntimeError('Invalid month %d' % (month))
+
+                return y, m
+
+            def last_day_of_month(self, year, month):
+                """
+                get the number of days in the month, with logic for
+                leap years
+                """
+                return \
+                    calendar_util.days_in_month(self.calendar,
+                                                self.units, year,
+                                                month)
+
+            def __iter__(self):
+                return self
+
+            def __next__(self):
+                """
+                return a pair of time steps bracketing the current month.
+                both returned time steps belong to the current month.
+                """
+                # get the end of the current season
+                ey, em, ed = self.get_season_end(self.year, self.month)
+
+                # verify that we have data for the current season
+                if ((ey > self.t1.year) or
+                    ((ey == self.t1.year) and (em > self.t1.month)) or
+                    ((ey == self.t1.year) and (em == self.t1.month) and
+                    (ed > self.t1.day))):
+                    raise StopIteration
+
+                # find the time step of the first day
+                sy = self.year
+                sm = self.month
+
+                t0 = '%04d-%02d-01 00:00:00' % (sy, sm)
+                i0 = coordinate_util.time_step_of(self.t, True, True,
+                                                  self.calendar,
+                                                  self.units, t0)
+
+                # find the time step of the last day
+                t1 = '%04d-%02d-%02d 23:59:59' % (ey, em, ed)
+                i1 = coordinate_util.time_step_of(self.t, True, True,
+                                                  self.calendar,
+                                                  self.units, t1)
+
+                # move to next season
+                self.year, self.month = \
+                     self.get_next_season(sy, sm)
+
+                return teca_temporal_reduction_internals.c_struct(
+                    time=self.t[i0], year=sy, month=sm,
+                    day=1, start_index=i0, end_index=i1)
+
+
         class month_iterator:
             """
             An iterator over all months between 2 time_point's. A pair
@@ -218,6 +382,11 @@ def __next__(self):
 
         @staticmethod
         def New(interval, t, units, calendar):
+            if interval == 'seasonal':
+
+                return teca_temporal_reduction_internals. \
+                    interval_iterator.season_iterator(t, units, calendar)
+
             if interval == 'monthly':
 
                 return teca_temporal_reduction_internals. \
@@ -369,7 +538,7 @@ class teca_temporal_reduction(teca_threaded_python_algorithm):
     Reduce a mesh across the time dimensions by a defined increment using
     a defined operation.
 
-        time increments: daily, monthly
+        time increments: daily, monthly, seasonal
         reduction operators: average, min, max
 
     The output time axis will be defined using the selected increment.
@@ -420,6 +589,12 @@ def set_interval(self, interval):
         """
         self.interval_name = interval
 
+    def set_interval_to_seasonal(self):
+        """
+        set the output interval to seasonal.
+        """
+        self.interval_name = 'seasonal'
+
     def set_interval_to_monthly(self):
         """
         set the output interval to monthly.
@@ -505,7 +680,7 @@ def report(self, port, md_in):
 
         t_units = t_atts['units']
 
-        # convert the time axis to a monthly delta t
+        # convert the time axis to the specified interval
         self.indices = [ii for ii in teca_temporal_reduction_internals.
                         interval_iterator.New(
                             self.interval_name, t, t_units, cal)]

apps/teca_temporal_reduction.in

@@ -26,7 +26,7 @@ parser.add_argument('--input_regex', type=str, required=True,
 
 parser.add_argument('--interval', type=str, default='monthly',
                     help='interval to reduce the time axis to. One of '
-                         'daily, or monthly (monthly)')
+                         'daily, monthly, or seasonal (monthly)')
 
 parser.add_argument('--operator', type=str, default='average',
                     help='reduction operator to use. One of minimum, '

test/apps/CMakeLists.txt

@@ -133,6 +133,22 @@ teca_add_test(test_event_filter_app
     ${CMAKE_BINARY_DIR}/${BIN_PREFIX} ${TECA_DATA_ROOT}
     REQ_TECA_DATA)
 
+teca_add_test(test_temporal_reduction_app_seasonal_average_thread
+    COMMAND ${CMAKE_CURRENT_SOURCE_DIR}/test_temporal_reduction_app.sh
+    ${CMAKE_BINARY_DIR}/${BIN_PREFIX} ${TECA_DATA_ROOT}
+    "prw_hus_day_MRI-CGCM3_historical_r1i1p1_19500101-19501231\\.nc" prw
+    seasonal average 7
+    FEATURES ${TECA_HAS_NETCDF} ${TECA_HAS_UDUNITS}
+    REQ_TECA_DATA)
+
+teca_add_test(test_temporal_reduction_app_seasonal_average_mpi_thread
+    COMMAND ${CMAKE_CURRENT_SOURCE_DIR}/test_temporal_reduction_app.sh
+    ${CMAKE_BINARY_DIR}/${BIN_PREFIX} ${TECA_DATA_ROOT}
+    "prw_hus_day_MRI-CGCM3_historical_r1i1p1_19500101-19501231\\.nc" prw
+    seasonal average 7 ${MPIEXEC} ${HALF_TEST_CORES}
+    FEATURES ${TECA_HAS_NETCDF_MPI} ${TECA_HAS_UDUNITS} ${TECA_HAS_MPI} ${MPI4Py_FOUND}
+    REQ_TECA_DATA)
+
 teca_add_test(test_temporal_reduction_app_monthly_average_thread
     COMMAND ${CMAKE_CURRENT_SOURCE_DIR}/test_temporal_reduction_app.sh
     ${CMAKE_BINARY_DIR}/${BIN_PREFIX} ${TECA_DATA_ROOT}

test/python/CMakeLists.txt

@@ -193,6 +193,21 @@ teca_add_test(py_test_python_reduce_mpi_thread
     FEATURES ${TECA_HAS_NETCDF} ${TECA_HAS_UDUNITS} ${TECA_HAS_MPI} ${MPI4Py_FOUND}
     REQ_TECA_DATA)
 
+teca_add_test(py_test_temporal_seasonal_average_thread
+    COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/test_temporal_reduction.py
+    "${TECA_DATA_ROOT}/prw_hus_day_MRI-CGCM3_historical_r1i1p1_19500101-19501231\\.nc"
+    "." "${TECA_DATA_ROOT}/test_temporal_reduction_prw" 7 2 seasonal average 0 prw
+    FEATURES ${TECA_HAS_NETCDF} ${TECA_HAS_UDUNITS}
+    REQ_TECA_DATA)
+
+teca_add_test(py_test_temporal_seasonal_average_mpi_thread
+    COMMAND ${MPIEXEC} -n ${HALF_TEST_CORES} ${PYTHON_EXECUTABLE}
+    ${CMAKE_CURRENT_SOURCE_DIR}/test_temporal_reduction.py
+    "${TECA_DATA_ROOT}/prw_hus_day_MRI-CGCM3_historical_r1i1p1_19500101-19501231\\.nc"
+    "." "${TECA_DATA_ROOT}/test_temporal_reduction_prw" 7 2 seasonal average 0 prw
+    FEATURES ${TECA_HAS_NETCDF} ${TECA_HAS_UDUNITS} ${TECA_HAS_MPI} ${MPI4Py_FOUND}
+    REQ_TECA_DATA)
+
 teca_add_test(py_test_temporal_monthly_average_thread
     COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/test_temporal_reduction.py
     "${TECA_DATA_ROOT}/prw_hus_day_MRI-CGCM3_historical_r1i1p1_19500101-19501231\\.nc"

test/python/test_temporal_reduction.py

@@ -49,7 +49,7 @@
 mav.set_operator(operator)
 mav.set_point_arrays(arrays)
 mav.set_use_fill_value(use_fill)
-mav.set_verbose(1)
+mav.set_verbose(2)
 mav.set_thread_pool_size(n_threads)
 mav.set_stream_size(2)
 
@@ -68,6 +68,7 @@
     diff = teca_dataset_diff.New()
     diff.set_input_connection(0, bcfr.get_output_port())
     diff.set_input_connection(1, mav.get_output_port())
+    diff.set_tolerance(1e-5)
     diff.set_executive(exe)
 
     diff.update()