Repository: metno/pyaerocom

Issue:
season definition for evaluations of type "long" in CAMS2_83
for the evaluations of type "long" in cams2_83 the periods are years, plus the full range, e.g. `"periods": ["2022-2025", "2022", "2023", "2024", "2025"]`  
and this results in these periods in the output map files, since we use (so far) `add_seasons=True` :
```
2022-2025-all
2022-2025-DJF
2022-2025-MAM
2022-2025-JJA
2022-2025-SON
2022-all
2022-DJF
2022-MAM
2022-JJA
2022-SON
```
the issue is that DJF is not defined in the meteorological sense, so for example the key `2022-DJF` will be populated also for a run that starts at March 1st 2022, and will contain statistics calculated on the base of Dec 2022, which is wrong, similarly `2023-DJF` will contain results on the base of Jan-Feb 2022 and Dec 2023  

all this comes from the fact that the season split for cams2_83 is not differentiated from the rest of pyaerocom:  
the loop on seasons is a loop on this set of strings: 'all', 'DJF', 'MAM', 'JJA', 'SON'     
for a given year period, say 2022, the function `_select_period_season_coldata` https://github.com/metno/pyaerocom/blob/temporary-cams2-83/pyaerocom/aeroval/coldatatojson_helpers.py#L837, which is called by `_proc_map_and_scat` for example, checks if a season, say 'DJF', is in the seasons array in the colocated data [see below for an example] after the year has been subset from the full time data, without distinguishing between the part of DJF that belongs to the year in question or the next year, as a meteorological seasons...   

