4D model data in colocator. requires preprocessing

pyaerocom/colocation_auto.py

@@ -1418,8 +1418,10 @@ def _prepare_colocation_args(self, model_var: str, obs_var: str):
             ts_type = self._get_colocation_ts_type(model_data.ts_type, obs_data.ts_type)
             args.update(ts_type=ts_type)
         if self.obs_is_vertical_profile:
-            args.update(colocation_layer_limits=self.colocation_layer_limits,
-                        profile_layer_limits=self.profile_layer_limits)
+            args.update(
+                colocation_layer_limits=self.colocation_layer_limits,
+                profile_layer_limits=self.profile_layer_limits,
+            )
         return args
 
     def _check_dimensionality(self, args):
@@ -1431,10 +1433,6 @@ def _check_dimensionality(self, args):
         if mdata.ndim == 4 and self.obs_vert_type == "Surface":
             mdata = mdata.extract_surface_level()
             args["data"] = mdata
-        elif mdata.ndim > 3:
-            raise DataDimensionError(
-                f"cannot co-locate model data with more than 3 dimensions: {mdata}"
-            )
 
         if isinstance(odata, GriddedData):
             if odata.ndim == 4 and self.obs_vert_type == "Surface":

pyaerocom/griddeddata.py

@@ -135,7 +135,6 @@ class GriddedData:
     def __init__(
         self, input=None, var_name=None, check_unit=True, convert_unit_on_init=True, **meta
     ):
-
         if input is None:
             input = iris.cube.Cube([])
 
@@ -742,7 +741,6 @@ def _check_invalid_unit_alias(self):
         """
         cube = self.grid
         if "invalid_units" in cube.attributes and cube.attributes["invalid_units"] in UALIASES:
-
             from_unit = cube.attributes["invalid_units"]
             to_unit = UALIASES[from_unit]
             logger.info(f"Updating invalid unit in {repr(cube)} from {from_unit} to {to_unit}")
@@ -825,7 +823,6 @@ def _try_convert_custom_unit(self, new_unit):
         self._apply_unit_mulfac(new_unit, mulfac)
 
     def _apply_unit_mulfac(self, new_unit, mulfac):
-
         if mulfac != 1:
             new_cube = self._grid * mulfac
             new_cube.attributes.update(self._grid.attributes)
@@ -1042,7 +1039,6 @@ def mean_at_coords(self, latitude=None, longitude=None, time_resample_kwargs=Non
         return np.nanmean(mean)
 
     def _coords_to_iris_sample_points(self, **coords):
-
         sample_points = []
         num = None
         for cname, vals in coords.items():
@@ -1057,7 +1053,7 @@ def _coords_to_iris_sample_points(self, **coords):
 
     def _iris_sample_points_to_coords(self, sample_points):
         lats, lons = None, None
-        for (name, vals) in sample_points:
+        for name, vals in sample_points:
             if isnumeric(vals):
                 vals = [vals]
             if name in ("lat", "latitude"):
@@ -1082,7 +1078,6 @@ def to_time_series(
         use_iris=False,
         **coords,
     ):
-
         """Extract time-series for provided input coordinates (lon, lat)
 
         Extract time series for each lon / lat coordinate in this cube or at
@@ -1176,7 +1171,6 @@ def to_time_series(
         )
 
     def _to_time_series_xarray(self, scheme="nearest", add_meta=None, ts_type=None, **coords):
-
         try:
             self.check_dimcoords_tseries()
         except DimensionOrderError:
