Collocation of emep-mscw output with gridded data in a projection

pyaerocom/colocation/colocation_utils.py

@@ -681,10 +681,12 @@
     """
     if filter_name is None:
         filter_name = const.DEFAULT_REG_FILTER
-    try:
-        data.check_dimcoords_tseries()
-    except DimensionOrderError:
-        data.reorder_dimensions_tseries()
+    if data.proj_info is None:
+        # reordering only implemented if no projection
+        try:
+            data.check_dimcoords_tseries()
+        except DimensionOrderError:
+            data.reorder_dimensions_tseries()
 
     var, var_aerocom = resolve_var_name(data)
     if var_ref is None:
@@ -747,12 +749,28 @@
 
     latitude = data.latitude.points
     longitude = data.longitude.points
-    lat_range = [np.min(latitude), np.max(latitude)]
-    lon_range = [np.min(longitude), np.max(longitude)]
-    # use only sites that are within model domain
+    if data.proj_info is None:
+        lat_range = [np.min(latitude), np.max(latitude)]
+        lon_range = [np.min(longitude), np.max(longitude)]
+        # use only sites that are within model domain
 
-    # filter_by_meta wipes is_vertical_profile
-    data_ref = data_ref.filter_by_meta(latitude=lat_range, longitude=lon_range)
+        # filter_by_meta wipes is_vertical_profile
+        data_ref = data_ref.filter_by_meta(latitude=lat_range, longitude=lon_range)
+    else:
+        # gridded data with projection,
+        # add x/y information to ungridded
+        for coord in data.cube.dim_coords:
+            if coord.var_name == data.proj_info.x_axis:
+                vals = coord.points
+                xrange = (np.min(vals), np.max(vals))
+            if coord.var_name == data.proj_info.y_axis:
+                vals = coord.points
+                yrange = (np.min(vals), np.max(vals))
+        if xrange is None or yrange is None:
+            raise VariableDefinitionError(
+                f"x/y axis not found in cube: {data.proj_info.x_axis}, {data.proj_info.y_axis}"
+            )
+        data_ref = data_ref.filter_by_projection(data.proj_info.to_proj, xrange, yrange)
 
     # get timeseries from all stations in provided time resolution
     # (time resampling is done below in main loop)
@@ -774,6 +792,7 @@
             f"Variable {var_ref} is not available in specified time interval ({start}-{stop})"
         )
 
+    # to read
     grid_stat_data = data.to_time_series(longitude=ungridded_lons, latitude=ungridded_lats)
 
     pd_freq = col_tst.to_pandas_freq()

pyaerocom/griddeddata.py

@@ -41,6 +41,7 @@
 )
 from pyaerocom.helpers_landsea_masks import load_region_mask_iris
 from pyaerocom.mathutils import estimate_value_range, exponent
+from pyaerocom.projection_information import ProjectionInformation
 from pyaerocom.region import Region
 from pyaerocom.stationdata import StationData
 from pyaerocom.time_config import IRIS_AGGREGATORS, TS_TYPE_TO_NUMPY_FREQ
@@ -130,10 +131,17 @@
         concatenated=False,
         region=None,
         reader=None,
+        proj_info=None,
     )
 
     def __init__(
-        self, input=None, var_name=None, check_unit=True, convert_unit_on_init=True, **meta
+        self,
+        input=None,
+        var_name=None,
+        check_unit=True,
+        convert_unit_on_init=True,
+        proj_info: ProjectionInformation | None = None,
+        **meta,
     ):
         if input is None:
             input = iris.cube.Cube([])
@@ -152,6 +160,8 @@
         self._coord_var_names = None
         self._coord_standard_names = None
         self._coord_long_names = None
+        # projection information
+        meta.update({"proj_info": proj_info})
 
         if input:
             self.load_input(input, var_name)
@@ -195,6 +205,10 @@
                 f"No default access available for variable {self.var_name}"
             )
 
+    @property
+    def proj_info(self) -> ProjectionInformation:
+        return self.metadata["proj_info"]
+
     @property
     def ts_type(self):
         """
@@ -1130,10 +1144,11 @@
             collapse_scalar = coords.pop("collapse_scalar")
         else:
             collapse_scalar = True
-        try:
-            self.check_dimcoords_tseries()
-        except DimensionOrderError:
-            self.reorder_dimensions_tseries()
+        if self.proj_info is None:
+            try:
+                self.check_dimcoords_tseries()
+            except DimensionOrderError:
+                self.reorder_dimensions_tseries()
         pinfo = False
         if np.prod(self.shape) > 5913000:  # (shape of 2x2 deg, daily data)
             pinfo = True
@@ -1176,10 +1191,11 @@
         )
 
     def _to_time_series_xarray(self, scheme="nearest", add_meta=None, ts_type=None, **coords):
