|
| 1 | +import numpy as np |
| 2 | +import xarray as xr |
| 3 | + |
| 4 | +from access_mopper.base import CMIP6_CMORiser |
| 5 | +from access_mopper.derivations import custom_functions, evaluate_expression |
| 6 | + |
| 7 | + |
| 8 | +class CMIP6_Atmosphere_CMORiser(CMIP6_CMORiser): |
| 9 | + """ |
| 10 | + Handles CMORisation of NetCDF datasets using CMIP6 metadata (Atmosphere/Land). |
| 11 | + """ |
| 12 | + |
| 13 | + def select_and_process_variables(self): |
| 14 | + # Find all required bounds variables |
| 15 | + bnds_required = [] |
| 16 | + bounds_rename_map = {} |
| 17 | + for dim, v in self.vocab.axes.items(): |
| 18 | + if v.get("must_have_bounds") == "yes": |
| 19 | + # Find the input dimension name that maps to this output name |
| 20 | + input_dim = None |
| 21 | + for k, val in self.mapping[self.cmor_name]["dimensions"].items(): |
| 22 | + if val == v["out_name"]: |
| 23 | + input_dim = k |
| 24 | + break |
| 25 | + if input_dim is None: |
| 26 | + raise KeyError( |
| 27 | + f"Can't find input dimension mapping for output dimension '{v['out_name']}'." |
| 28 | + ) |
| 29 | + bnds_var = input_dim + "_bnds" |
| 30 | + bounds_rename_map[bnds_var] = v["out_name"] + "_bnds" |
| 31 | + bnds_required.append(bnds_var) |
| 32 | + |
| 33 | + # Select input variables |
| 34 | + input_vars = self.mapping[self.cmor_name]["model_variables"] |
| 35 | + calc = self.mapping[self.cmor_name]["calculation"] |
| 36 | + |
| 37 | + required_vars = set(input_vars + bnds_required) |
| 38 | + self.load_dataset(required_vars=required_vars) |
| 39 | + |
| 40 | + # Handle the calculation type |
| 41 | + if calc["type"] == "direct": |
| 42 | + # If the calculation is direct, just rename the variable |
| 43 | + self.ds = self.ds.rename({input_vars[0]: self.cmor_name}) |
| 44 | + elif calc["type"] == "formula": |
| 45 | + # If the calculation is a formula, evaluate it |
| 46 | + context = {var: self.ds[var] for var in input_vars} |
| 47 | + context.update(custom_functions) |
| 48 | + self.ds[self.cmor_name] = evaluate_expression(calc, context) |
| 49 | + # Drop the original input variables, except the CMOR variable and keep bounds |
| 50 | + self.ds = self.ds.drop_vars( |
| 51 | + [ |
| 52 | + var |
| 53 | + for var in input_vars |
| 54 | + if var != self.cmor_name and var not in bnds_required |
| 55 | + ], |
| 56 | + errors="ignore", |
| 57 | + ) |
| 58 | + else: |
| 59 | + raise ValueError(f"Unsupported calculation type: {calc['type']}") |
| 60 | + |
| 61 | + # Rename dimensions according to the CMOR vocabulary |
| 62 | + dim_rename = self.mapping[self.cmor_name]["dimensions"] |
| 63 | + dims_to_rename = {k: v for k, v in dim_rename.items() if k in self.ds.dims} |
| 64 | + self.ds = self.ds.rename(dims_to_rename) |
| 65 | + |
| 66 | + # Also rename coordinates if needed |
| 67 | + coords_to_rename = {k: v for k, v in dim_rename.items() if k in self.ds.coords} |
| 68 | + if coords_to_rename: |
| 69 | + self.ds = self.ds.rename(coords_to_rename) |
| 70 | + |
| 71 | + # Rename bounds variables |
| 72 | + for bnds_var, out_bnds_name in bounds_rename_map.items(): |
| 73 | + if bnds_var in self.ds: |
| 74 | + self.ds = self.ds.rename({bnds_var: out_bnds_name}) |
| 75 | + elif bnds_var in self.ds.coords: |
| 76 | + self.ds = self.ds.rename({bnds_var: out_bnds_name}) |
| 77 | + |
| 78 | + # Update "bounds" attribute in all variables and coordinates |
| 79 | + for var in list(self.ds.variables) + list(self.ds.coords): |
| 80 | + bounds_attr = self.ds[var].attrs.get("bounds") |
| 81 | + if bounds_attr and bounds_attr in bounds_rename_map: |
| 82 | + self.ds[var].attrs["bounds"] = bounds_rename_map[bounds_attr] |
| 83 | + |
| 84 | + # Transpose the data variable according to the CMOR dimensions |
| 85 | + cmor_dims = self.vocab.variable["dimensions"].split() |
| 86 | + transpose_order = [ |
| 87 | + self.vocab.axes[dim]["out_name"] |
| 88 | + for dim in cmor_dims |
| 89 | + if "value" not in self.vocab.axes[dim] |
| 90 | + ] |
| 91 | + # Squeeze singleton dimensions if they are not in the transpose order |
| 92 | + for dim in self.ds[self.cmor_name].dims: |
| 93 | + if dim not in transpose_order and self.ds[self.cmor_name][dim].size == 1: |
| 94 | + self.ds[self.cmor_name] = self.ds[self.cmor_name].squeeze(dim) |
| 95 | + |
| 96 | + self.ds[self.cmor_name] = self.ds[self.cmor_name].transpose(*transpose_order) |
| 97 | + |
| 98 | + def update_attributes(self): |
| 99 | + self.ds.attrs = { |
| 100 | + k: v |
| 101 | + for k, v in self.vocab.get_required_global_attributes().items() |
| 102 | + if v not in (None, "") |
| 103 | + } |
| 104 | + |
| 105 | + required_coords = { |
| 106 | + v["out_name"] for v in self.vocab.axes.values() if "value" in v |
| 107 | + }.union({v["out_name"] for v in self.vocab.axes.values()}) |
| 108 | + self.ds = self.ds.drop_vars( |
| 109 | + [c for c in self.ds.coords if c not in required_coords], errors="ignore" |
| 110 | + ) |
| 111 | + |
| 112 | + cmor_attrs = self.vocab.variable |
| 113 | + self._check_units(self.cmor_name, cmor_attrs.get("units")) |
| 114 | + |
| 115 | + self.ds[self.cmor_name].attrs.update( |
| 116 | + {k: v for k, v in cmor_attrs.items() if v not in (None, "")} |
| 117 | + ) |
| 118 | + var_type = cmor_attrs.get("type", "double") |
| 119 | + self.ds[self.cmor_name] = self.ds[self.cmor_name].astype( |
| 120 | + self.type_mapping.get(var_type, np.float64) |
| 121 | + ) |
| 122 | + |
| 123 | + try: |
| 124 | + if cmor_attrs.get("valid_min") not in (None, "") and cmor_attrs.get( |
| 125 | + "valid_max" |
| 126 | + ) not in (None, ""): |
| 127 | + vmin = self.type_mapping.get(var_type, np.float64)( |
| 128 | + cmor_attrs["valid_min"] |
| 129 | + ) |
| 130 | + vmax = self.type_mapping.get(var_type, np.float64)( |
| 131 | + cmor_attrs["valid_max"] |
| 132 | + ) |
| 133 | + self._check_range(self.cmor_name, vmin, vmax) |
| 134 | + except ValueError as e: |
| 135 | + raise ValueError( |
| 136 | + f"Failed to validate value range for {self.cmor_name}: {e}" |
| 137 | + ) |
| 138 | + |
| 139 | + for dim, meta in self.vocab.axes.items(): |
| 140 | + name = meta["out_name"] |
| 141 | + dtype = self.type_mapping.get(meta.get("type", "double"), np.float64) |
| 142 | + if name in self.ds: |
| 143 | + self._check_units(name, meta.get("units", "")) |
| 144 | + if meta.get("standard_name") == "time": |
| 145 | + self._check_calendar(name) |
| 146 | + original_units = self.ds[name].attrs.get("units", "") |
| 147 | + coord_attrs = { |
| 148 | + k: v |
| 149 | + for k, v in { |
| 150 | + "standard_name": meta.get("standard_name"), |
| 151 | + "long_name": meta.get("long_name"), |
| 152 | + "units": meta.get("units"), |
| 153 | + "axis": meta.get("axis"), |
| 154 | + "positive": meta.get("positive"), |
| 155 | + "valid_min": dtype(meta["valid_min"]) |
| 156 | + if "valid_min" in meta |
| 157 | + else None, |
| 158 | + "valid_max": dtype(meta["valid_max"]) |
| 159 | + if "valid_max" in meta |
| 160 | + else None, |
| 161 | + }.items() |
| 162 | + if v is not None |
| 163 | + } |
| 164 | + if coord_attrs.get( |
| 165 | + "units" |
| 166 | + ) == "days since ?" and original_units.lower().startswith("days since"): |
| 167 | + coord_attrs["units"] = original_units |
| 168 | + updated = self.ds[name].astype(dtype) |
| 169 | + updated.attrs.update(coord_attrs) |
| 170 | + self.ds[name] = updated |
| 171 | + elif "value" in meta: |
| 172 | + val = meta["value"] |
| 173 | + # Handle character type (e.g., string coordinate) |
| 174 | + if meta["type"] == "character": |
| 175 | + arr = xr.DataArray( |
| 176 | + np.array( |
| 177 | + val, dtype="S" |
| 178 | + ), # ensure type is character (byte string) |
| 179 | + dims=(), |
| 180 | + attrs={ |
| 181 | + k: v |
| 182 | + for k, v in { |
| 183 | + "standard_name": meta.get("standard_name"), |
| 184 | + "long_name": meta.get("long_name"), |
| 185 | + "units": meta.get("units"), |
| 186 | + "axis": meta.get("axis"), |
| 187 | + "positive": meta.get("positive"), |
| 188 | + "valid_min": meta.get("valid_min"), |
| 189 | + "valid_max": meta.get("valid_max"), |
| 190 | + }.items() |
| 191 | + if v is not None |
| 192 | + }, |
| 193 | + ) |
| 194 | + else: |
| 195 | + arr = xr.DataArray( |
| 196 | + dtype(val), |
| 197 | + dims=(), |
| 198 | + attrs={ |
| 199 | + k: v |
| 200 | + for k, v in { |
| 201 | + "standard_name": meta.get("standard_name"), |
| 202 | + "long_name": meta.get("long_name"), |
| 203 | + "units": meta.get("units"), |
| 204 | + "axis": meta.get("axis"), |
| 205 | + "positive": meta.get("positive"), |
| 206 | + "valid_min": dtype(meta["valid_min"]) |
| 207 | + if "valid_min" in meta |
| 208 | + else None, |
| 209 | + "valid_max": dtype(meta["valid_max"]) |
| 210 | + if "valid_max" in meta |
| 211 | + else None, |
| 212 | + }.items() |
| 213 | + if v is not None |
| 214 | + }, |
| 215 | + ) |
| 216 | + self.ds = self.ds.assign_coords({name: arr}) |
0 commit comments