@@ -1223,7 +1217,6 @@ def _to_time_series_xarray(self, scheme="nearest", add_meta=None, ts_type=None,
         lats = subset[lat_id].data
         lons = subset[lon_id].data
         for sidx in range(subset.shape[-1]):
-
             data = StationData(
                 latitude=lats[sidx],
                 longitude=lons[sidx],
@@ -1320,7 +1313,6 @@ def _to_timeseries_2D(
     def _to_timeseries_3D(
         self, sample_points, scheme, collapse_scalar, vert_scheme, add_meta=None
     ):
-
         # Data contains vertical dimension
         data = self._apply_vert_scheme(sample_points, vert_scheme)
 
@@ -1413,8 +1405,7 @@ def _infer_index_surface_level(self):
                 return np.argmin(self.grid.dim_coords[3].points)
             elif coord.attributes["positive"] == "down":
                 return np.argmax(self.grid.dim_coords[3].points)
-
-        breakpoint()
+
         try:
             coord = vc.VerticalCoordinate(cname)
             if coord.lev_increases_with_alt:
@@ -1647,15 +1638,16 @@ def _resample_time_iris(self, to_ts_type):
         return data
 
     def _resample_time_xarray(self, to_ts_type, how, min_num_obs):
-
         arr = xr.DataArray.from_iris(self.cube)
         from_ts_type = self.ts_type
         try:
             rs = TimeResampler(arr)
             arr_out = rs.resample(
                 to_ts_type, from_ts_type=from_ts_type, how=how, min_num_obs=min_num_obs
             )
-        except ValueError:  # likely non-standard datetime objects in array (cf https://github.com/pydata/xarray/issues/3426)
+        except (
+            ValueError
+        ):  # likely non-standard datetime objects in array (cf https://github.com/pydata/xarray/issues/3426)
             arr["time"] = self.time_stamps()
             rs = TimeResampler(arr)
             arr_out = rs.resample(
@@ -2088,7 +2080,6 @@ def _check_meta_netcdf(self):
         self.cube.attributes = meta_out
 
     def _to_netcdf_aerocom(self, out_dir, **kwargs):
-
         years = self.years_avail()
         outpaths = []
         for subset in self.split_years(years):
@@ -2563,7 +2554,6 @@ def delete_aux_vars(self):
         """Delete auxiliary variables and iris AuxFactories"""
         c = self.cube
         for aux_fac in c.aux_factories:
-
             c.remove_aux_factory(aux_fac)
 
         for coord in c.coords():

pyaerocom/io/iris_io.py

@@ -75,6 +75,7 @@ def load_cubes_custom(files, var_name=None, file_convention=None, perform_fmt_ch
     """
     cubes = []
     loaded_files = []
+
     for i, _file in enumerate(files):
         try:
             cube = load_cube_custom(
@@ -122,6 +123,7 @@ def load_cube_custom(file, var_name=None, file_convention=None, perform_fmt_chec
         file = str(file)  # iris load does not like PosixPath
     if perform_fmt_checks is None:
         perform_fmt_checks = const.GRID_IO.PERFORM_FMT_CHECKS
+
     cube_list = iris.load(file)
     cube = None
     if var_name is None:
@@ -213,7 +215,6 @@ def check_and_regrid_lons_cube(cube):
 
 
 def check_dim_coord_names_cube(cube):
-
     from pyaerocom import const
 
     coords = dict(
@@ -536,7 +537,6 @@ def correct_time_coord(cube, ts_type, year):
 
 
 def _check_correct_dtypes_timedim_cube_list(cubes):
-
     try:
         dtypes = np.unique([cube.coord("time").points.dtype for cube in cubes])
     except iris.exceptions.CoordinateNotFoundError:

pyaerocom/vert_coords.py

@@ -185,7 +185,6 @@ def is_supported(standard_name):
 
 
 class VerticalCoordinate:
-
     NAMES_SUPPORTED = {
         "altitude": "z",
         "air_pressure": "pres",
@@ -216,7 +215,9 @@ class VerticalCoordinate:
     FUNS_YIELD = {"asc": "air_pressure", "ahspc": "air_pressure", "gph": "altitude"}
 
     _LEV_INCREASES_WITH_ALT = dict(
-        atmosphere_sigma_coordinate=False, atmosphere_hybrid_sigma_pressure_coordinate=False
+        atmosphere_sigma_coordinate=False,
+        atmosphere_hybrid_sigma_pressure_coordinate=False,
+        altitude=True,
     )
 
     def __init__(self, name=None):
@@ -287,7 +288,6 @@ def calc_pressure(self, lev, **kwargs):
         """
 
         if not self.var_name in self.NAMES_SUPPORTED:
-
             raise CoordinateNameError(
                 f"Variable {self.var_name} cannot be converted to pressure levels. "
                 f"Conversion is only possible for supported variables:\n{self.vars_supported_str}"
@@ -319,7 +319,6 @@ def pressure2altitude(self, p, **kwargs):
 
 
 class AltitudeAccess:
-
     #: Additional variable names (in AEROCOM convention) that are used
     #: to search for additional files that can be used to access or compute
     #: the altitude levels at each grid point
@@ -479,7 +478,6 @@ def _check_vars_in_data_obj(self):
 
     # ToDo: check alias names
     def _check_var_in_data_obj(self, var_name):
-
         c = VerticalCoordinate(var_name)
 
         if c.var_name in self.data_obj:
@@ -515,7 +513,6 @@ def check_altitude_access(self, **coord_info):
         return False
 
     def _check_altitude_access_helper(self, coord_name, **coord_info):
-
         cstd_name = const.COORDINFO[coord_name].standard_name
 
         if not self.search_aux_coords(coord_name):