-        try:
-            self.check_dimcoords_tseries()
-        except DimensionOrderError:
-            self.reorder_dimensions_tseries()
+        if self.proj_info is None:
+            try:
+                self.check_dimcoords_tseries()
+            except DimensionOrderError:
+                self.reorder_dimensions_tseries()
 
         arr = self.to_xarray()
 
@@ -1198,7 +1214,21 @@
                 lon = vals
         if lat is None or lon is None:
             raise ValueError("Please provide latitude and longitude coords")
-        subset = extract_latlon_dataarray(arr, lat, lon, method=scheme, new_index_name="latlon")
+        if self.proj_info is None:
+            subset = extract_latlon_dataarray(
+                arr, lat, lon, method=scheme, new_index_name="latlon"
+            )
+        else:
+            x, y = self.proj_info.to_proj(lat, lon)
+            subset = extract_latlon_dataarray(
+                arr,
+                lon=x,
+                lat=y,
+                lon_dimname=self.proj_info.x_axis,
+                lat_dimname=self.proj_info.y_axis,
+                method=scheme,
+                new_index_name="latlon",
+            )
 
         lat_id = subset.attrs["lat_dimname"]
         lon_id = subset.attrs["lon_dimname"]
@@ -1219,8 +1249,11 @@
         result = []
         subset = subset.compute()
         data_np = subset.data
-        lats = subset[lat_id].data
-        lons = subset[lon_id].data
+        if self.proj_info is None:
+            lats = subset[lat_id].data
+            lons = subset[lon_id].data
+        else:
+            lats, lons = self.proj_info.to_latlon(subset[lon_id].data, subset[lat_id].data)
         for sidx in range(subset.shape[-1]):
             data = StationData(
                 latitude=lats[sidx],
@@ -2603,8 +2636,8 @@
                     f"Skipping assignment of var_name from metadata in GriddedData, "
                     f"since attr. needs to be str and is {val}"
                 )
-                continue
-            self._grid.attributes[key] = val
+            else:
+                self._grid.attributes[key] = val
 
     def delete_all_coords(self, inplace=True):
         """Deletes all coordinates (dimension + auxiliary) in this object"""

pyaerocom/helpers.py

@@ -119,6 +119,7 @@ def extract_latlon_dataarray(
     check_domain=True,
 ):
     """Extract individual lat / lon coordinates from `DataArray`
+    lon/lat can also be x/y coordinates if the `DataArray` has only projected axes.
 
     Parameters
     ----------

pyaerocom/io/mscw_ctm/reader.py

@@ -11,6 +11,7 @@
 from pyaerocom import const
 from pyaerocom.exceptions import VarNotAvailableError
 from pyaerocom.griddeddata import GriddedData
+from pyaerocom.projection_information import ProjectionInformation
 from pyaerocom.units_helpers import UALIASES
 
 from .additional_variables import (
@@ -599,19 +600,21 @@ def _compute_var(self, var_name_aerocom, ts_type):
 
         Returns
         -------
-        GriddedData
+        xarray.DataArray
             loaded data object
+
         """
 
         temp_arrs = []
         req = self.AUX_REQUIRES[var_name_aerocom]
         aux_func = self.AUX_FUNS[var_name_aerocom]
         logger.info(f"computing {var_name_aerocom} from {req} using {aux_func}")
+        proj_info = None
         for aux_var in self.AUX_REQUIRES[var_name_aerocom]:
-            arr = self._load_var(aux_var, ts_type)
+            arr, proj_info = self._load_var(aux_var, ts_type)
             temp_arrs.append(arr)
 
-        return aux_func(*temp_arrs)
+        return aux_func(*temp_arrs), proj_info
 
     def _load_var(self, var_name_aerocom, ts_type):
         """
@@ -636,6 +639,8 @@ def _load_var(self, var_name_aerocom, ts_type):
         -------
         xarray.DataArray
             loaded data
+        ProjectionInformation
+            projection of variable
 
         """
         if var_name_aerocom in self.var_map:  # can be read
@@ -677,7 +682,7 @@ def read_var(self, var_name, ts_type=None, **kwargs):
 
         ts_type = self.ts_type
 
-        arr = self._load_var(var_name_aerocom, ts_type)
+        arr, proj_info = self._load_var(var_name_aerocom, ts_type)
         if arr.units in UALIASES:
             arr.attrs["units"] = UALIASES[arr.units]
         try:
@@ -693,6 +698,7 @@ def read_var(self, var_name, ts_type=None, **kwargs):
             ts_type=ts_type,
             check_unit=True,
             convert_unit_on_init=True,
+            proj_info=proj_info,
         )
 
         #!obsolete