here for example a `subset` object returned by this call https://github.com/metno/pyaerocom/blob/temporary-cams2-83/pyaerocom/aeroval/coldatatojson_helpers.py#L837 when the period is '2022' and the season is 'DJF', for a 'long' run covering 20220101 20230228:
```
(Pdb) per
'2022'
(Pdb) season
'DJF'
(Pdb) subset.coords['time'].values
array(['2022-01-01T12:00:00.000000000', '2022-01-02T12:00:00.000000000',
       '2022-01-03T12:00:00.000000000', '2022-01-04T12:00:00.000000000',
       '2022-01-05T12:00:00.000000000', '2022-01-06T12:00:00.000000000',
       '2022-01-07T12:00:00.000000000', '2022-01-08T12:00:00.000000000',
       '2022-01-09T12:00:00.000000000', '2022-01-10T12:00:00.000000000',
       '2022-01-11T12:00:00.000000000', '2022-01-12T12:00:00.000000000',
       '2022-01-13T12:00:00.000000000', '2022-01-14T12:00:00.000000000',
       '2022-01-15T12:00:00.000000000', '2022-01-16T12:00:00.000000000',
       '2022-01-17T12:00:00.000000000', '2022-01-18T12:00:00.000000000',
       '2022-01-19T12:00:00.000000000', '2022-01-20T12:00:00.000000000',
       '2022-01-21T12:00:00.000000000', '2022-01-22T12:00:00.000000000',
       '2022-01-23T12:00:00.000000000', '2022-01-24T12:00:00.000000000',
       '2022-01-25T12:00:00.000000000', '2022-01-26T12:00:00.000000000',
       '2022-01-27T12:00:00.000000000', '2022-01-28T12:00:00.000000000',
       '2022-01-29T12:00:00.000000000', '2022-01-30T12:00:00.000000000',
       '2022-01-31T12:00:00.000000000', '2022-02-01T12:00:00.000000000',
       '2022-02-02T12:00:00.000000000', '2022-02-03T12:00:00.000000000',
       '2022-02-04T12:00:00.000000000', '2022-02-05T12:00:00.000000000',
       '2022-02-06T12:00:00.000000000', '2022-02-07T12:00:00.000000000',
       '2022-02-08T12:00:00.000000000', '2022-02-09T12:00:00.000000000',
       '2022-02-10T12:00:00.000000000', '2022-02-11T12:00:00.000000000',
       '2022-02-12T12:00:00.000000000', '2022-02-13T12:00:00.000000000',
       '2022-02-14T12:00:00.000000000', '2022-02-15T12:00:00.000000000',
       '2022-02-16T12:00:00.000000000', '2022-02-17T12:00:00.000000000',
       '2022-02-18T12:00:00.000000000', '2022-02-19T12:00:00.000000000',
       '2022-02-20T12:00:00.000000000', '2022-02-21T12:00:00.000000000',
       '2022-02-22T12:00:00.000000000', '2022-02-23T12:00:00.000000000',
       '2022-02-24T12:00:00.000000000', '2022-02-25T12:00:00.000000000',
       '2022-02-26T12:00:00.000000000', '2022-02-27T12:00:00.000000000',
       '2022-02-28T12:00:00.000000000', '2022-12-01T12:00:00.000000000',
       '2022-12-02T12:00:00.000000000', '2022-12-03T12:00:00.000000000',
       '2022-12-04T12:00:00.000000000', '2022-12-05T12:00:00.000000000',
       '2022-12-06T12:00:00.000000000', '2022-12-07T12:00:00.000000000',
       '2022-12-08T12:00:00.000000000', '2022-12-09T12:00:00.000000000',
       '2022-12-10T12:00:00.000000000', '2022-12-11T12:00:00.000000000',
       '2022-12-12T12:00:00.000000000', '2022-12-13T12:00:00.000000000',
       '2022-12-14T12:00:00.000000000', '2022-12-15T12:00:00.000000000',
       '2022-12-16T12:00:00.000000000', '2022-12-17T12:00:00.000000000',
       '2022-12-18T12:00:00.000000000', '2022-12-19T12:00:00.000000000',
       '2022-12-20T12:00:00.000000000', '2022-12-21T12:00:00.000000000',
       '2022-12-22T12:00:00.000000000', '2022-12-23T12:00:00.000000000',
       '2022-12-24T12:00:00.000000000', '2022-12-25T12:00:00.000000000',
       '2022-12-26T12:00:00.000000000', '2022-12-27T12:00:00.000000000',
       '2022-12-28T12:00:00.000000000', '2022-12-29T12:00:00.000000000',
       '2022-12-30T12:00:00.000000000', '2022-12-31T12:00:00.000000000'],
      dtype='datetime64[ns]')
```

[example of colocated data] the colocated data object itself has seasons as a coordinate where only the 3 letter string is recorded 
 ```
 (Pdb) cd.coords['season'].values
 array(['DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF',
       'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF',
       'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF',
       'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF',
       'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF',
       'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF',
       'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'MAM', 'MAM', 'MAM', 'MAM',
       'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM',
       'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM',
       'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM',
       'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM',
       'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM',
       'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM',
       'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM',
       'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM',
       'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM',
       'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'JJA', 'JJA',
       'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA',
       'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA',
       'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA',
       'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA',
       'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA',
       'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA',
       'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA',
       'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA',
       'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA',
       'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA', 'JJA',
       'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON',
       'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON',
       'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON',
       'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON',
       'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON',
       'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON',
       'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON',
       'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON',
       'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON',
       'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON', 'SON',
       'SON', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF',
       'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF',
       'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF',
       'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF',
       'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF',
       'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF',
       'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF',
       'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF',
       'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF',
       'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF',
       'DJF', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM',
       'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM'], dtype='<U3')
```
so slicing by year will keep all these within that year regardless of the meteorological sense

Discussion:
@charlienegri: basically we cannot rely on pyaerocom's season definition,  but circumventing that implies also that the periods have probably to include the seasons directly, which means there will be no more keys like 'MAM-2023' in the map files...

to clarify:`"periods": ["2022-2025", "2022", "2023", "2024", "2025"]` and `add_seasons=True`  
results in the keys  
```
2022-2025-all
2022-2025-DJF
2022-2025-MAM
2022-2025-JJA
2022-2025-SON
2022-all
2022-DJF
2022-MAM
2022-JJA
2022-SON
```

