Aggregate Scores RMSE over multiple cases and produce Hinton plot

[jwarner8]

Something else

Note, this will be taking the average RMSE timeseries over multiple cases. There are different ways to aggregate RMSE. This is mostly to test the plumbing of aggregation, and into a Hinton plot.

[jwarner8]

constraint:
    operator: constraints.generate_var_constraint
    varname:
    - observed_('temperature_at_screen_level',)
    - &id001
      - temperature_at_screen_level
- operator: filters.filter_multiple_cubes
  constraint:
    operator: constraints.generate_cell_methods_constraint
    cell_methods: []

producing weird varname parsing that needs fixing...

[jwarner8]

There are various challenges presented here

1. Parsing of varnames seems to pass a tuple rather than a string, causing issues with filtering.
2. Ensuring we have forecast_reference_time as a promotable dimension in model and obs, currently after reading for two models, for two forecasts, applying fix to above varnames;

0: air_temperature / (K)               (forecast_period: 49; forecast_reference_time: 2; grid_latitude: 342; grid_longitude: 347)
1: observed_temperature_at_screen_level / (K) (time: 49; station: 198)
2: observed_temperature_at_screen_level / (K) (time: 49; station: 198)
3: temperature_at_screen_level / (K)   (time: 48; grid_latitude: 342; grid_longitude: 347)
4: temperature_at_screen_level / (K)   (time: 48; grid_latitude: 342; grid_longitude: 347)

3. Ensuring extract_common_points and interpolating onto point cube can handle multiple models/cube times.

[jwarner8]

We will need an operator that checks that stations are present throughout all forecasts

2026-08-21 14:29:29,928 CSET.operators DEBUG first_arg in kwargs.
2026-08-21 14:29:29,930 CSET.operators.aggregate DEBUG Buckets:
0: observed_temperature_at_screen_level / (K) (time: 49; station: 198)
1: observed_temperature_at_screen_level / (K) (time: 49; station: 198)
failed to merge into a single cube.
  Coordinates in cube.dim_coords differ: station.
  Coordinates in cube.aux_coords (non-scalar) differ: Station_Name, latitude, longitude.

[jwarner8]

note if two case studies are retrieved, the station numbers increase sequentially making it impossible to determine which ones are consistent i.e.

>>> print(cubes_obs[0].coord('station_0').points)
[198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215
 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233
 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251
 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269
 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287
 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305
 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323
 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341
 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359
 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377
 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395]
>>> print(cubes_obs[1].coord('station').points)
[  0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17
  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35
  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53
  54  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71
  72  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89
  90  91  92  93  94  95  96  97  98  99 100 101 102 103 104 105 106 107
 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125
 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143
 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161
 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179
 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197]
>>> 

They do however have unique latitudes and longitudes, which we will need to rebuild a concatenated obs cube.

[jwarner8]

Plan
Takes a CubeList of observation cubes, each corresponding to a different forecast_reference_time.
Uses Station_Name as the station identity.
Finds the intersection of stations present in all cubes.
Drops stations not present in every forecast.
Reorders stations consistently.
Converts time - forecast_reference_time into a forecast_period coordinate.
Stacks the cubes into a single cube with dimensions

forecast_reference_time
forecast_period
station

and keeps time as an auxiliary coordinate on (forecast_reference_time, forecast_period).

[jwarner8]

Have some basic code working for this, with output

2026-08-21 14:54:29,319 CSET.operators INFO Recipe output:
observed_temperature_at_screen_level / (K) (forecast_reference_time: 2; forecast_period: 49; station: 198)
    Dimension coordinates:
        forecast_reference_time                                    x                   -            -
        forecast_period                                            -                   x            -
        station                                                    -                   -            x
    Auxiliary coordinates:
        time                                                       x                   x            -
        Station_Name                                               -                   -            x
        latitude                                                   -                   -            x
        longitude                                                  -                   -            x
    Attributes:
        Conventions                        'CF-1.7'
        cset_comparison_base               1
        model_name                         'OBS'
2026-08-21 14:54:29,320 CSET.operators INFO Creating diagnostic archive.

[jwarner8]

Progress:

