Post-regrids Field.cyclic() retains removed domain axes #758
Description
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.