Repository: metno/pyaerocom

Issue:
fairmode output
the quantities used to far for the "line" charts at the bottom of the targets tab are not what they are meant to be   

all that is needed by the frontend to create these lines can and will be added to the fairmode files

Discussion:

Solution patch:
diff --git a/pyaerocom/aeroval/fairmode_statistics.py b/pyaerocom/aeroval/fairmode_statistics.py
index e24ac6ef..4f6421dd 100644
--- a/pyaerocom/aeroval/fairmode_statistics.py
+++ b/pyaerocom/aeroval/fairmode_statistics.py
@@ -12,10 +12,25 @@ logger = logging.getLogger(__name__)
 
 
 SPECIES = dict(
-    concno2=dict(UrRV=0.24, RV=200, alpha=0.2, freq=TsType("hourly"), percentile=99.8),
-    conco3mda8=dict(UrRV=0.18, RV=120, alpha=0.79, freq=TsType("daily"), percentile=92.9),
-    concpm10=dict(UrRV=0.28, RV=50, alpha=0.25, freq=TsType("daily"), percentile=90.1),
-    concpm25=dict(UrRV=0.36, RV=25, alpha=0.5, freq=TsType("daily"), percentile=90.1),
+    concno2=dict(
+        UrRV=0.24, RV=200, alpha=0.2, freq=TsType("hourly"), percentile=99.8, Np=5.2, Nnp=5.5
+    ),
+    conco3mda8=dict(
+        UrRV=0.18, RV=120, alpha=0.79, freq=TsType("daily"), percentile=92.9, Np=11.0, Nnp=3.0
+    ),
+    concpm10=dict(
+        UrRV=0.28, RV=50, alpha=0.25, freq=TsType("daily"), percentile=90.1, Np=20.0, Nnp=1.5
+    ),
+    concpm25=dict(
+        UrRV=0.36, RV=25, alpha=0.5, freq=TsType("daily"), percentile=90.1, Np=20.0, Nnp=1.5
+    ),
+)
+
+EXC_THRESHOLDS = dict(  # we assume all the units are ug/m3
+    concpm25=25.0,
+    concpm10=50.0,
+    conco3mda8=120.0,
+    concno2=200.0,
 )
 
 
@@ -70,9 +85,10 @@ class FairmodeStatistics:
         obsvals = data.data[0]
         modvals = data.data[1]
 
-        mask = ~np.isnan(obsvals) * ~np.isnan(modvals)
+        mask = np.isfinite(obsvals) & np.isfinite(modvals)
 
-        obsmean = np.nanmean(obsvals, axis=0)
+        obsmean = np.mean(obsvals, axis=0, where=mask)
+        modmean = np.mean(modvals, axis=0, where=mask)
 
         obsstd = np.std(obsvals, axis=0, where=mask)
         modstd = np.std(modvals, axis=0, where=mask)
@@ -83,15 +99,26 @@ class FairmodeStatistics:
         rms = np.sqrt(np.nanmean(diffsquare, axis=0, where=mask))
         bias = np.nanmean(diff, axis=0, where=mask)
 
+        βRMSUt = self._βRMSU_t(obsvals, beta=1, var_name=var_name, mask=mask)
+        βRMSUs = self._βRMSU_s(obsmean, beta=1, var_name=var_name)
+
+        NMB = self._NMB(modvals, obsvals)
         R = self.pearson_R(obsvals, modvals)
-        rmsu = self._RMSU(obsmean, obsstd, var_name)
         sign = self._fairmode_sign(modstd, obsstd, R)
         crms = self._crms(modstd, obsstd, R)
-        mqi = self._mqi(rms, rmsu, beta=1)
-        mb = self._mb(bias, rmsu, beta=1)
+        mqi = self._mqi(rms, βRMSUt, beta=1)
+        mb = self._mb(bias, βRMSUt, beta=1)
         beta_Hperc = self._beta_Hperc(obsvals, modvals, var_name)
+        exceedances = self._exceedances(data=data, var_name=var_name)
+
+        MPI_bias_t = self._MPI_bias_t(obsmean, modmean, βRMSUt)
+        MPI_R_t = self._MPI_R_t(obsstd, modstd, R, βRMSUt)
+        MPI_std_t = self._MPI_std_t(obsstd, modstd, βRMSUt)
+
+        MPI_R_s = self._MPI_R_s(obsmean, modmean, βRMSUs)
+        MPI_std_s = self._MPI_std_s(obsmean, modmean, βRMSUs)
 
-        assert len(rmsu) == len(stations)
+        assert len(βRMSUt) == len(stations)
         assert len(sign) == len(stations)
         assert len(crms) == len(stations)
         assert len(bias) == len(stations)