Have got

  - operator: misc.combine_cubes_into_cubelist
    first:
      operator: aggregate.combine_obs_across_forecasts
      cubes:
        operator: filters.filter_multiple_cubes
        constraint:
          operator: constraints.generate_var_constraint
          varname: observed_temperature_at_screen_level

    second:     
      operator: aggregate.ensure_aggregatable_across_cases
      cubes:
        operator: filters.filter_multiple_cubes
        constraint:
          operator: constraints.generate_var_constraint
          varname: temperature_at_screen_level

Producing

2026-08-21 15:27:59,651 CSET.operators DEBUG step_input: 0: observed_temperature_at_screen_level / (K) (forecast_reference_time: 2; forecast_period: 49; station: 198)
1: air_temperature / (K)               (forecast_period: 49; forecast_reference_time: 2; grid_latitude: 342; grid_longitude: 347)
2: temperature_at_screen_level / (K)   (forecast_period: 48; forecast_reference_time: 2; grid_latitude: 342; grid_longitude: 347)

Now just to fix some issues with extracting common time points, currently using forecast_period rather than time, as time is now 2D and varies with forecast_period and forecast_reference_time. Some issues in interpolate_to_point_cube with time being 2D, need to investigate.

[jwarner8]

There is a large chunk of code in interpolate_to_point_cube which is doing time checking, not grid checking, and if we are using ensure_common_points operator ahead of this, it is probably redundant

    # Iterate over all cubes and regrid.
    for cube in iter_maybe(fld):
        # Ensure matching times in fld cube and point_cube
        base_time_coord = point_cube.coord("time")
        other_time_coord = cube.coord("time")
        base_times = base_time_coord.units.num2date(base_time_coord.points)
        other_times = other_time_coord.units.num2date(other_time_coord.points)
        shared_times = set.intersection(set(base_times), set(other_times))
        logger.debug("Shared times: %s", shared_times)
        time_constraint = iris.Constraint(
            coord_values={
                "time": lambda cell, shared_times=shared_times: (
                    cell.point in shared_times
                )
            }
        )

        # Extract points matching the shared times.
        cube = cube.extract(time_constraint)
        point_cube = point_cube.extract(time_constraint)
        if cube is None or point_cube is None:
            raise ValueError("No common time points found!")

[jwarner8]

Needs some broader thinking about creating an operator to ensure time is a dimension, with forecast_reference_times also promoted to a time dimension, and forecast_period a 2D aux coord, as forecast_reference_time will always exist.

In meantime is forecast_reference_time is a dim, and there are two or more forecasts, could iterate over and extract common time points.

[jwarner8]

Revisiting this again with fresh perspective; after filtering cubes within the read operator, we are presented;

0: air_temperature / (K)               (forecast_period: 49; forecast_reference_time: 2; grid_latitude: 342; grid_longitude: 347)
1: observed_temperature_at_screen_level / (K) (time: 49; station: 198)
2: observed_temperature_at_screen_level / (K) (time: 49; station: 198)
3: temperature_at_screen_level / (K)   (time: 48; grid_latitude: 342; grid_longitude: 347)
4: temperature_at_screen_level / (K)   (time: 48; grid_latitude: 342; grid_longitude: 347)

Cube 0 automatically merged two forecasts into one. Cubes 3 and 4 can be dealt with by filtering for variable without observed_, and then passing it through ensure_aggregatable_across_cases. We will then end up with;

0: air_temperature / (K)               (forecast_period: 49; forecast_reference_time: 2; grid_latitude: 342; grid_longitude: 347)
1: observed_temperature_at_screen_level / (K) (time: 49; station: 198)
2: observed_temperature_at_screen_level / (K) (time: 49; station: 198)
3: temperature_at_screen_level / (K)   (time: 48; forecast_reference_time: 2; grid_latitude: 342; grid_longitude: 347)

Now, the models have a mix of either time or forecast_period as the lead dimension/auxcoord, which should be easily translatable. For each station, we need to extract all the valid times, and check each model to see if there is a matching time, and if not, drop it. As part of this, the observation cubes can also be merged into a cube that looks like a forecast (i.e. create a forecast_reference_time dimension). We will need to adapt the regrid.interpolate_to_point_cube as this does spatial regridding as well as time matching (two purposes when it should only be doing one).