@@ -732,14 +738,16 @@ def _read_var_from_file(self, var_name_aerocom, ts_type):
         -------
         xarray.DataArray
             loaded data
+        ProjectionInformation
+            projection of variable
 
         """
         emep_var = self.var_map[var_name_aerocom]
 
         try:
             filedata = self.filedata
             data = filedata[emep_var]
-
+            proj_info = ProjectionInformation.from_xarray(filedata, emep_var)
         except KeyError:
             raise VarNotAvailableError(
                 f"{var_name_aerocom} ({emep_var}) not available in {self.filename}"
@@ -749,7 +757,7 @@ def _read_var_from_file(self, var_name_aerocom, ts_type):
         data.time.attrs["standard_name"] = "time"
         prefix = emep_var.split("_")[0]
         data.attrs["units"] = self.preprocess_units(data.units, prefix)
-        return data
+        return data, proj_info
 
     @staticmethod
     def preprocess_units(units, prefix):

pyaerocom/projection_information.py

@@ -0,0 +1,89 @@
+import xarray
+from pyproj import CRS, transformer
+
+
+class ProjectionInformationException(Exception):
+    pass
+
+
+class ProjectionInformation:
+    def __init__(self):
+        self._crs = None
+        self._x_axis = None
+        self._y_axis = None
+        self._units = None
+
+    @property
+    def x_axis(self):
+        return self._x_axis
+
+    @property
+    def y_axis(self):
+        return self._y_axis
+
+    def to_proj(self, latitude, longitude):
+        """convert latitude and longitude to x and y
+
+        :param latitude: latitude values
+        :param longitude: longitude values
+        :return: tuple of x and y coordinates
+        """
+        trans = transformer.Transformer.from_crs(4326, self._crs)
+        return trans.transform(latitude, longitude)
+
+    def to_latlon(self, x, y):
+        """convert x and y to latitude and longitude
+
+        :param x: x values
+        :param y: y values
+        :return: tuple of latitude and longitude coordinates
+        """
+        trans = transformer.Transformer.from_crs(self._crs, 4326)
+        return trans.transform(x, y)
+
+    @staticmethod
+    def from_xarray(ds: xarray.Dataset, var: str):
+        """initialize the ProjectionInformation from an xarray variable
+
+        returns None if no projection exists for the variable or ProjectionInformation
+        """
+        da = ds[var]
+        if not "grid_mapping" in da.attrs:
+            return None
+        pi = ProjectionInformation()
+        pi._crs = CRS.from_cf(ds[da.grid_mapping].attrs)
+
+        for c in da.coords:
+            if len(da.coords[c].dims) != 1:
+                continue
+            if "standard_name" in da.coords[c].attrs:
+                if da.coords[c].standard_name in (
+                    "longitude",
+                    "grid_longitude",
+                    "projection_x_coordinate",
+                ):
+                    pi._x_axis = c
+                    break
+            if "axis" in da.coords[c].attrs:
+                if da.coords[c].axis in ("x", "X"):
+                    pi._x_axis = c
+                    break
+        for c in da.coords:
+            if len(da.coords[c].dims) != 1:
+                continue
+            if "standard_name" in da.coords[c].attrs:
+                if da.coords[c].standard_name in (
+                    "latitude",
+                    "grid_latitude",
+                    "projection_y_coordinate",
+                ):
+                    pi._y_axis = c
+                    break
+            if "axis" in da.coords[c].attrs:
+                if da.coords[c].axis in ("y", "Y"):
+                    pi._y_axis = c
+                    break
+        if pi._x_axis is None or pi._y_axis is None:
+            raise ProjectionInformationException(f"no x or y axis found for variable '{var}'")
+        pi._units = da.coords[pi._y_axis].units
+        return pi

pyaerocom/sample_data_access/minimal_dataset.py

@@ -8,7 +8,7 @@
 logger = logging.getLogger(__name__)
 
 #: tarfile to download
-TESTDATA_FILE = "testdata-minimal.tar.gz.20231116"
+TESTDATA_FILE = "testdata-minimal.tar.gz.20240722"
 
 minimal_dataset = pooch.create(
     path=const.OUTPUTDIR,  # ~/MyPyaerocom/
@@ -21,6 +21,7 @@
         "testdata-minimal.tar.gz.20231017": "md5:705d91e01ca7647b4c93dfe67def661f",
         "testdata-minimal.tar.gz.20231019": "md5:f8912ee83d6749fb2a9b1eda1d664ca2",
         "testdata-minimal.tar.gz.20231116": "md5:5da747f6596817295ba7affe3402b722",
+        "testdata-minimal.tar.gz.20240722": "md5:7d933901c6d273d012f132c60df086cc",
     },
 )