Merge pull request #934 from metno/fix930

ReadICOS: CH4 and CO

pyaerocom/aeroval/glob_defaults.py

@@ -99,6 +99,22 @@
         "scale": [400.0, 405.0, 410.0, 415.0, 420.0, 425.0, 430.0, 435.0, 440.0, 445.0, 450.0],
         "colmap": "coolwarm",
     },
+    "vmrch4": {
+        "scale": [
+            1700,
+            1750,
+            1800,
+            1850,
+            1900,
+            1950,
+            2000,
+            2050,
+            2100,
+            2150,
+            2200,
+        ],
+        "colmap": "coolwarm",
+    },
     "concco": {
         "scale": [100.0, 125.0, 150.0, 175.0, 200.0, 225.0, 250.0, 275.0, 300.0],
         "colmap": "coolwarm",
@@ -358,6 +374,7 @@
     concco=["CO", "3D", "Particle concentration"],
     vmrco=["CO", "3D", "Volume mixing ratios"],
     vmrco2=["CO2", "3D", "Volume mixing ratios"],
+    vmrch4=["CH4", "3D", "Volume mixing ratios"],
     # PMs
     concpm10=["PM10", "3D", "Particle concentrations"],
     concpm25=["PM2.5", "3D", "Particle concentrations"],

pyaerocom/colocation_auto.py

@@ -1438,7 +1438,9 @@ def _prepare_colocation_args(self, model_var: str, obs_var: str):
         if self.obs_is_ungridded:
             ts_type = self._get_colocation_ts_type(model_data.ts_type)
             args.update(
-                ts_type=ts_type, var_ref=obs_var, use_climatology_ref=self.obs_use_climatology
+                ts_type=ts_type,
+                var_ref=obs_var,
+                use_climatology_ref=self.obs_use_climatology,
             )
         else:
             ts_type = self._get_colocation_ts_type(model_data.ts_type, obs_data.ts_type)

pyaerocom/data/variables.ini

@@ -3065,6 +3065,17 @@ minimum = 0
 maximum = 1000
 dimensions = time,lat,lon
 
+[vmrch4]
+var_name = vmrch4
+description = CH4 Volume Mixing Ratio
+standard_name = mole_fraction_of_methane_in_air
+var_type = volume mixing ratios
+unit = ppb
+minimum = 0
+maximum = 10000000
+dimensions = time,lat,lon
+
+
 [vmrhno3]
 var_name = vmrhno3
 description = HNO3 Volume Mixing Ratio
@@ -3763,6 +3774,41 @@ maximum = 0.1
 dimensions = time,lat,lon
 comments_and_purpose = Speciation of size partitioning the total aerosol ; Sources attribution and impact analysis
 
+
+[mmrco2]
+var_name = mmrco2
+description = Carbon Dioxide
+standard_name = mass_fraction_of_carbon_dioxide_in_air
+var_type = mass mixing ratio
+unit = kg kg-1
+minimum = 0
+maximum = 1000000
+dimensions = time,lat,lon
+
+[mmrch4]
+var_name = mmrco2
+description = Methane
+standard_name = mass_fraction_of_methane_in_air
+var_type = mass mixing ratio
+unit = kg kg-1
+minimum = 0
+maximum = 1000000
+dimensions = time,lat,lon
+
+[mmrco]
+var_name = mmrco
+description = Carbon Monoxide
+standard_name = mass_fraction_of_carbon_monoxide_in_air
+var_type = mass mixing ratio
+unit = kg kg-1
+minimum = 0
+maximum = 1000000
+dimensions = time,lat,lon
+
+
+
+
+
 [vmrhg0]
 var_name = vmrhg0
 description = Hg0(g) Volume Mixing Ratio

pyaerocom/molmasses.py

@@ -16,6 +16,7 @@
     "glyox": 58.036,
     "hcho": 30.026,
     "co2": 44.0095,
+    "ch4": 16.04,
 }
 
 

pyaerocom/plugins/icos/reader.py

@@ -2,9 +2,18 @@
 
 import logging
 import re
+from collections import defaultdict
+from collections.abc import Iterable
+from functools import cached_property, lru_cache
+from pathlib import Path
+
+import numpy as np
+import xarray as xr
 
 from pyaerocom import const
 from pyaerocom.plugins.mep.reader import ReadMEP
+from pyaerocom.stationdata import StationData
+from pyaerocom.ungriddeddata import UngriddedData
 
 logger = logging.getLogger(__name__)
 
