Post-`regrids` `Field.cyclic()` retains removed domain axes
On the current main branch, a Field after a regrids can return with Field.cyclic() a (cyclic) domain axes which no longer exists on the Field because it has been removed via the regridding operation. I believe this is the underlying reason that we see the issue of #756, though I think adding in a tweak to the logic there as suggested in #757 makes it more robust, so think it would be good to add regardless and will leave that PR open for consideration.
We may want to prioritise a fix for this to put in before our imminent release, because it is probably a issue with the new feature of "Added ... the option to regrid the vertical axis in logarithmic coordinates to cf.Field.regrids and cf.Field.regridc " (since it shows up in a vertical regrid, as per the details below).
Example case and behaviour
I've spent a lot of time trying to get a minimal (or at least, pared down) working example of this to share here, but I am struggling to find something so have nearly given up on a small reproducer, in favour of moving onto getting a fix since I feel like I have isolated the problem enough to do that now.
But, for detail and context, the issue has emerged as part of the breaking of our VISION end-to-end script recently as noticed during #756 (see the report there for the underlying error emerging). So I'll share relevant details of that here now which can hopefully shine sufficient light on the problem, since I am unsure whether I'll have time before going on leave after tomorrow to get the reduced example up.
Before the regridding we have a field a with print(a, a.constructs(), a.cyclic()) being (I added spaces between lines for clarity):
Field: id%UM_m01s51i010_vn1105 (ncvar%UM_m01s51i010_vn1105) ----------------------------------------------------------- Data : id%UM_m01s51i010_vn1105(time(5), air_pressure(19), latitude(144), longitude(192)) 1 Cell methods : time(5): point Dimension coords: time(5) = [2017-07-03 11:00:00, ..., 2017-07-03 15:00:00] gregorian : air_pressure(19) = [1000.0, ..., 100.0] hPa : latitude(144) = [-89.375, ..., 89.375] degrees_north : longitude(192) = [0.9375, ..., 359.0625] degrees_east Constructs: {'cellmethod0': <CF CellMethod: domainaxis0: point>, 'dimensioncoordinate0': <CF DimensionCoordinate: time(5) days since 2017-1-1 gregorian>, 'dimensioncoordinate1': <CF DimensionCoordinate: air_pressure(19) hPa>, 'dimensioncoordinate2': <CF DimensionCoordinate: latitude(144) degrees_north>, 'dimensioncoordinate3': <CF DimensionCoordinate: longitude(192) degrees_east>, 'domainaxis0': <CF DomainAxis: size(5)>, 'domainaxis1': <CF DomainAxis: size(19)>, 'domainaxis2': <CF DomainAxis: size(144)>, 'domainaxis3': <CF DomainAxis: size(192)>} {'domainaxis3'}
and with a field b of:
Field: mole_fraction_of_ozone_in_air (ncvar%O3_TECO) ---------------------------------------------------- Data : mole_fraction_of_ozone_in_air(ncdim%obs(11160)) ppb Auxiliary coords: time(ncdim%obs(11160)) = [2017-07-03 11:15:07, ..., 2017-07-03 14:21:06] standard : altitude(ncdim%obs(11160)) = [2577.927001953125, ..., 151.16905212402344] m : air_pressure(ncdim%obs(11160)) = [751.6758422851562, ..., 1006.53076171875] hPa : latitude(ncdim%obs(11160)) = [52.56147766113281, ..., 52.0729866027832] degree_north : longitude(ncdim%obs(11160)) = [0.3171832859516144, ..., -0.6249311566352844] degree_east : cf_role=trajectory_id(cf_role=trajectory_id(1)) = [STANCO]
after an operation of:
c = a.regrids(b, method="linear", z="air_pressure", ln_z=True,)
we get, for print(c, c.constructs(), c.cyclic()):
AFTER Field: id%UM_m01s51i010_vn1105 (ncvar%UM_m01s51i010_vn1105) ----------------------------------------------------------- Data : id%UM_m01s51i010_vn1105(time(5), ncdim%obs(11160)) 1 Cell methods : time(5): point Dimension coords: time(5) = [2017-07-03 11:00:00, ..., 2017-07-03 15:00:00] gregorian Auxiliary coords: time(ncdim%obs(11160)) = [2017-07-03 11:15:07, ..., 2017-07-03 14:21:06] standard : altitude(ncdim%obs(11160)) = [2577.927001953125, ..., 151.16905212402344] m : air_pressure(ncdim%obs(11160)) = [751.6758422851562, ..., 1006.53076171875] hPa : latitude(ncdim%obs(11160)) = [52.56147766113281, ..., 52.0729866027832] degree_north : longitude(ncdim%obs(11160)) = [0.3171832859516144, ..., -0.6249311566352844] degree_east Constructs: {'auxiliarycoordinate0': <CF AuxiliaryCoordinate: time(11160) days since 1900-01-01 00:00:00 standard>, 'auxiliarycoordinate1': <CF AuxiliaryCoordinate: altitude(11160) m>, 'auxiliarycoordinate2': <CF AuxiliaryCoordinate: air_pressure(11160) hPa>, 'auxiliarycoordinate3': <CF AuxiliaryCoordinate: latitude(11160) degree_north>, 'auxiliarycoordinate4': <CF AuxiliaryCoordinate: longitude(11160) degree_east>, 'cellmethod0': <CF CellMethod: domainaxis0: point>, 'dimensioncoordinate0': <CF DimensionCoordinate: time(5) days since 2017-1-1 gregorian>, 'domainaxis0': <CF DomainAxis: size(5)>, 'domainaxis4': <CF DomainAxis: size(11160)>} {'domainaxis3'}
Notice how 'domainaxis3' is no longer a construct on the Field c, but it remains listed in c.cyclic(), when it should have been removed.