@@ -102,13 +129,21 @@ class FairmodeStatistics:
 
         stats_list: dict[str, dict[str, float]] = {
             stations[i]: dict(
-                RMSU=rmsu[i],
+                exceedances_obs=int(exceedances[0][i]),
+                MPI_mean=obsmean[i],
+                MPI_R_t=MPI_R_t[i],
+                MPI_bias_t=MPI_bias_t[i],
+                MPI_std_t=MPI_std_t[i],
+                MPI_R_s=MPI_R_s,
+                MPI_std_s=MPI_std_s,
+                MPI_Hperc=beta_Hperc[i],
+                bias=bias[i],
+                NMB=NMB[i],
+                RMSU=βRMSUt[i],
                 sign=[sign[i]],
                 crms=crms[i],
-                bias=bias[i],
                 rms=[rms[i]],
                 beta_mqi=[mqi[i]],
-                Hperc=beta_Hperc[i],
                 persistence_model=False,
                 station_type=station_types[i],
                 **{k: (str(v) if k == "freq" else v) for (k, v) in SPECIES[var_name].items()},
@@ -118,9 +153,31 @@ class FairmodeStatistics:
 
         return stats_list
 
+    @staticmethod
+    def _exceedances(data: xr.DataArray, var_name: str) -> list[np.array]:
+        if var_name == "concno2":
+            new_data = data.resample(time="1D", skipna=True).max()
+        else:
+            new_data = data
+
+        obsvals = new_data.data[0]
+        modvals = new_data.data[1]
+
+        mask = np.isfinite(obsvals) & np.isfinite(modvals)
+        obsex = np.sum(obsvals > EXC_THRESHOLDS[var_name], axis=0, where=mask)
+        modex = np.sum(modvals > EXC_THRESHOLDS[var_name], axis=0, where=mask)
+
+        return [obsex, modex]
+
+    @staticmethod
+    def _NMB(x: np.ndarray, y: np.ndarray) -> np.ndarray:
+        num = np.nansum(x - y, axis=0)
+        denum = np.nansum(x, axis=0)
+        return np.where(denum == 0, np.nan, num / denum)
+
     @staticmethod
     def pearson_R(x: np.ndarray, y: np.ndarray) -> np.ndarray:
-        mask = ~np.isnan(x) * ~np.isnan(y)
+        mask = np.isfinite(x) & np.isfinite(y)
 
         xmean = np.mean(x, axis=0, where=mask)
         ymean = np.mean(y, axis=0, where=mask)
@@ -137,20 +194,105 @@ class FairmodeStatistics:
 
         return r
 
+    # in the following _s and _t are for the spatial(ly averaged) and temporal(ly-only averaged) quantities,
+    # we define a _βRMSU_t to be used for the calculation of (following MeteoFrance/evaltools naming)
+    # Time Bias Norm (_MPI_bias_t), Time Corr Norm (_MPI_r_t) and Time StdDev Norm (_MPI_std_t)
+    # and a _βRMSU_s to be used for the calculation of Space Corr Norm (_MPI_R_s) and Space StDev Norm (_MPI_std_s)
+
+    @staticmethod
+    def _βRMSU_t(
+        obsvals: np.ndarray,
+        beta: float,
+        var_name: str,
+        mask: np.ndarray,
+    ) -> np.ndarray:
+        def obsuncertainty(obs: np.ndarray, spec: str) -> np.ndarray:
+            """eq. 7 and 37 here https://fairmode.jrc.ec.europa.eu/document/fairmode/WG1/Guidance_MQO_Bench_vs3.3_20220519.pdf"""
+
+            if spec not in SPECIES:
+                raise ValueError(f"Unsupported {spec=}")
+
+            UrRV = SPECIES[spec]["UrRV"]
+            RV = SPECIES[spec]["RV"]
+            alpha = SPECIES[spec]["alpha"]
+
+            in_sqrt = (1 - alpha**2) * (obs**2) + alpha**2 * RV**2
+
+            return UrRV * np.sqrt(in_sqrt)
+
+        return beta * np.sqrt(
+            np.nanmean(np.square(obsuncertainty(obsvals, var_name), where=mask), axis=0)
+        )
+
+    @staticmethod
+    def _βRMSU_s(obsmean: np.ndarray, beta: float, var_name: str) -> float:
+        def obsuncertainty(obs: np.ndarray, spec: str) -> np.ndarray:
+            """eq. 39 here https://fairmode.jrc.ec.europa.eu/document/fairmode/WG1/Guidance_MQO_Bench_vs3.3_20220519.pdf"""
+
+            if spec not in SPECIES:
+                raise ValueError(f"Unsupported {spec=}")
+
+            UrRV = SPECIES[spec]["UrRV"]
+            RV = SPECIES[spec]["RV"]
+            alpha = SPECIES[spec]["alpha"]
+            N_p = SPECIES[spec]["Np"]
+            N_np = SPECIES[spec]["Nnp"]
+
+            in_sqrt = (1 - alpha**2) / N_p * (obs**2) + alpha**2 * RV**2 / N_np
+
+            return UrRV * np.sqrt(in_sqrt)
+
+        return beta * np.sqrt(np.nanmean(np.square(obsuncertainty(obsmean, var_name))))
+
     @staticmethod
-    def _RMSU(mean: float, std: float, spec: str) -> float:
-        """RMSU is the Root Mean Squared Uncertainty associated with the uncertainty of the observations, U(O_i)."""
+    def _MPI_R_t(obsstd: np.array, modstd: np.array, R: float, βRMSUt: np.array) -> np.array:
+        # TIME Corr Norm: 2 sigma_O sigma_M(1-R) / (beta^2 RMS_U^2)
+        # ----------------------------------------------------------
+        # MF eq.: (2.*scores['obs_std']*scores['sim_std']*(1. - scores['PearsonR'])) / (beta*rmsu)**2
+        return np.where(βRMSUt == 0, np.nan, 2 * modstd * obsstd * (1 - R) / βRMSUt**2)
 
-        if spec not in SPECIES:
-            raise ValueError(f"Unsupported {spec=}")
+    @staticmethod
+    def _MPI_bias_t(obsmean: np.array, modmean: np.array, βRMSUt: np.array) -> np.array:
+        # TIME Bias Norm: |BIAS| / (beta RMS_U)
+        # -------------------------------------------
+        # MF eq.: scores['MeanBias']/(beta*rmsu)
+        return np.where(
+            βRMSUt == 0, np.nan, (modmean - obsmean) / βRMSUt
+        )  # check the abs here, why NMF does not have it?
 
-        UrRV = SPECIES[spec]["UrRV"]
-        RV = SPECIES[spec]["RV"]
-        alpha = SPECIES[spec]["alpha"]
+    @staticmethod
+    def _MPI_std_t(obsstd: np.array, modstd: np.array, βRMSUt: np.array) -> np.array:
+        # TIME StDev Norm: (sigma_M-sigma_O) / (beta RMS_U)
+        # ---------------------------------------------------
+        # MF eq.: (scores['sim_std']-scores['obs_std'])/(beta*rmsu)
+        return np.where(βRMSUt == 0, np.nan, (modstd - obsstd) / βRMSUt)
 
-        in_sqrt = (1 - alpha**2) * (mean**2 + std**2) + alpha**2 * RV**2
+    @staticmethod
+    def _MPI_R_s(obsmean: np.array, modmean: np.array, βRMSUs: float) -> float:
+        # SPACE Corr Norm: 2 sigma_bar{O} sigma_bar{M} (1-R) / (beta^2 RMS_bar{U}^2)
+        # ----------------------------------------------------------------------------
+        # MF eq.: ((2.*np.nanstd(obs)*np.nanstd(sim)*(1. - corr)) / (beta*rmsu_)**2)
+        # where sim = scores['sim_mean'], obs = scores['obs_mean']
+        # corr = ((np.nanmean((obs-np.nanmean(obs))*(sim-np.nanmean(sim)))) / (np.nanstd(obs)*np.nanstd(sim)))
+        corr = (
+            0.0
+            if np.nanstd(obsmean) * np.nanstd(modmean) == 0
+            else np.nanmean((obsmean - np.nanmean(obsmean)) * (modmean - np.nanmean(modmean)))
+            / (np.nanstd(obsmean) * np.nanstd(modmean))
+        )
+        return (
+            2 * np.nanstd(obsmean) * np.nanstd(modmean) * (1.0 - corr) / βRMSUs**2
+            if βRMSUs != 0
+            else np.nan
+        )
 
-        return UrRV * np.sqrt(in_sqrt)
+    @staticmethod
+    def _MPI_std_s(obsmean: np.array, modmean: np.array, βRMSUs: float) -> float:
+        # SPACE StDev Norm: (sigma_bar{M}-sigma_bar{O}) / ( beta RMS_bar{U})
+        # -------------------------------------------------------------------
+        # MF eq.: (np.nanstd(sim)-np.nanstd(obs))/(beta*rmsu_)
+        # where sim = scores['sim_mean'], obs = scores['obs_mean']
+        return (np.nanstd(modmean) - np.nanstd(obsmean)) / βRMSUs if βRMSUs != 0 else np.nan
 
     @staticmethod
     def _fairmode_sign(mod_std: float, obs_std: float, R: float) -> float:
@@ -187,4 +329,4 @@ class FairmodeStatistics:
             (1 - SPECIES[var_name]["alpha"] ** 2) * Operc**2 + factor
         )
 
-        return (Operc - Mperc) / (beta * uncertainty_Operc)
+        return (Mperc - Operc) / (beta * uncertainty_Operc)
diff --git a/pyaerocom/colocation/colocator.py b/pyaerocom/colocation/colocator.py
index 03799af1..426cd755 100644
--- a/pyaerocom/colocation/colocator.py
+++ b/pyaerocom/colocation/colocator.py
@@ -508,7 +508,6 @@ class Colocator:
             filter_post=obs_filters_post,
             **self.colocation_setup.read_opts_ungridded,
         )