@@ -18,7 +27,7 @@ class ReadICOS(ReadMEP):
     """
 
     #: Mask for identifying datafiles
-    _FILEMASK = "icos-co2-*.nc"
+    _FILEMASK = "icos-*.nc"
 
     #: Version log of this class (for caching)
     __version__ = "0.01"
@@ -49,10 +58,19 @@ class ReadICOS(ReadMEP):
     #: functions used to convert variables that are computed
     AUX_FUNS = {}
 
-    VAR_MAPPING = {"vmrco2": "CO2_volume_mixing_ratio"}
+    VAR_MAPPING = {
+        "vmrco2": "CO2_volume_mixing_ratio",
+        "vmrch4": "CH4_volume_mixing_ratio",
+        "vmrco": "CO_volume_mixing_ratio",
+    }
+
+    FILEMASK_MAPPING = {
+        "vmrco2": "icos-co2-*.nc",
+        "vmrch4": "icos-ch4-*.nc",
+        "vmrco": "icos-co-*.nc",
+    }
 
-    # STATION_REGEX = re.compile("icos-co2-nrt-(.*)-.*-.*-.*.nc")
-    STATION_REGEX = re.compile("icos-co2-(.*)-.*-.*.nc")
+    STATION_REGEX = re.compile("icos-.*-(.*)-.*-.*.nc")  # co2, ch4, co agnostic
 
     DEFAULT_VARS = list(VAR_MAPPING)
 
@@ -65,4 +83,115 @@ def __init__(self, data_id: str | None = None, data_dir: str | None = None):
             data_dir = const.OBSLOCS_UNGRIDDED[const.ICOS_NAME]
 
         super().__init__(data_id=data_id, data_dir=data_dir)
-        self.files = sorted(map(str, self.FOUND_FILES))
+
+    def read_file(
+        self, filename: str | Path, vars_to_retrieve: Iterable[str] | None = None
+    ) -> UngriddedData:
+        """Reads data for a single year for one component"""
+        if not isinstance(filename, Path):
+            filename = Path(filename)
+        if not filename.is_file():
+            raise ValueError(f"missing {filename}")
+        return self.read(vars_to_retrieve, (filename,))
+
+    def read(
+        self,
+        var_to_retrieve: str | None = None,
+        files: Iterable[str | Path] | None = None,
+        first_file: int | None = None,
+        last_file: int | None = None,
+    ) -> UngriddedData:
+        if len(var_to_retrieve) > 1:
+            raise NotImplementedError("Unnskyld, ReadICOS can only read one variable at a time...")
+
+        if var_to_retrieve is None:
+            var_to_retrieve = self.DEFAULT_VARS[0]
+
+        if unsupported := set(var_to_retrieve) - set(self.PROVIDES_VARIABLES):
+            raise ValueError(f"Unsupported variables: {', '.join(sorted(unsupported))}")
+
+        if isinstance(var_to_retrieve, list):
+            var_to_retrieve = var_to_retrieve[0]
+
+        if files is not None and not isinstance(files, tuple):
+            files = tuple(files)
+
+        if files is not None and first_file is not None:
+            files = files[first_file:]
+
+        if files is not None and last_file is not None:
+            files = files[:last_file]
+
+        if files is None:
+            # # Files will depend on the var
+            stations: list[StationData] = []
+            this_var_files = sorted(
+                Path(self.data_dir).rglob(self.FILEMASK_MAPPING[var_to_retrieve])
+            )
+
+            for station_name, paths in self.stations(files).items():
+                paths_to_read = list(set(paths) & set(this_var_files))
+                if not paths_to_read:
+                    continue
+
+                logger.debug(f"Reading station {station_name}")
+                ds = self._read_dataset(paths_to_read)
+                ds = ds.rename({self.VAR_MAPPING[var_to_retrieve]: var_to_retrieve})
+                ds = ds.assign(
+                    time=self._dataset_time(ds),
+                    **{name: func(ds) for name, func in self.AUX_FUNS.items()},
+                )
+                ds.set_coords(("latitude", "longitude", "altitude"))
+                stations.append(self.to_stationdata(ds, station_name))
+            return UngriddedData.from_station_data(stations)
+        else:
+            stations: list[StationData] = []
+            for station_name, paths in self.stations(files).items():
+                logger.debug(f"Reading station {station_name}")
+                ds = self._read_dataset(paths)[var_to_retrieve]
+                ds = ds.rename({self.VAR_MAPPING[var_to_retrieve]: var_to_retrieve})
+                ds = ds.assign(
+                    time=self._dataset_time(ds),
+                    **{name: func(ds) for name, func in self.AUX_FUNS.items()},
+                )
+                ds.set_coords(("latitude", "longitude", "altitude"))
+                stations.append(self.to_stationdata(ds, station_name))
+
+            return UngriddedData.from_station_data(stations)
+
+    def _read_dataset(self, paths: list[Path]) -> xr.Dataset:
+        return xr.open_mfdataset(
+            sorted(paths),
+            concat_dim="time",
+            combine="nested",
+            parallel=True,
+            decode_cf=True,
+        )
+
+    @classmethod
+    def to_stationdata(cls, ds: xr.Dataset, station_name: str) -> StationData:
+        station = StationData()
+        station.data_id = cls.DATA_ID
+        station.dataset_name = cls.DATASET_NAME
+        station.station_id = station_name
+        station.station_name = station_name
+        station.latitude = float(ds["latitude"][0])
+        station.longitude = float(ds["longitude"][0])
+        station.altitude = float(ds["altitude"][0])
+
+        station.station_coords = {
+            "latitude": station.latitude,
+            "longitude": station.longitude,
+            "altitude": station.altitude,
+        }
+
+        station.ts_type = "hourly"
+        station["dtime"] = ds["time"].values
+
+        for var in ds.data_vars:
+            if not var in cls.PROVIDES_VARIABLES:
+                continue
+            station[var] = ds[var].to_series()
+            station["var_info"][var] = {"units": ds[var].units}
+
+        return station

tests/plugins/icos/test_reader.py

@@ -13,7 +13,7 @@
         allow_module_level=True,
     )
 
-VARS_DEFAULT = {"vmrco2"}
+VARS_DEFAULT = {"vmrco2", "vmrch4", "vmrco"}
 VARS_PROVIDED = VARS_DEFAULT  # | {} add more if ever needed
 
 station_names = pytest.mark.parametrize("station", ("bir", "gat", "hpb"))