From 440c62be82b239095ceadaa037d241aea4f9843b Mon Sep 17 00:00:00 2001 From: lewisblake Date: Tue, 14 Nov 2023 08:48:52 +0000 Subject: [PATCH] WIP for now. revisit when have time --- pyaerocom/aeroval/helpers.py | 14 +++-- pyaerocom/colocation_3d.py | 5 +- pyaerocom/colocation_auto.py | 1 + pyaerocom/helpers.py | 100 +++++++++++++++++++++++++++++++++++ 4 files changed, 116 insertions(+), 4 deletions(-) diff --git a/pyaerocom/aeroval/helpers.py b/pyaerocom/aeroval/helpers.py index 5f192ac8c..72ce8ec46 100644 --- a/pyaerocom/aeroval/helpers.py +++ b/pyaerocom/aeroval/helpers.py @@ -17,7 +17,8 @@ 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 @@ -25,6 +26,8 @@ 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) diff --git a/pyaerocom/colocation_3d.py b/pyaerocom/colocation_3d.py index 82a0b96b4..4efd80341 100644 --- a/pyaerocom/colocation_3d.py +++ b/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 diff --git a/pyaerocom/colocation_auto.py b/pyaerocom/colocation_auto.py index 9d9f45658..765b09691 100644 --- a/pyaerocom/colocation_auto.py +++ b/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 = {} diff --git a/pyaerocom/helpers.py b/pyaerocom/helpers.py index 858a754e2..a3ebd648d 100644 --- a/pyaerocom/helpers.py +++ b/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)