Add obs-only visualization to vertical profiles

pyaerocom/aeroval/helpers.py

@@ -17,14 +17,17 @@
     get_highest_resolution,
     get_max_period_range,
     make_dummy_cube,
-    start_stop_str,
+    # start_stop_str,
+    make_dummy_cube_with_altitude,
 )
 from pyaerocom.io import ReadGridded
 from pyaerocom.tstype import TsType
 from pyaerocom.variable import Variable
 
 logger = logging.getLogger(__name__)
 
+COLOCATION_LAYER_LIMITS_NAME = "colocation_layer_limts"
+
 
 def check_var_ranges_avail(model_data, var_name):
     """
@@ -170,10 +173,15 @@ def make_dummy_model(obs_list: list, cfg) -> str:
     tmp_var_obj = Variable()
     # Loops over variables in obs
     for obs in obs_list:
+        vert_code = cfg.obs_cfg[obs]["obs_vert_type"]
         for var in cfg.obs_cfg[obs]["obs_vars"]:
             # Create dummy cube
-
-            dummy_cube = make_dummy_cube(var, start_yr=start, stop_yr=stop, freq=freq)
+            if hasattr(cfg.obs_cfg[obs], COLOCATION_LAYER_LIMITS_NAME):
+                dummy_cube = make_dummy_cube_with_altitude(
+                    var, start_yr=start, stop_yr=stop, freq=freq
+                )
+            else:
+                dummy_cube = make_dummy_cube(var, start_yr=start, stop_yr=stop, freq=freq)
 
             # Converts cube to GriddedData
             dummy_grid = GriddedData(dummy_cube)

pyaerocom/colocation_3d.py

@@ -336,6 +336,9 @@ def colocate_vertical_profile_gridded(
     except DimensionOrderError:
         data.reorder_dimensions_tseries()
 
+    if not hasattr(data, "altitude"):
+        raise AttributeError(f"GriddedData object must have altitude as a coordinate")
+
     var, var_aerocom = resolve_var_name(data)
     if var_ref is None:
         var_ref = var_aerocom
@@ -372,7 +375,7 @@ def colocate_vertical_profile_gridded(
             )
         else:
             data_ref_meta_idxs_with_var_info.append(i)
-
+    breakpoint()
     if any(
         data.altitude.units != Unit(data_ref.metadata[i]["var_info"]["altitude"]["units"])
         for i in data_ref_meta_idxs_with_var_info

pyaerocom/colocation_auto.py

@@ -1099,6 +1099,7 @@ def _print_processing_status(self):
         oname = self.get_obs_name()
         logger.info(f"Colocation processing status for {mname} vs. {oname}")
         logger.info(self.processing_status)
+        breakpoint()
 
     def _filter_var_matches_var_name(self, var_matches, var_name):
         filtered = {}

pyaerocom/helpers.py

@@ -68,6 +68,7 @@ def varlist_aerocom(varlist):
     elif not isinstance(varlist, list):
         raise ValueError("Need string or list")
     output = []
+
     for var in varlist:
         try:
             _var = const.VARS[var].var_name_aerocom
@@ -1765,3 +1766,102 @@ def make_dummy_cube(
     for coord in dummy.coords():
         coord.points = coord.points.astype(dtype)
     return dummy
+
+
+def make_dummy_cube_with_altitude(
+    var_name: str, start_yr: int = 2000, stop_yr: int = 2020, freq: str = "daily", dtype=float
+) -> iris.cube.Cube:
+    startstr = f"days since {start_yr}-01-01 00:00"
+
+    if freq not in TS_TYPE_TO_PANDAS_FREQ.keys():
+        raise ValueError(f"{freq} not a recognized frequency")
+
+    start_str = f"{start_yr}-01-01 00:00"
+    stop_str = f"{stop_yr}-12-31 00:00"
+    times = pd.date_range(start_str, stop_str, freq=TS_TYPE_TO_PANDAS_FREQ[freq])
+
+    days_since_start = (times - times[0]).days
+    unit = get_variable(var_name).units
+
+    lat_range = (-90, 90)
+    lon_range = (-180, 180)
+    alt_range = (0, 10000)
+    lat_res_deg = 45
+    lon_res_deg = 90
+    alt_res_deg = 2500
+    time_unit = Unit(startstr, calendar="gregorian")
+
+    lons = np.arange(
+        lon_range[0] + (lon_res_deg / 2), lon_range[1] + (lon_res_deg / 2), lon_res_deg
+    )
+    lats = np.arange(
+        lat_range[0] + (lat_res_deg / 2), lat_range[1] + (lat_res_deg / 2), lat_res_deg
+    )
+    alts = np.arange(
+        alt_range[0] + (alt_res_deg / 2), alt_range[1] + (alt_res_deg / 2), alt_res_deg
+    )
+
+    latdim = iris.coords.DimCoord(
+        lats,
+        var_name="lat",
+        standard_name="latitude",
+        long_name="Center coordinates for latitudes",
+        circular=False,
+        units=Unit("degrees"),
+    )
+
+    londim = iris.coords.DimCoord(
+        lons,
+        var_name="lon",
+        standard_name="longitude",
+        long_name="Center coordinates for longitudes",
+        circular=False,
+        units=Unit("degrees"),
+    )
+
+    altdim = iris.coords.DimCoord(
+        alts,
+        var_name="alt",
+        standard_name="altitude",
+        long_name="Center coordinates for altitudes",
+        circular=False,
+        units=Unit("meters"),
+    )
+    timedim = iris.coords.DimCoord(
+        days_since_start, var_name="time", standard_name="time", long_name="Time", units=time_unit
+    )
+
+    latdim.guess_bounds()
+    londim.guess_bounds()
+    dummy = iris.cube.Cube(np.ones((len(times), len(lats), len(lons), len(alts))), units=unit)
+
+    dummy.add_dim_coord(latdim, 1)
+    dummy.add_dim_coord(londim, 2)
+    dummy.add_dim_coord(altdim, 3)
+    dummy.add_dim_coord(timedim, 0)
+    dummy.var_name = var_name
+    dummy.ts_type = freq
+
+    dummy.data = dummy.data.astype(dtype)
+    for coord in dummy.coords():
+        coord.points = coord.points.astype(dtype)
+    return dummy
+
+
+# def add_alt_dim_to_dummy_cube(dummy_cube: iris.cube.Cube) -> iris.cube.Cube:
+#     breakpoint()
+#     alt_range = (0, 10000)
+#     alt_res_deg = 2500
+#     alts = np.arange(
+#         alt_range[0] + (alt_res_deg / 2), alt_range[1] + (alt_res_deg / 2), alt_res_deg
+#     )
+#     altdim = iris.coords.DimCoord(
+#         alts,
+#         var_name="alt",
+#         standard_name="altitude",
+#         long_name="Center coordinates for altitudes",
+#         circular=False,
+#         units=Unit("meters"),
+#     )
+
+#     dummy_cube.add_dim_coord(altdim, 3)