-        obs_data.check_unit(var_name)
 
         if self.colocation_setup.obs_remove_outliers:
             oor = self.colocation_setup.obs_outlier_ranges
diff --git a/pyaerocom/io/cachehandler_ungridded.py b/pyaerocom/io/cachehandler_ungridded.py
index 87042b3b..b610eefe 100644
--- a/pyaerocom/io/cachehandler_ungridded.py
+++ b/pyaerocom/io/cachehandler_ungridded.py
@@ -14,6 +14,7 @@ from pyaerocom.exceptions import CacheReadError, CacheWriteError
 from pyaerocom.ungridded_data_container import UngriddedDataContainer
 from pyaerocom.ungriddeddata import UngriddedData
 from pyaerocom.ungriddeddata_structured import UngriddedDataStructured
+from pyaerocom.io.pyaro.read_pyaro import ReadPyaro
 
 logger = logging.getLogger(__name__)
 
@@ -57,7 +58,7 @@ class CacheHandlerUngridded:
     def __init__(self, reader=None, cache_dir=None, **kwargs):
         self._reader = None
         if reader is not None:
-            self.reader = reader
+            self._reader = reader
 
         self.loaded_data = {}
         self._cache_dir = cache_dir
@@ -207,6 +208,8 @@ class CacheHandlerUngridded:
         current["ungridded_data_class"] = data_classname
         current["ungridded_data_version"] = dataclass.__version__
         current["cacher_version"] = self.__version__
+        if isinstance(self._reader, ReadPyaro):
+            current["pyaro_config"] = self.reader.config
         return current
 
     def check_and_load(self, var_or_file_name, force_use_outdated=False, cache_dir=None):
diff --git a/pyaerocom/io/mscw_ctm/reader.py b/pyaerocom/io/mscw_ctm/reader.py
index 894122b7..ec272adc 100755
--- a/pyaerocom/io/mscw_ctm/reader.py
+++ b/pyaerocom/io/mscw_ctm/reader.py
@@ -13,6 +13,7 @@ from pyaerocom.griddeddata import GriddedData
 from pyaerocom.io.gridded_reader import GriddedReader
 from pyaerocom.projection_information import ProjectionInformation
 from pyaerocom.units.helpers import get_standard_unit
