Repository: metno/pyaerocom

Issue:
Revise FAIRMODE plots
The FAIRMODE plots should be based on hourly mean data for NO2 (and not on daily or monthly means as they are now), and for ozone they should be based on the 8h moving average daily maximum. 

Related issue: #863

Hints from issue comments:
@heikoklein: @charlienegri is trying to run the cams2_83 collocation with hourly statistics@lewisblake: Blocked by #1277@heikoklein: mda8 is implemented (@charlienegri )

no2 hourly currently takes too long time due to a performance bug in regional statistics #1277 blocking the implementation@heikoklein: @charlienegri made a new test of fairmode / hourly data:

- reading observations is no problem
- current production runs about 12h, all 6-7 variables (O3/O3MDA), 16.6G memory-consumption, only daily 
[cams283_analysis_long_1_JJA2024.log](https://github.com/user-attachments/files/17679528/cams283_analysis_long_1_JJA2024.log)
  - statistics comparinga against one model took 320s, equal share of about 1min for each of the individual statistical tasks
  - data-output (json files) is much larger, e.g. ts 15G, hm/ts 3.5G (only no2-files)
- test-run, only no2 but hourly runs about 12.5h, 17.1G max memory consumption
[cams283_test_forecast_long_no2only_hourly.log](https://github.com/user-attachments/files/17679530/cams283_test_forecast_long_no2only_hourly.log)
  - statistics comparing against a single model took ~5000s (most in Processing statistics timeseries for all regions 50min but also 20min in Processing individual site timeseries data)
  - data-output (json files) is much smaller, e.g. ts 0.8G, hm/ts 0.2G (only no2-files)

Speed is partly connected to the output-size, so just disabling the hourly timeseries (only daily) might speed things up a lot. @michaelgau We should sit together and discuss, which output needs to be available hourly. For (most) all others plots, we should stick to daily. Also plotting hourly data might get very slow of so huge datasets.
@michaelgau: ok, thanks! Will look at this when I'm back from the CAMS2_61 meeting (14 Nov)@heikoklein: After discussion with @michaelgau

For the "last seasons" plot, hourly output is in general not needed. Only the Fairmode NO2 plots should be based on hourly data.
In particular:
Fairmode plots are only defined for **one** time resolution:

- daily:
   - O3-MDA
   - PM10
   - PM2.5
- hourly:
   - NO2

e.g. monthly PM10 plots are wrong because they are based on threshold values which are only valid for daily data etc. This means, a time-resolution selector for Target-plots does not make sense, all undefined plots should be removed (as they are already for CO/SO2).  @charlienegri: yeah Michael has explained us this before, my understanding was that it was up to the users to select the frequency that makes sense, if available. if meaningless selections can be hidden front-end side, great.  

about NO2, so far there is no way of running hourly for fairmode statistics only @michaelgau: just for the record: FAIRMODE plots do exist also for annual time resolutions, but these are not relevant for us because we don't need to make annual evaluations on CAMS2_83 aeroval.@lewisblake: Addressed in #1551 

Solution patch:
diff --git a/pyaerocom/aeroval/coldatatojson_helpers.py b/pyaerocom/aeroval/coldatatojson_helpers.py
index 70951f58..3a739c7b 100644
--- a/pyaerocom/aeroval/coldatatojson_helpers.py
+++ b/pyaerocom/aeroval/coldatatojson_helpers.py
@@ -1001,7 +1001,7 @@ def _process_map_and_scat(
                         if use_fairmode and freq != "yearly" and not np.isnan(obs_vals).all():
                             stats["mb"] = np.nanmean(mod_vals - obs_vals)
 
-                            stats["fairmode"] = fairmode_stats(obs_var, stats)
+                            stats["fairmode"] = fairmode_stats(obs_var, stats, freq)
 
                         #  Code for the calculation of trends
                         if add_trends and freq != "daily":
diff --git a/pyaerocom/aeroval/fairmode_stats.py b/pyaerocom/aeroval/fairmode_stats.py
index bef60d99..24966b68 100644
--- a/pyaerocom/aeroval/fairmode_stats.py
+++ b/pyaerocom/aeroval/fairmode_stats.py
@@ -14,7 +14,6 @@ import numpy as np
 
 SPECIES = dict(
     concno2=dict(UrRV=0.24, RV=200, alpha=0.2),
-    conco3=dict(UrRV=0.18, RV=120, alpha=0.79),
     conco3mda8=dict(UrRV=0.18, RV=120, alpha=0.79),
     concpm10=dict(UrRV=0.28, RV=50, alpha=0.25),
     concpm25=dict(UrRV=0.36, RV=25, alpha=0.5),
@@ -53,10 +52,18 @@ def _mqi(rms: float, rmsu: float, *, beta: float) -> float:
     return rms / (rmsu * beta)
 
 
-def fairmode_stats(obs_var: str, stats: dict) -> dict:
+def fairmode_stats(obs_var: str, stats: dict, freq: str) -> dict:
     if obs_var not in SPECIES or np.isnan(list(stats.values())).any():
         return {}
 
+    # compute only what it makes sense to compute
+    if obs_var == "concno2":
+        if freq != "hourly":
+            return {}
+    else:
+        if freq != "daily":
+            return {}
+
     mean = stats["refdata_mean"]
     obs_std = stats["refdata_std"]
     mod_std = stats["data_std"]


Test patch:
diff --git a/tests/aeroval/test_fairmode_stats.py b/tests/aeroval/test_fairmode_stats.py
index e7de440d..89294f17 100644
--- a/tests/aeroval/test_fairmode_stats.py
+++ b/tests/aeroval/test_fairmode_stats.py
@@ -9,67 +9,93 @@ FAIRMODE_KEYS = {"RMSU", "sign", "crms", "bias", "rms", "alpha", "UrRV", "RV", "
 
 
 @pytest.mark.parametrize(
-    "obs_var,stats",
+    "obs_var,stats,freq",
     [
         pytest.param(
             "concno2",
             dict(refdata_mean=0, refdata_std=1, data_std=1, R=1, mb=0, rms=0),
-            id="dummy vars",
-        )
+            "hourly",
+            id="dummy_no2",
+        ),
+        pytest.param(
+            "conco3mda8",
+            dict(refdata_mean=0, refdata_std=1, data_std=1, R=1, mb=0, rms=0),
+            "daily",
+            id="dummy_o3mda8",
+        ),
     ],
 )
-def test_fairmode_stats(obs_var: str, stats: dict):
-    fairmode = fairmode_stats(obs_var, stats)
+def test_fairmode_stats(obs_var: str, stats: dict, freq: str):
+    fairmode = fairmode_stats(obs_var, stats, freq)
     assert set(fairmode) == FAIRMODE_KEYS
 
 
 @pytest.mark.parametrize(
-    "obs_var,stats",
+    "obs_var,stats,freq",
     [
         pytest.param(
-            "not_a_species",
+            "conco3",
+            dict(refdata_mean=0, refdata_std=1, data_std=1, R=1, mb=0, rms=0),
+            "daily",
+            id="not a valid species",
+        ),
+        pytest.param(
+            "concno2",
             dict(refdata_mean=0, refdata_std=1, data_std=1, R=1, mb=0, rms=0),
-            id="unknown obs variable",
+            "daily",
+            id="wrong frequency no2",
         ),
         pytest.param(
-            "vmro3",
+            "concpm10",
+            dict(refdata_mean=0, refdata_std=1, data_std=1, R=1, mb=0, rms=0),
+            "hourly",
+            id="wrong frequency pm10",
+        ),
+        pytest.param(
+            "concpm10",
             dict(refdata_mean=np.nan, refdata_std=1, data_std=1, R=1, mb=0, rms=0),
+            "daily",
             id="NaN mean",
         ),
         pytest.param(
-            "conco3",
+            "concpm10",
             dict(refdata_mean=0, refdata_std=np.nan, data_std=1, R=1, mb=0, rms=0),
+            "daily",
             id="NaN obs_std",
         ),
         pytest.param(
-            "conco3",
+            "concpm10",
             dict(refdata_mean=0, refdata_std=1, data_std=np.nan, R=1, mb=0, rms=0),
+            "daily",
             id="NaN mod_std",
         ),
         pytest.param(
-            "conco3",
+            "concpm10",
             dict(refdata_mean=0, refdata_std=1, data_std=1, R=np.nan, mb=0, rms=0),
+            "daily",
             id="NaN R",
         ),
         pytest.param(
-            "conco3",
+            "concpm10",
             dict(refdata_mean=0, refdata_std=1, data_std=1, R=1, mb=np.nan, rms=0),
+            "daily",
             id="NaN bias",
         ),
         pytest.param(
-            "vmro3",
+            "concpm10",
             dict(refdata_mean=0, refdata_std=1, data_std=1, R=1, mb=1, rms=np.nan),
+            "daily",
             id="NaN rms",
         ),
     ],
 )
-def test_empty_stats(obs_var: str, stats: dict):
-    fairmode = fairmode_stats(obs_var, stats)
+def test_empty_stats(obs_var: str, stats: dict, freq: str):
+    fairmode = fairmode_stats(obs_var, stats, freq)
     assert not fairmode
 
 
 @pytest.mark.parametrize(
-    "obs_var,stats,error",
+    "obs_var,stats,freq,error",
     [
         pytest.param(
             "concpm10",
@@ -81,6 +107,7 @@ def test_empty_stats(obs_var: str, stats: dict):
                 "mb": 1,
                 "rms": 1,
             },
+            "daily",
             "negative obs_std=-1",
             id="obs_std",
         ),
@@ -94,6 +121,7 @@ def test_empty_stats(obs_var: str, stats: dict):
                 "mb": 1,
                 "rms": 1,
             },
+            "daily",
             "negative mod_std=-1",
             id="mod_std",
         ),
@@ -107,12 +135,13 @@ def test_empty_stats(obs_var: str, stats: dict):
                 "mb": 1,
                 "rms": 1,
             },
+            "daily",
             "out of range R=10",
             id="R",
         ),
     ],
 )
-def test_fairmode_error(obs_var: str, stats: dict, error: str):
+def test_fairmode_error(obs_var: str, stats: dict, freq: str, error: str):
     with pytest.raises(AssertionError) as e:
-        fairmode_stats(obs_var, stats)
+        fairmode_stats(obs_var, stats, freq)
     assert str(e.value) == error