Post-`regrids` `Field.cyclic()` retains removed domain axes

[sadielbartholomew]

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 https://github.com/NCAS-CMS/cf-python/issues/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.