while having periods like ["20211201-20220228","20220301-20220531",.., "20221201-20230228"] and `add_seasons=False`  
would result in keys like: 
```
2021/12/01-2022/02/28-all
2022/03/01-2022/05/31-all
...
2022/12/01-2023/02/28-all
```
the menus on aeroval.met.no would not have a [all|DJF|MAM|JJA|SON] selection anymore, each season would be listed as a single date range in the drop down, like "2022/03/01-2022/05/31"   
the frontend would not complain, but it would not be possible to show the results as they are now directly tho

another option is to try to fix the bug in the ` _select_period_season_coldata` function if we wanna keep the current structure... but differentiating the behavior for cams2_83 only is not straightforward@charlienegri: the first option is implemented in a branch and gives something like this https://aeroval-test.met.no/charlien/pages/evaluation/?project=cams2-83&experiment=test_long_seasonasperiods_2 (only EMEP model, no median scores)  
if the different look on the website is acceptable, I can open a PR@charlienegri: after discussing w Daniel, we think we can actually fix the behavior of  `select_period_season_coldata` on-demand (to be discussed if it should be an option or what) so that DJF is handled correctly
@heikoklein: @charlienegri Did #1540 fix this issue? Can we thus close it?@charlienegri: yeah

Solution patch:
diff --git a/docs/_static/aeroval/cfg_examples_example1.py b/docs/_static/aeroval/cfg_examples_example1.py
index 19ff104a..7893aa46 100644
--- a/docs/_static/aeroval/cfg_examples_example1.py
+++ b/docs/_static/aeroval/cfg_examples_example1.py
@@ -147,8 +147,19 @@ GLOB_CFG = dict(
     main_freq="monthly",
     # Time periods for which statistical parameters are computed
     periods=["2010"],
-    # Whether or not to add seasonal statistics
+    # Whether or not to add seasonal statistics.
+    # If True seasons are ['all', 'DJF', 'MAM', 'JJA', 'SON'],
+    # if False, just ['all']
     add_seasons=True,
+    # Wheter or not statistics are based on the meteorological definition of seasons.
+    # This is relevant for periods that are a single year. So if use_meteorological_seasons=True
+    # and add_seasons=True, for a given year ['DJF'] will refer to data from Dec of the previous year
+    # (if available) and Jan/Feb of the same year, while if use_meteorological_seasons=False,
+    # it will refer to data from Jan/Feb and December of the same year. Similarly, and independently
+    # of the value of add_seasons, if use_meteorological_seasons=True, ['all'] (whole year) will refer
+    # to data from Dec of the previous year to Nov of the same year, while if False, it will refer to data
+    # from Jan to Dec of the same year.
+    use_meteorological_seasons=False,
     # Whether or not to add trends output to the analysis. Trends analysis
     # needs at least 7 years of data, so this is skipped here for this
     # single year experiment
diff --git a/pyaerocom/aeroval/coldatatojson_engine.py b/pyaerocom/aeroval/coldatatojson_engine.py
index f03bf01f..aa405051 100644
--- a/pyaerocom/aeroval/coldatatojson_engine.py
+++ b/pyaerocom/aeroval/coldatatojson_engine.py
@@ -86,6 +86,7 @@ class ColdataToJsonEngine(ProcessingEngine):
         freqs = self.cfg.time_cfg.freqs
         periods = self.cfg.time_cfg.periods
         seasons = self.cfg.time_cfg.get_seasons()
+        use_meteorological_seasons = self.cfg.time_cfg.use_meteorological_seasons
         main_freq = self.cfg.time_cfg.main_freq
         annual_stats_constrained = self.cfg.statistics_opts.annual_stats_constrained
 
@@ -228,6 +229,7 @@ class ColdataToJsonEngine(ProcessingEngine):
                     stats_min_num=stats_min_num,
                     use_fairmode=use_fairmode,
                     avg_over_trends=avg_over_trends,
+                    use_meteorological_seasons=use_meteorological_seasons,
                 )
             if coldata.ts_type == "hourly" and use_diurnal:
                 logger.info("Processing diurnal profiles")