+from pyaerocom.units.units import Unit
 
 from .additional_variables import (
     add_dataarrays,
@@ -42,7 +43,6 @@ from .additional_variables import (
 from .model_variables import emep_variables
 import pathlib
 
-
 logger = logging.getLogger(__name__)
 
 
@@ -781,7 +781,8 @@ class ReadMscwCtm(GriddedReader):
         ts_type = self._ts_type
 
         arr, proj_info = self._load_var(var_name_aerocom, ts_type)
-        arr.attrs["units"] = arr.units
+        if arr.units in Unit._UALIASES:
+            arr.attrs["units"] = Unit._UALIASES[arr.units]
         try:
             cube = arr.to_iris()
         except MemoryError as e:  # pragma: no cover
diff --git a/pyaerocom/io/pyaro/postprocess.py b/pyaerocom/io/pyaro/postprocess.py
index 8d6aa65e..9d3c76bf 100644
--- a/pyaerocom/io/pyaro/postprocess.py
+++ b/pyaerocom/io/pyaro/postprocess.py
@@ -120,6 +120,14 @@ class PostProcessingReaderData(Data):
     def __len__(self):
         return self.data.__len__()
 
+    @property
+    def variable(self) -> str:
+        return self._variable
+
+    @property
+    def units(self) -> str:
+        return self._units
+
     @property
     def values(self):
         if self.scaling is None:
@@ -204,7 +212,9 @@ class PostProcessingReader(Reader):
         if isinstance(transform, VariableScaling):
             data = self.reader.data(transform.REQ_VAR)
             scaling = transform.SCALING_FACTOR * get_unit_conversion_fac(
-                from_unit=data.units, to_unit=transform.IN_UNIT, var_name=transform.REQ_VAR
+                from_unit=data.units,
+                to_unit=transform.IN_UNIT,
+                var_name=transform.REQ_VAR,
             )
             return PostProcessingReaderData(
                 data, variable=varname, units=transform.OUT_UNIT, scaling=scaling
@@ -350,13 +360,21 @@ class DictBackedData(Data):
     def slice(self, index):
         return DictBackedData(
             data=self._data()[index],
-            variable=self._variable,
+            variable=self.variable,
             units=self._units,
         )
 
     def __len__(self):
         return len(self.values)
 
+    @property
+    def variable(self) -> str:
+        return self._variable
+
+    @property
+    def units(self) -> str:
+        return self._units
+
     @property
     def values(self):
         return self._data["values"]
diff --git a/pyaerocom/io/readungridded.py b/pyaerocom/io/readungridded.py
index 05912046..a49744e6 100755
--- a/pyaerocom/io/readungridded.py
+++ b/pyaerocom/io/readungridded.py
@@ -513,7 +513,7 @@ class ReadUngridded:
                 f"supported by {data_id} interface"
             )
         cache = CacheHandlerUngridded(reader)
-        if not self.ignore_cache and not isinstance(cache.reader, ReadPyaro):
+        if not self.ignore_cache:
             # initiate cache handler
             for var in vars_available:
                 try:
@@ -538,7 +538,7 @@ class ReadUngridded:
 
             for var in vars_to_read:
                 # write the cache file
-                if not self.ignore_cache and not isinstance(cache.reader, ReadPyaro):
+                if not self.ignore_cache:
                     try:
                         cache.write(data_read, var)
                     except Exception as e:
diff --git a/pyaerocom/sample_data_access/minimal_dataset.py b/pyaerocom/sample_data_access/minimal_dataset.py
index 0e9e5808..bf305bc0 100644
--- a/pyaerocom/sample_data_access/minimal_dataset.py
+++ b/pyaerocom/sample_data_access/minimal_dataset.py
@@ -11,7 +11,7 @@ logger = logging.getLogger(__name__)
 __all__ = ["download_minimal_dataset"]
 
 #: tarfile to download
-DEFAULT_TESTDATA_FILE = "testdata-minimal.tar.gz.20250506"
+DEFAULT_TESTDATA_FILE = "testdata-minimal.tar.gz.20250512"
 
 minimal_dataset = pooch.create(
     path=const.OUTPUTDIR,  # ~/MyPyaerocom/
@@ -28,6 +28,7 @@ minimal_dataset = pooch.create(
         "testdata-minimal.tar.gz.20241120": "md5:4d2bc1782b1f468321817139d327e014",
         "testdata-minimal.tar.gz.20250425": "md5:23f5b2e34f294c3232248a2b9b779864",
         "testdata-minimal.tar.gz.20250506": "md5:aab174c263d350e9c6120614a0bda8a5",
+        "testdata-minimal.tar.gz.20250512": "md5:00d2f7cf41e6303bad33e822d25fe960",
     },
 )
 
diff --git a/pyaerocom/scripts/cams2_83/engine.py b/pyaerocom/scripts/cams2_83/engine.py
index 44e0200b..516e2846 100644
--- a/pyaerocom/scripts/cams2_83/engine.py
+++ b/pyaerocom/scripts/cams2_83/engine.py
@@ -80,6 +80,7 @@ class CAMS2_83_Engine(ProcessingEngine):
         use_weights = self.cfg.statistics_opts.weighted_stats
         forecast_days = self.cfg.statistics_opts.forecast_days
         periods = self.cfg.time_cfg.periods
+        min_num_obs = self.cfg.colocation_opts.min_num_obs
 
         # use_fairmode = self.cfg.statistics_opts.use_fairmode
         use_fairmode = self.cfg.cams2_83_cfg.use_cams2_83_fairmode
@@ -98,6 +99,7 @@ class CAMS2_83_Engine(ProcessingEngine):
                 persistence_coldata = persistence_coldata.resample_time(
                     SPECIES[var_name]["freq"],
                     settings_from_meta=True,
+                    min_num_obs=min_num_obs,
                 )
 
         if "var_name_input" in coldata[0].metadata:
@@ -202,6 +204,7 @@ class CAMS2_83_Engine(ProcessingEngine):
                             fairmode_subset = fairmode_subset.resample_time(
                                 SPECIES[var_name]["freq"],
                                 settings_from_meta=True,
+                                min_num_obs=min_num_obs
                             )
 
                         results_fairmode[f"{regname}"][f"{perstr}"] = (
@@ -219,7 +222,7 @@ class CAMS2_83_Engine(ProcessingEngine):
 
                                 
                                 mqi_results = self._calc_forecast_target_MQI_vectorized(
-                                    ds, ds_p, var_name, day
+                                    ds, ds_p, var_name, day, min_num_obs
                                 )
 
                                 results_mqi.append(mqi_results)
@@ -352,6 +355,7 @@ class CAMS2_83_Engine(ProcessingEngine):
         persistence_coldata: ColocatedData,
         var_name: str,
         forecast_day: int,
+        min_num_obs: dict,
     ) -> dict[str, float]:
 
         results = {}
@@ -361,6 +365,7 @@ class CAMS2_83_Engine(ProcessingEngine):
             coldata = coldata.resample_time(
                 SPECIES[var_name]["freq"],
                 settings_from_meta=True,
+                min_num_obs=min_num_obs
             )
 
         # Creation of mask of shared stations between normal data and persistence data
diff --git a/pyaerocom/stationdata.py b/pyaerocom/stationdata.py
index 4197c5b5..f56e6f68 100644
--- a/pyaerocom/stationdata.py
+++ b/pyaerocom/stationdata.py
@@ -980,6 +980,7 @@ class StationData(StationMetaData):
         var_name: str,
         low: float | None = None,
         high: float | None = None,
+        unit_ref=None,
         check_unit: bool = True,
     ):
         """Remove outliers from one of the variable timeseries
diff --git a/pyaerocom/ungridded_data_container.py b/pyaerocom/ungridded_data_container.py
index bb78d804..d7f2af44 100644
--- a/pyaerocom/ungridded_data_container.py
+++ b/pyaerocom/ungridded_data_container.py
@@ -477,6 +477,7 @@ class UngriddedDataContainer(abc.ABC):
         inplace=False,
         low=None,
         high=None,
+        unit_ref=None,
         move_to_trash=True,
     ):
         """Method that can be used to remove outliers from data
diff --git a/pyaerocom/ungriddeddata.py b/pyaerocom/ungriddeddata.py
index 4295ba45..40a79e2b 100644
--- a/pyaerocom/ungriddeddata.py
+++ b/pyaerocom/ungriddeddata.py
@@ -1023,6 +1023,7 @@ class UngriddedData(UngriddedDataMetadata):
         inplace=False,
         low=None,
         high=None,
+        unit_ref=None,
         move_to_trash=True,
     ):
         """Method that can be used to remove outliers from data
@@ -1069,7 +1070,7 @@ class UngriddedData(UngriddedDataMetadata):
         else:
             new = self.copy()
 
-        new.check_unit(var_name)
+        new.check_convert_var_units(var_name, to_unit=unit_ref)
 
         if low is None:
             low = const.VARS[var_name].minimum
diff --git a/pyaerocom/ungriddeddata_structured.py b/pyaerocom/ungriddeddata_structured.py
index 1cb098da..38544cf7 100644
--- a/pyaerocom/ungriddeddata_structured.py
+++ b/pyaerocom/ungriddeddata_structured.py
@@ -14,6 +14,7 @@ from pyaerocom.dynamic_rec_array import DynamicRecArray
 from pyaerocom.exceptions import (
     DataCoverageError,
     DataExtractionError,
+    MetaDataError,
     VarNotAvailableError,
 )
 from pyaerocom.helpers import merge_station_data, start_stop
@@ -22,7 +23,7 @@ from pyaerocom.stationdata import StationData
 from pyaerocom.ungridded_data_container import UngriddedDataContainer
 from pyaerocom.ungridded_data_metadata import UngriddedDataMetadata
 from pyaerocom.units.datetime import TsType
-from pyaerocom.units.helpers import get_standard_unit
+from pyaerocom.units.units_helpers import get_unit_conversion_fac
 
 if sys.version_info >= (3, 12):
     from typing import override
@@ -457,6 +458,39 @@ class UngriddedDataStructured(UngriddedDataMetadata):
                 _iter.append(stat_name)
         return _iter
 
+    def check_convert_var_units(self, var_name, to_unit=None, inplace=True):
+        obj = self if inplace else self.copy()
+
+        # get the unit
+        if to_unit is None:
+            to_unit = const.VARS[var_name]["units"]
+
+        for meta_idx, meta in obj.metadata.items():
+            if var_name in meta["var_info"]:
+                try:
+                    unit = meta["var_info"][var_name]["units"]
+                except KeyError:
+                    add_str = ""
+                    if "unit" in meta["var_info"][var_name]:
+                        add_str = (
+                            "Corresponding var_info dict contains "
+                            'attr. "unit", which is deprecated, please '
+                            "check corresponding reading routine. "
+                        )
+                    raise MetaDataError(
+                        f"Failed to access unit information for variable {var_name} "
+                        f"in metadata block {meta_idx}. {add_str}"
+                    )
+                fac = get_unit_conversion_fac(unit, to_unit, var_name)
+                if fac != 1:
+                    idx = (obj._dra.data["meta_id"] == meta_idx) & (
+                        obj._dra.data["var_id"] == self.var_idx[var_name]
+                    )
+                    obj._dra.data["data"][idx] *= fac
+                meta["var_info"][var_name]["units"] = to_unit
+
+        return obj
+
     @override
     def remove_outliers(
         self,
@@ -464,6 +498,7 @@ class UngriddedDataStructured(UngriddedDataMetadata):
         inplace=False,
         low=None,
         high=None,
+        unit_ref=None,
         move_to_trash=True,
     ):
         """see super.remove_outliers move_to_trash is ignored for now"""
@@ -472,7 +507,7 @@ class UngriddedDataStructured(UngriddedDataMetadata):
         else:
             new = self.copy()
 
-        new.check_unit(var_name)
+        new.check_convert_var_units(var_name, to_unit=unit_ref)
 
         if low is None:
             low = const.VARS[var_name].minimum
@@ -605,12 +640,6 @@ class UngriddedDataStructured(UngriddedDataMetadata):
                 self.metadata[meta_idx]["var_info"] = {}
                 self.metadata[meta_idx]["var_info"][var] = {}
 
-                if "units" in self.metadata[meta_idx]["var_info"][var]:
-                    u1 = self.metadata[meta_idx]["var_info"][var]
-                    # u2 = station_data["var_info"][var]["units"]
-                    station_data.convert_unit(var, u1)
-                else:
-                    station_data.convert_unit(var, get_standard_unit(var))
                 self.metadata[meta_idx]["var_info"][var].update(station_data["var_info"][var])
                 for x in ("longitude", "latitude", "altitude"):
                     if x not in self.metadata[meta_idx]["var_info"][var]:
diff --git a/pyaerocom_env.yml b/pyaerocom_env.yml
index ff9a36f4..dcdff243 100644
--- a/pyaerocom_env.yml
+++ b/pyaerocom_env.yml
@@ -33,7 +33,7 @@ dependencies:
   - pip:
     - geocoder_reverse_natural_earth >= 0.0.2
     - aerovaldb >= 0.4.6, < 0.5.0
-    - pyaro >= 0.1.3
+    - pyaro >= 0.2.0, < 0.3
 ## testing
   - pytest >=7.4
   - pytest-dependency
diff --git a/pyproject.toml b/pyproject.toml
index 5cdb34fe..05d5a2ad 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -49,7 +49,7 @@ dependencies = [
     'cf-units>=3.3.0',
     "pydantic~=2.7",
     "pyproj>=3.0.0",
-    "pyaro>=0.1.3",
+    "pyaro>=0.2.0, <0.3.0",
     "pooch>=1.7.0",
     "psutil>=5.0.0",
     "contourpy>=1.3.0",
@@ -242,7 +242,7 @@ deps =
     typing-extensions ==4.6.1; python_version < "3.11"
     cf-units ==3.3.0; python_version < "3.11"
     pydantic ==2.7.0; python_version < "3.11"
-    pyaro == 0.1.3; python_version < "3.11"
+    pyaro == 0.2.0; python_version < "3.11"
     pooch ==1.7.0; python_version < "3.11"
     xarray ==2022.12.0; python_version < "3.11"
     pandas ==1.5.3; python_version < "3.11"


Test patch:
diff --git a/tests/aeroval/test_coldatatojson_helpers2.py b/tests/aeroval/test_coldatatojson_helpers2.py
index 6ef1a66c..959d12b8 100644
--- a/tests/aeroval/test_coldatatojson_helpers2.py
+++ b/tests/aeroval/test_coldatatojson_helpers2.py
@@ -433,7 +433,7 @@ def test_calculate_fairmode(eval_config: dict, caplog):
     )
     assert all(
         item in results["ALL"][f"{period}-{season}"]["Agoufou"]
-        for item in ["RMSU", "sign", "beta_mqi", "Hperc", "crms", "bias", "rms"]
+        for item in ["RMSU", "sign", "beta_mqi", "MPI_Hperc", "crms", "rms"]
     )
 
     wrongperiod = "2025"
diff --git a/tests/aeroval/test_fairmode_statistics.py b/tests/aeroval/test_fairmode_statistics.py
index f018c93f..24009321 100644
--- a/tests/aeroval/test_fairmode_statistics.py
+++ b/tests/aeroval/test_fairmode_statistics.py
@@ -1,5 +1,6 @@
 import numpy as np
 import pytest
+import xarray as xr
 
 # from pyaerocom import ColocatedData, Colocator
 from pyaerocom.aeroval import EvalSetup  # , ExperimentProcessor
@@ -59,12 +60,28 @@ def test_fairmode_statistics(fairmode_statistics, dummy_coldata_to_fairmode_stat
     )
     assert all(
         item in fm_stats["Agoufou"]
-        for item in ["RMSU", "sign", "beta_mqi", "Hperc", "crms", "bias", "rms"]
+        for item in [
+            "exceedances_obs",
+            "MPI_mean",
+            "MPI_R_t",
+            "MPI_bias_t",
+            "MPI_std_t",
+            "MPI_R_s",
+            "MPI_std_s",
+            "MPI_Hperc",
+            "bias",
+            "NMB",
+            "RMSU",
+            "sign",
+            "beta_mqi",
+            "crms",
+            "rms",
+        ]
     )
 
 
 def test_fairmode_statistics_wrongspecies(
-    fairmode_statistics, dummy_coldata_to_fairmode_statistics, caplog
+    fairmode_statistics, dummy_coldata_to_fairmode_statistics
 ):
     wrongspec = "concco"
     example_coldata = dummy_coldata_to_fairmode_statistics
@@ -75,19 +92,46 @@ def test_fairmode_statistics_wrongspecies(
     assert f"Unsupported spec='{wrongspec}'" in str(e.value)
 
 
+def test_βRMSUt_wrongspecies(fairmode_statistics, dummy_coldata_to_fairmode_statistics):
+    obsvals = dummy_coldata_to_fairmode_statistics.data[0]
+    wrongspec = "concso2"
+
+    with pytest.raises(ValueError) as e:
+        fairmode_statistics._βRMSU_t(obsvals, beta=1, var_name=wrongspec, mask=None)
+    assert f"Unsupported spec='{wrongspec}'" in str(e.value)
+
+
+def test_βRMSUs_wrongspecies(fairmode_statistics, dummy_coldata_to_fairmode_statistics):
+    obsvals = dummy_coldata_to_fairmode_statistics.data[0]
+    obsmean = np.nanmean(obsvals, axis=0)
+    wrongspec = "concso2"
+
+    with pytest.raises(ValueError) as e:
+        fairmode_statistics._βRMSU_s(obsmean, beta=1, var_name=wrongspec)
+    assert f"Unsupported spec='{wrongspec}'" in str(e.value)
+
+
 @pytest.fixture
 def fairmode_stats_example() -> dict:
     return {
         "ALL": {
             "2010-DJF": {
                 "Alta_Floresta": {
-                    "RMSU": np.float64(9.600056476313366),
+                    "exceedances_obs": 0,
+                    "MPI_mean": np.float64(0.040463086151092224),
+                    "MPI_R_t": np.float64(0.0005913779651058249),
+                    "MPI_bias_t": np.float64(9.168897156799856e-05),
+                    "MPI_std_t": np.float64(0.01621269226072656),
+                    "MPI_R_s": np.float64(0.001440664233921059),
+                    "MPI_std_s": np.float64(-0.031303579243870915),
+                    "MPI_Hperc": np.float64(0.054898295733501835),
+                    "bias": np.float64(0.0008803075523111548),
+                    "NMB": np.float64(0.0020848533921114204),
+                    "RMSU": np.float64(9.601018936702271),
                     "sign": [np.float64(1.0)],
-                    "crms": np.float64(0.02482145009443063),
-                    "bias": np.float64(-0.0593675912149444),
-                    "rms": [np.float64(0.064347612787539)],
-                    "beta_mqi": [np.float64(0.006702836899585607)],
-                    "Hperc": np.float64(0.00965331794743008),
+                    "crms": np.float64(0.2806109229399387),
+                    "rms": [np.float64(0.2806123037476991)],
+                    "beta_mqi": [np.float64(0.029227346138750866)],
                     "persistence_model": False,
                     "station_type": np.str_("bla"),
                     "UrRV": 0.24,
@@ -95,15 +139,25 @@ def fairmode_stats_example() -> dict:
                     "alpha": 0.2,
                     "freq": "hourly",
                     "percentile": 99.8,
+                    "Np": 5.2,
+                    "Nnp": 5.5,
                 },
                 "Thessaloniki": {
-                    "RMSU": np.float64(9.600986774272528),
+                    "exceedances_obs": 0,
+                    "MPI_mean": np.float64(0.019649444323128436),
+                    "MPI_R_t": np.float64(1.910032059129425e-05),
+                    "MPI_bias_t": np.float64(0.005486454983897092),
+                    "MPI_std_t": np.float64(0.001172539992895984),
+                    "MPI_R_s": np.float64(0.001440664233921059),
+                    "MPI_std_s": np.float64(-0.031303579243870915),
+                    "MPI_Hperc": np.float64(0.0019844616931031385),
+                    "bias": np.float64(0.05267057812092559),
+                    "NMB": np.float64(0.20555936048513765),
+                    "RMSU": np.float64(9.600111233121442),
                     "sign": [np.float64(-1.0)],
-                    "crms": np.float64(0.4285870343218295),
-                    "bias": np.float64(-0.23278371875543458),
-                    "rms": [np.float64(0.4877244157374022)],
-                    "beta_mqi": [np.float64(0.05079940501994467)],
-                    "Hperc": np.float64(0.08138205388790398),
+                    "crms": np.float64(0.043440021333975364),
+                    "rms": [np.float64(0.06827316642055471)],
+                    "beta_mqi": [np.float64(0.0071117057670128615)],
                     "persistence_model": False,
                     "station_type": np.str_("bla"),
                     "UrRV": 0.24,
@@ -111,15 +165,25 @@ def fairmode_stats_example() -> dict:
                     "alpha": 0.2,
                     "freq": "hourly",
                     "percentile": 99.8,
+                    "Np": 5.2,
+                    "Nnp": 5.5,
                 },
                 "Trelew": {
-                    "RMSU": np.float64(9.600004401086988),
+                    "exceedances_obs": 0,
+                    "MPI_mean": np.float64(0.004412910794973912),
+                    "MPI_R_t": np.float64(2.1497603568997187e-06),
+                    "MPI_bias_t": np.float64(7.3362264675617376e-06),
+                    "MPI_std_t": np.float64(0.0001217991618271935),
+                    "MPI_R_s": np.float64(0.001440664233921059),
+                    "MPI_std_s": np.float64(-0.031303579243870915),
+                    "MPI_Hperc": np.float64(-0.0005195246019656068),
+                    "bias": np.float64(7.04278139020833e-05),
+                    "NMB": np.float64(0.0013019203413407127),
+                    "RMSU": np.float64(9.600005426971212),
                     "sign": [np.float64(-1.0)],
-                    "crms": np.float64(0.006734189622538359),
-                    "bias": np.float64(0.010339278114124648),
-                    "rms": [np.float64(0.012338961941489254)],
-                    "beta_mqi": [np.float64(0.0012853079463267891)],
-                    "Hperc": np.float64(-0.0009102268365390729),
+                    "crms": np.float64(0.01412406941853019),
+                    "rms": [np.float64(0.014124245006952786)],
+                    "beta_mqi": [np.float64(0.0014712746898320208)],
                     "persistence_model": False,
                     "station_type": np.str_("bla"),
                     "UrRV": 0.24,
@@ -127,6 +191,8 @@ def fairmode_stats_example() -> dict:
                     "alpha": 0.2,
                     "freq": "hourly",
                     "percentile": 99.8,
+                    "Np": 5.2,
+                    "Nnp": 5.5,
                 },
             }
         }
@@ -156,3 +222,42 @@ def test_save_fairmode_stats(
         / f"{fairmode_exp_output.cfg.proj_id}/{fairmode_exp_output.cfg.exp_id}/fairmode/{list(fairmode_stats_example.keys())[0]}_{obs_name}_{var_name_web}_{vert_code}.json"
     )
     assert fileout.is_file()
+
+
+@pytest.mark.parametrize(
+    "var, val1, val2",
+    [
+        pytest.param(
+            "concno2",
+            24,
+            1,
+            id="no2",
+        ),
+        pytest.param(
+            "concpm10",
+            1,
+            0,
+            id="pm10",
+        ),
+    ],
+)
+def test_exceedances(fairmode_statistics, dummy_coldata_to_fairmode_statistics, var, val1, val2):
+    # reindex fake data hourly and assign new fake values all above threshold (for concno2 threshold is 200)
+    start = dummy_coldata_to_fairmode_statistics.data["time"].values[0]
+    end = dummy_coldata_to_fairmode_statistics.data["time"].values[-1]
+    dummy_coldata_to_fairmode_statistics.data = dummy_coldata_to_fairmode_statistics.data.reindex(
+        {"time": xr.date_range(start, end, freq="h")}
+    )
+
+    nhours = len(xr.date_range(start, end, freq="h"))
+    assert dummy_coldata_to_fairmode_statistics.shape == (2, nhours, 8)
+    dummy_coldata_to_fairmode_statistics.data[1] = dummy_coldata_to_fairmode_statistics.data[
+        1
+    ].where(False, 250.0)
+    dummy_coldata_to_fairmode_statistics.data[0] = dummy_coldata_to_fairmode_statistics.data[
+        0
+    ].where(False, 250.0)
+    [exco, excm] = fairmode_statistics._exceedances(dummy_coldata_to_fairmode_statistics.data, var)
+
+    assert all(exco == nhours // val1 + val2)
+    assert all(excm == nhours // val1 + val2)
diff --git a/tests/cams2_83/test_cams2_83_engine.py b/tests/cams2_83/test_cams2_83_engine.py
index 7d716010..67eb0bb4 100644
--- a/tests/cams2_83/test_cams2_83_engine.py
+++ b/tests/cams2_83/test_cams2_83_engine.py
@@ -46,7 +46,11 @@ def test__calc_forecast_target_MQI_vectorized(eval_config: dict):
     c23engine = CAMS2_83_Engine(setup)
 
     mqi_vectorized = c23engine._calc_forecast_target_MQI_vectorized(
-        example_coldata, example_persistence_coldata, "concno2", 1
+        example_coldata,
+        example_persistence_coldata,
+        "concno2",
+        1,
+        min_num_obs=eval_config["min_num_obs"],
     )
 
     stations_list = example_coldata.coords["station_name"].values.tolist()
diff --git a/tests/io/mscw_ctm/test_reader.py b/tests/io/mscw_ctm/test_reader.py
index 554f71a6..4d7cbaac 100644
--- a/tests/io/mscw_ctm/test_reader.py
+++ b/tests/io/mscw_ctm/test_reader.py
@@ -600,3 +600,10 @@ def test_reader_regexp(tmp_path: Path):
     reader._search_all_files()
 
     assert len(reader._private.filepaths) == 10
+
+
+def test_reader_read_deposition():
+    # https://github.com/metno/pyaerocom/issues/1624
+    file_path = Path("~/MyPyaerocom/testdata-minimal/modeldata/EMEP_DEPOSITION/").expanduser()
+    reader = ReadMscwCtm(data_dir=str(file_path))
+    reader.read_var("dryrdn", ts_type="monthly")
diff --git a/tests/io/test_cachehandler_ungridded.py b/tests/io/test_cachehandler_ungridded.py
index a1291f65..5fd63bac 100644
--- a/tests/io/test_cachehandler_ungridded.py
+++ b/tests/io/test_cachehandler_ungridded.py
@@ -6,6 +6,8 @@ from pyaerocom import UngriddedData
 from pyaerocom.io import ReadAeronetSunV3
 from pyaerocom.io.cachehandler_ungridded import CacheHandlerUngridded
 from tests.conftest import lustre_avail
+from pyaerocom.io.pyaro.pyaro_config import PyaroConfig
+from pyaerocom.io.pyaro.read_pyaro import ReadPyaro
 
 
 @pytest.fixture(scope="module")
@@ -48,3 +50,49 @@ def test_reload(
     reloaded = cache_handler.loaded_data["od550aer"]
     assert isinstance(reloaded, UngriddedData)
     assert reloaded.shape == subset.shape
+
+
+def test_reload_pyaro(cache_handler: CacheHandlerUngridded, pyaro_test_data_file):
+    config0 = PyaroConfig(
+        reader_id="csv_timeseries",
+        filename_or_obj_or_url=str(pyaro_test_data_file),
+        name="test",
+        name_map={
+            "NO": "concno",
+        },
+        post_processing=[
+            "concNno_from_concno",
+        ],
+        filters=dict(),
+    )
+
+    reader0 = ReadPyaro(config0)
+    cache_handler.reader = reader0
+    has_data = cache_handler.check_and_load("concno")
+    assert not has_data, "Data was present"
+
+    data = reader0.read("concno")
+    _outpath = cache_handler.write(data, var_or_file_name="concno")
+
+    has_data = cache_handler.check_and_load("concno")
+    assert has_data, "Data should not be present"
+    assert cache_handler.loaded_data.get("concno") is not None
+
+    reader1 = ReadPyaro(config0)
+    cache_handler.reader = reader1
+    has_data = cache_handler.check_and_load("concno")
+    assert has_data, "Data should be present"
+
+    # Enable a filter which should invalidate the cache
+    config2 = config0.model_copy()
+    config2.filters["countries"] = {"include": "DE"}
+
+    reader2 = ReadPyaro(config2)
+    cache_handler.reader = reader2
+    has_data = cache_handler.check_and_load("concno")
+    assert not has_data, "Data should not be present"
+
+    data = reader2.read("concno")
+    _outpath = cache_handler.write(data, var_or_file_name="concno")
+    has_data = cache_handler.check_and_load("concno")
+    assert has_data, "Data should be present"