@@ -249,6 +251,7 @@ class ColdataToJsonEngine(ProcessingEngine):
                 seasons=seasons,
                 obs_name=obs_name,
                 var_name_web=var_name_web,
+                use_meteorological_seasons=use_meteorological_seasons,
             )
 
         logger.info(
@@ -293,6 +296,7 @@ class ColdataToJsonEngine(ProcessingEngine):
         seasons: tuple[str, ...] = None,
         obs_name: str = None,
         var_name_web: str = None,
+        use_meteorological_seasons: bool = False,
     ):
         if region_names is None and station_names is None:
             raise ValueError("Both region_id and station_name can not both be None")
@@ -305,6 +309,7 @@ class ColdataToJsonEngine(ProcessingEngine):
                 use_country=use_country,
                 periods=periods,
                 seasons=seasons,
+                use_meteorological_seasons=use_meteorological_seasons,
             )
             location = region_names[regid]
             self.exp_output.add_profile_entry(
@@ -325,6 +330,7 @@ class ColdataToJsonEngine(ProcessingEngine):
                 use_country=use_country,
                 periods=periods,
                 seasons=seasons,
+                use_meteorological_seasons=use_meteorological_seasons,
             )
 
             self.exp_output.add_profile_entry(
@@ -363,6 +369,7 @@ class ColdataToJsonEngine(ProcessingEngine):
         stats_min_num: int = 1,
         use_fairmode: bool = False,
         avg_over_trends: bool = False,
+        use_meteorological_seasons: bool = False,
     ):
         input_freq = self.cfg.statistics_opts.stats_tseries_base_freq
 
@@ -414,6 +421,7 @@ class ColdataToJsonEngine(ProcessingEngine):
             add_trends,
             trends_min_yrs,
             avg_over_trends,
+            use_meteorological_seasons,
         )
 
         for freq, hm_data in hm_all.items():
@@ -456,6 +464,7 @@ class ColdataToJsonEngine(ProcessingEngine):
                 use_fairmode,
                 obs_var,
                 drop_stats,
+                use_meteorological_seasons,
             )
 
             with self.avdb.lock():
diff --git a/pyaerocom/aeroval/coldatatojson_helpers.py b/pyaerocom/aeroval/coldatatojson_helpers.py
index 42cc49dc..e759ea9b 100644
--- a/pyaerocom/aeroval/coldatatojson_helpers.py
+++ b/pyaerocom/aeroval/coldatatojson_helpers.py
@@ -766,7 +766,7 @@ def process_trends(
     freq: str,
     use_weights: bool,
 ) -> dict:
-    # subset = _select_period_season_coldata(coldata, per, season)
+    # subset = _select_period_season_coldata(coldata, per, season, use_meteorological_seasons)
     stats = {}
     trends_successful = False
     mean_trends_successful = False
@@ -976,6 +976,7 @@ def _process_map_and_scat(
     use_fairmode,
     obs_var,
     drop_stats,
+    use_meteorological_seasons,
 ):
     stats_dummy = _init_stats_dummy(drop_stats=drop_stats)
     scat_data = {}
@@ -986,7 +987,9 @@ def _process_map_and_scat(
                 use_dummy = cd is None
                 if not use_dummy:
                     try:
-                        subset = _select_period_season_coldata(cd, per, season)
+                        subset = _select_period_season_coldata(
+                            cd, per, season, use_meteorological_seasons
+                        )
                         jsdate = subset.data.jsdate.values.tolist()
                     except (DataCoverageError, TemporalResolutionError):
                         use_dummy = True
@@ -1253,8 +1256,14 @@ def _calc_temporal_corr(coldata):
         return (np.nanmean(corr_time.data), np.nanmedian(corr_time.data))
 
 
-def _select_period_season_coldata(coldata, period, season):
+def _select_period_season_coldata(coldata, period, season, use_meteorological_seasons):
     tslice = _period_str_to_timeslice(period)
+    if use_meteorological_seasons and len(period) == 4:
+        # relevant only for single years
+        # for period = '2022' tslice needs to be slice('2021-12','2022-11')
+        yeardt = datetime.strptime(period, "%Y")
+        tslice = slice(f"{yeardt.year - 1}-12", f"{yeardt.year}-11")
+        logger.debug(f"Using meteorological year slicing for {period}: {tslice}")
     # expensive, try use solution with numpy indexing directly...
     # also, keep an eye on: https://github.com/pydata/xarray/issues/2799
     arr = coldata.data.sel(time=tslice)
@@ -1285,6 +1294,7 @@ def _process_heatmap_data(
     add_trends,
     trends_min_yrs,
     avg_over_trends,
+    use_meteorological_seasons,
 ):
     output = {}
     stats_dummy = _init_stats_dummy(drop_stats=drop_stats)
@@ -1300,7 +1310,9 @@ def _process_heatmap_data(
                         stats = stats_dummy
                     else:
                         try:
-                            subset = _select_period_season_coldata(coldata, per, season)
+                            subset = _select_period_season_coldata(
+                                coldata, per, season, use_meteorological_seasons
+                            )
 
                             if add_trends and freq != "daily":
                                 trend_stats = process_trends(
@@ -1423,8 +1435,7 @@ def _process_statistics_timeseries(
     # input frequency is lower resolution than output frequency
     if TsType(data_freq) < TsType(freq):
         raise TemporalResolutionError(
-            f"Desired input frequency {data_freq} is lower than desired "
-            f"output frequency {freq}"
+            f"Desired input frequency {data_freq} is lower than desired output frequency {freq}"
         )
 
     output = {}
@@ -1520,6 +1531,7 @@ def process_profile_data_for_regions(
     use_country: bool,
     periods: list[str],
     seasons: list[str],
+    use_meteorological_seasons: bool,
 ) -> dict:  # pragma: no cover
     """
     This method populates the json files in data/profiles which are use for visualization.
@@ -1558,7 +1570,9 @@ def process_profile_data_for_regions(
                     output["mod"][freq][perstr] = np.nan
                 else:
                     try:
-                        per_season_subset = _select_period_season_coldata(coldata, per, season)
+                        per_season_subset = _select_period_season_coldata(
+                            coldata, per, season, use_meteorological_seasons
+                        )
 
                         subset = per_season_subset.filter_region(
                             region_id=region_id, check_country_meta=use_country
@@ -1587,6 +1601,7 @@ def process_profile_data_for_stations(
     use_country: bool,
     periods: list[str],
     seasons: list[str],
+    use_meteorological_seasons: bool,
 ) -> dict:  # pragma: no cover
     """
     This method populates the json files in data/profiles which are use for visualization.
@@ -1625,7 +1640,9 @@ def process_profile_data_for_stations(
                     output["mod"][freq][perstr] = np.nan
                 else:
                     try:
-                        per_season_subset = _select_period_season_coldata(coldata, per, season)
+                        per_season_subset = _select_period_season_coldata(
+                            coldata, per, season, use_meteorological_seasons
+                        )
 
                         subset = per_season_subset.data[
                             :,
@@ -1672,7 +1689,7 @@ def _remove_less_covered(
 
     new_stations = data.data.station_name.data
 
-    logger.info(f"Removed {len(stations)-len(new_stations)} stations")
+    logger.info(f"Removed {len(stations) - len(new_stations)} stations")
     if len(new_stations) == 0:
         logger.warning(
             f"No stations left after removing stations with fewer than {min_yrs} years!"
diff --git a/pyaerocom/aeroval/setup_classes.py b/pyaerocom/aeroval/setup_classes.py
index 28e77f0a..5cb59eed 100644
--- a/pyaerocom/aeroval/setup_classes.py
+++ b/pyaerocom/aeroval/setup_classes.py
@@ -1,3 +1,4 @@
+import datetime
 import logging
 import os
 import sys
@@ -6,7 +7,6 @@ from functools import cached_property
 from getpass import getuser
 from pathlib import Path
 from typing import Annotated, Literal
-import datetime
 
 from pyaerocom.aeroval.glob_defaults import VarWebInfo, VarWebScaleAndColormap
 from pyaerocom.aeroval.obsentry import ObsEntry
@@ -22,8 +22,8 @@ from pydantic import (
     BaseModel,
     ConfigDict,
     Field,
-    PositiveInt,
     NonNegativeInt,
+    PositiveInt,
     computed_field,
     field_serializer,
     field_validator,
@@ -35,13 +35,13 @@ from pyaerocom.aeroval.aux_io_helpers import ReadAuxHandler
 from pyaerocom.aeroval.collections import ModelCollection, ObsCollection
 from pyaerocom.aeroval.exceptions import ConfigError
 from pyaerocom.aeroval.helpers import (
+    BoundingBox,
     _check_statistics_periods,
     _get_min_max_year_periods,
     check_if_year,
-    BoundingBox,
 )
-from pyaerocom.aeroval.modelmaps_helpers import CONTOUR, OVERLAY
 from pyaerocom.aeroval.json_utils import read_json, set_float_serialization_precision
+from pyaerocom.aeroval.modelmaps_helpers import CONTOUR, OVERLAY
 from pyaerocom.colocation.colocation_setup import ColocationSetup
 
 logger = logging.getLogger(__name__)
@@ -246,12 +246,31 @@ class StatisticsSetup(BaseModel, extra="allow"):
 
 
 class TimeSetup(BaseModel):
+    """
+    Time setup options
+
+    Attributes
+    ----------
+    add_seasons : bool, default True
+        if True, seasons will be ['all', 'DJF', 'MAM', 'JJA', 'SON'], if False, just ['all'].
+    use_meteorological_seasons : bool, default False
+        if True, then statistics are based on the meteorological definition of seasons. This is relevant
+        for periods that are a single year. So if :attr:`add_seasons` is True, for a given year ['DJF'] will
+        refer to data from Dec of the previous year (if available) and Jan/Feb of the same year, while if
+        :attr:`use_meteorological_seasons` is False, it will be based on data from Jan/Feb and December
+        of the same year. Similarly, and weather or not :attr:`add_seasons` is True,
+        if :attr:`use_meteorological_seasons` is True, ['all'] (whole year) will refer to data from Dec of
+        the previous year to Nov of the same year, while if False, it will refer to data from Jan to Dec
+        of the same year.
+    """
+
     DEFAULT_FREQS: Literal["monthly", "yearly"] = "monthly"
     SEASONS: list[str] = ["all", "DJF", "MAM", "JJA", "SON"]
     main_freq: str = "monthly"
     freqs: list[str] = ["monthly", "yearly"]
     periods: list[str] = Field(default_factory=list)
     add_seasons: bool = True
+    use_meteorological_seasons: bool = False
 
     def get_seasons(self):
         """
diff --git a/pyaerocom/scripts/cams2_83/cli_mos.py b/pyaerocom/scripts/cams2_83/cli_mos.py
index f141047c..68a09d7d 100644
--- a/pyaerocom/scripts/cams2_83/cli_mos.py
+++ b/pyaerocom/scripts/cams2_83/cli_mos.py
@@ -26,6 +26,7 @@ def make_config_mos(
     id: str,
     name: str,
     description: str,
+    add_seasons: bool,
 ) -> dict:
     logger.info("Making the configuration")
 
@@ -57,6 +58,9 @@ def make_config_mos(
 
     cfg.update(use_fairmode=True)
 
+    if add_seasons:
+        cfg.update(add_seasons=True)
+
     return cfg
 
 
@@ -90,6 +94,7 @@ def main(
     id: str = typer.Option(CFG["exp_id"], help="experiment ID"),
     name: str = typer.Option(CFG["exp_name"], help="experiment name"),
     description: str = typer.Option(CFG["exp_descr"], help="experiment description"),
+    add_seasons: bool = typer.Option(False, "--addseasons", help="set add_seasons"),
     pool: int = typer.Option(
         1,
         "--pool",
@@ -113,6 +118,7 @@ def main(
         id,
         name,
         description,
+        add_seasons,
     )
 
     # we do not want the cache produced in previous runs to be silently cleared
diff --git a/pyaerocom/scripts/cams2_83/config.py b/pyaerocom/scripts/cams2_83/config.py
index 496f1e19..c789289e 100644
--- a/pyaerocom/scripts/cams2_83/config.py
+++ b/pyaerocom/scripts/cams2_83/config.py
@@ -44,6 +44,7 @@ GLOBAL_CONFIG = dict(
     ],  # Periodes, can be single years or range, e.g. 2010-2015. EMEP only supports single years as of now
     main_freq="hourly",  # default frequency to use. This will be overwritten in most of the observation options (see below)
     add_seasons=False,
+    use_meteorological_seasons=True,
     # This has to be true for the web interface to show diurnal evaluation
     use_diurnal=False,
     # O3 is special, since we want to look at daily max
diff --git a/pyaerocom/scripts/cams2_83/engine.py b/pyaerocom/scripts/cams2_83/engine.py
index e396a662..ae45046a 100644
--- a/pyaerocom/scripts/cams2_83/engine.py
+++ b/pyaerocom/scripts/cams2_83/engine.py
@@ -39,6 +39,7 @@ class CAMS2_83_Engine(ProcessingEngine):
         out_dirs = self.cfg.path_manager.get_json_output_dirs(True)
         forecast_days = self.cfg.statistics_opts.forecast_days
         periods = self.cfg.time_cfg.periods
+        use_meteorological_seasons = self.cfg.time_cfg.use_meteorological_seasons
 
         if "var_name_input" in coldata[0].metadata:
             obs_var = coldata[0].metadata["var_name_input"][0]
@@ -97,7 +98,7 @@ class CAMS2_83_Engine(ProcessingEngine):
 
                     try:
                         subset = [
-                            _select_period_season_coldata(col, per, season)
+                            _select_period_season_coldata(col, per, season, use_meteorological_seasons)
                             for col in subset_region
                         ]
                     except (DataCoverageError, UnknownRegion) as e:
diff --git a/pyaerocom/ungriddeddata.py b/pyaerocom/ungriddeddata.py
index 1f558da6..ee0d9e00 100644
--- a/pyaerocom/ungriddeddata.py
+++ b/pyaerocom/ungriddeddata.py
@@ -1303,7 +1303,9 @@ class UngriddedData:
                 if hasattr(data, "station_type"):
                     out_data["station_type"].append(data["station_type"])
                 else:
-                    logger.info("No station_type found in StationData, station_type will be blank")
+                    logger.debug(
+                        "No station_type found in StationData, station_type will be blank"
+                    )
                 out_data["stats"].append(data)
 
             # catch the exceptions that are acceptable


Test patch:
diff --git a/tests/aeroval/test_coldatatojson_helpers2.py b/tests/aeroval/test_coldatatojson_helpers2.py
index ae24cb89..07337782 100644
--- a/tests/aeroval/test_coldatatojson_helpers2.py
+++ b/tests/aeroval/test_coldatatojson_helpers2.py
@@ -16,6 +16,7 @@ from pyaerocom.aeroval.coldatatojson_helpers import (
     _map_indices,
     _process_statistics_timeseries,
     _remove_less_covered,
+    _select_period_season_coldata,
 )
 from pyaerocom.aeroval.exceptions import TrendsError
 from pyaerocom.exceptions import TemporalResolutionError, UnknownRegion
@@ -339,3 +340,43 @@ def test__remove_less_covered_all():
     cd = COLDATA["fake_3d_partial_trends_coltime"]()
     new_cd = _remove_less_covered(cd, 1000)
     assert new_cd.data.station_name.shape[0] == 0
+
+
+@pytest.mark.parametrize(
+    "period, use_meteorological_seasons, resultfirst, resultlast",
+    [
+        (
+            "2004",
+            True,
+            np.datetime64("2003-12-01T00:00:00.000000000"),
+            np.datetime64("2004-02-29T00:00:00.000000000"),
+        ),
+        (
+            "2003-2004",
+            True,
+            np.datetime64("2003-01-01T00:00:00.000000000"),
+            np.datetime64("2004-12-31T00:00:00.000000000"),
+        ),
+        (
+            "2004",
+            False,
+            np.datetime64("2004-01-01T00:00:00.000000000"),
+            np.datetime64("2004-12-31T00:00:00.000000000"),
+        ),
+    ],
+)
+@pytest.mark.parametrize("coldataset", ["fake_3d_partial_trends"])
+def test__select_period_season_coldata(
+    coldata: ColocatedData,
+    period: str,
+    use_meteorological_seasons: bool,
+    resultfirst: np.datetime64,
+    resultlast: np.datetime64,
+):
+    coldatamod = coldata.copy()
+    coldatamod.data = coldatamod.data.assign_coords(
+        season=("time", coldata.coords["time.season"].values)
+    )
+    cd = _select_period_season_coldata(coldatamod, period, "DJF", use_meteorological_seasons)
+    assert cd.coords["time"].values[0] == resultfirst
+    assert cd.coords["time"].values[-1] == resultlast
diff --git a/tests/cams2_83/test_cams2_83_cli_mos.py b/tests/cams2_83/test_cams2_83_cli_mos.py
index 4017b87d..d0a8bb9b 100644
--- a/tests/cams2_83/test_cams2_83_cli_mos.py
+++ b/tests/cams2_83/test_cams2_83_cli_mos.py
@@ -23,9 +23,10 @@ def test_eval_mos_dummy(
     tmp_path: Path,
     caplog,
 ):
-    options = f"season 2024-03-01 2024-05-12 --data-path {tmp_path} --coldata-path {tmp_path} --name 'Test'"
+    options = f"season 2024-03-01 2024-05-12 --data-path {tmp_path} --coldata-path {tmp_path} --name 'Test' --addseasons"
     result = runner.invoke(app, options.split())
     assert result.exit_code == 0
+    assert "'add_seasons': True," in caplog.text
     assert "no output available" in caplog.text
 
 
diff --git a/tests/io/pyaro/test_read_pyaro.py b/tests/io/pyaro/test_read_pyaro.py
index 9ac8c865..eaa7bc7a 100644
--- a/tests/io/pyaro/test_read_pyaro.py
+++ b/tests/io/pyaro/test_read_pyaro.py
@@ -9,7 +9,7 @@ from pyaerocom.io import ReadPyaro, PyaroConfig
 from pyaerocom.io.pyaro.read_pyaro import PyaroToUngriddedData
 from pyaerocom.io.pyaro.postprocess import matching_indices
 
-from tests.conftest import lustre_unavail
+from tests.conftest import lustre_unavail, __package_installed
 
 
 def test_testfile(pyaro_test_data_file):
@@ -138,6 +138,10 @@ def test_matching_indices():
 
 
 @lustre_unavail
+@pytest.mark.skipif(
+    not __package_installed("pyaro_readers"),
+    reason="reader_id=eeareder requires pyaro-readers to be installed",
+)
 def test_vmrox():
     config = PyaroConfig.from_dict(
         {