Derived coordinates are not loaded anymore?

[valeriupredoi]

📰 Custom Issue

Hey guys, one more for me today, I am using this file and loading it simple:

import iris

cs = iris.load("/home/valeriu/ESMValCore/tests/integration/cmor/_fixes/test_data/common_cl_ap.nc")
for c in cs:
    print(c)

with iris=3 I am not seeing the derived coordinates no more:

• with iris=2.4:

b_bnds / (1)                        (atmosphere_hybrid_sigma_pressure_coordinate: 2; -- : 2)
     Dimension coordinates:
          atmosphere_hybrid_sigma_pressure_coordinate                           x       -
     Auxiliary coordinates:
          ap                                                                    x       -
          b                                                                     x       -
     Derived coordinates:
          air_pressure                                                          x       -
ap_bnds / (1)                       (atmosphere_hybrid_sigma_pressure_coordinate: 2; -- : 2)
     Dimension coordinates:
          atmosphere_hybrid_sigma_pressure_coordinate                           x       -
     Auxiliary coordinates:
          ap                                                                    x       -
          b                                                                     x       -
     Derived coordinates:
          air_pressure                                                          x       -
cloud_area_fraction_in_atmosphere_layer    (time: 1; atmosphere_hybrid_sigma_pressure_coordinate: 2; latitude: 3; longitude: 4)
     Dimension coordinates:
          time                                            x                                               -            -             -
          atmosphere_hybrid_sigma_pressure_coordinate     -                                               x            -             -
          latitude                                        -                                               -            x             -
          longitude                                       -                                               -            -             x
     Auxiliary coordinates:
          surface_air_pressure                            x                                               -            x             x
          ap                                              -                                               x            -             -
          b                                               -                                               x            -             -
     Derived coordinates:
          air_pressure                                    x                                               x            x             x
surface_air_pressure / (Pa)         (time: 1; latitude: 3; longitude: 4)
     Dimension coordinates:
          time                           x            -             -
          latitude                       -            x             -
          longitude                      -            -             x
     Attributes:
          additional_attribute: xyz

• with iris=3.0:

b_bnds / (unknown)                  (atmosphere_hybrid_sigma_pressure_coordinate: 2; -- : 2)
     Dimension coordinates:
          atmosphere_hybrid_sigma_pressure_coordinate                           x       -
     Auxiliary coordinates:
          ap                                                                    x       -
          b                                                                     x       -
ap_bnds / (unknown)                 (atmosphere_hybrid_sigma_pressure_coordinate: 2; -- : 2)
     Dimension coordinates:
          atmosphere_hybrid_sigma_pressure_coordinate                           x       -
     Auxiliary coordinates:
          ap                                                                    x       -
          b                                                                     x       -
cloud_area_fraction_in_atmosphere_layer       (time: 1; atmosphere_hybrid_sigma_pressure_coordinate: 2; latitude: 3; longitude: 4)
     Dimension coordinates:
          time                                            x                                               -            -             -
          atmosphere_hybrid_sigma_pressure_coordinate     -                                               x            -             -
          latitude                                        -                                               -            x             -
          longitude                                       -                                               -            -             x
     Auxiliary coordinates:
          surface_air_pressure                            x                                               -            x             x
          ap                                              -                                               x            -             -
          b                                               -                                               x            -             -
surface_air_pressure / (Pa)         (time: 1; latitude: 3; longitude: 4)
     Dimension coordinates:
          time                           x            -             -
          latitude                       -            x             -
          longitude                      -            -             x
     Attributes:
          additional_attribute: xyz

where did the Derived coordinates go?

[bjlittle]

@valeriupredoi Ohh very interesting, I'll take a peek for you...

[bjlittle]

@valeriupredoi Okay, so this took a wee bit of digging, but the mystery is thankfully solved 🕵️

Firstly, I can recreate your issue in iris 3.0.0, and I can also confirm that iris 2.4.0 does indeed realise derived coordinates for your example dataset 👍

This change in behaviour between iris 3.0.0 and 2.4.0 directly relates to PR #3795, where the units of an AuxCood or DimCoord that is not explicitly specified defaults to unknown rather then 1 (dimensionless).

We thought about this change in behaviour long and hard, as the treatment of default units was inconsistent across iris for different objects, and the CF Metadata standard did not add any clarity (no surprises there, indeed it was inconsistent), but in the end we unified the behaviour for consistency sake. That is, erring on the side that it correct and safer to assign units as unknown in the case when they are not specified, whereas it is potentially incorrect and not safe to assume default dimensionless units instead.

This change in behaviour is mentioned in the iris 3.0.0 release notes in the Incompatible Changes section, fifth bullet point.

Where this is stinging you is for the missing units for the sigma variable associated with the atmosphere_hybrid_sigma_pressure_coordinate i.e., given formula_terms = "ap: ap_bnds b: b_bnds ps: ps", sigma in this case is the NetCDF b variable, and this is defined within the example file simply as double b(lev) ;

As the b variable has no units associated with it iris 3.0.0 assigns the default units of unknown (which they are). However, during the creation of the auxiliary factory that creates the derived coordinate it enforces the contract that the sigma variable must be dimensionless, which in this case it's not, as it's unknown. Note that, enforcing that sigma is dimensionless has not changed from 2.4.0 to 3.0.0.

This failure for sigma units results in iris raising a warning and not creating the associated HybridPressureFactory auxiliary factory i.e., during loading you see the following warning:

UserWarning: Invalid units: sigma must be dimensionless.

Also, note the example given in the CF Metadata Conventions, where their sigma variable has explicit units = "1".

So in order to fix this, my advise is to add the correct units to your sigma variable in the file. In addition to this, you could also add the missing bounds attribute for the sigma (b) and delta (ap) variables e.g.,

netcdf common_cl_ap {
dimensions:
	time = 1 ;
	lev = 2 ;
	lat = 3 ;
	lon = 4 ;
	bnds = 2 ;
variables:
	double time(time) ;
		time:standard_name = "time" ;
		time:units = "days since 6543-2-1" ;
	double lev(lev) ;
		lev:bounds = "lev_bnds" ;
		lev:standard_name = "atmosphere_hybrid_sigma_pressure_coordinate" ;
		lev:units = "1" ;
		lev:formula_terms = "ap: ap b: b ps: ps" ;
	double lev_bnds(lev, bnds) ;
		lev_bnds:standard_name = "atmosphere_hybrid_sigma_pressure_coordinate" ;
		lev_bnds:units = "1" ;
		lev_bnds:formula_terms = "ap: ap_bnds b: b_bnds ps: ps" ;
	double lat(lat) ;
		lat:standard_name = "latitude" ;
		lat:units = "degrees_north" ;
	double lon(lon) ;
		lon:standard_name = "longitude" ;
		lon:units = "degrees_east" ;
	double b(lev) ;
                b:units = "1" ;                                                 <== add "units" for sigma
                b:bounds = "b_bnds" ;                                           <== add "bounds" for sigma
	double b_bnds(lev, bnds) ;
	double ps(time, lat, lon) ;
		ps:standard_name = "surface_air_pressure" ;
		ps:units = "Pa" ;
		ps:additional_attribute = "xyz" ;
	float cl(time, lev, lat, lon) ;
		cl:standard_name = "cloud_area_fraction_in_atmosphere_layer" ;
		cl:units = "%" ;
	double ap(lev) ;
		ap:units = "Pa" ;
                ap:bounds = "ap_bnds" ;                                         <== add "bounds" for delta
	double ap_bnds(lev, bnds) ;
}

Hope this helps 😄

[bjlittle]

The alternative is that we relax the contract in the auxiliary factories to accept coordinates that are dimensionless or unknown. However that could cause other undesirable side effects or complications e.g., The computation of the derived coordinate would break due to the use of sigma with unknown units (I think) or at least you'd end up with a derived coordinate also with units of unknown.

However, iris could override the unknown units to be dimensionless. This may affect load/save round-tripping... that would need to be investigated. Regardless, this would be an iris 3.0.1 release

[bjlittle]

@valeriupredoi Done a quick hack to see, if in principle, this promotion of units of unknown to dimensionless (1) within auxiliary factories works...

Python 3.7.8 | packaged by conda-forge | (default, Nov 17 2020, 23:45:15) 
Type 'copyright', 'credits' or 'license' for more information
IPython 7.12.0 -- An enhanced Interactive Python. Type '?' for help.
PyDev console: using IPython 7.12.0
Python 3.7.8 | packaged by conda-forge | (default, Nov 17 2020, 23:45:15) 
[GCC 7.5.0] on linux
>>> import iris
>>> iris.__version__
'3.0.0'
>>> fname = "/downloads/common_cl_ap.nc"
>>> cubes = iris.load(fname)
/net/home/h05/itwl/projects/git/iris/lib/iris/fileformats/cf.py:1207: UserWarning: Ignoring formula terms variable 'ps' referenced by data variable 'ap_bnds' via variable 'lev': Dimensions ('time', 'lat', 'lon') do not span ('lev', 'bnds')
  warnings.warn(msg)
/net/home/h05/itwl/projects/git/iris/lib/iris/fileformats/cf.py:1207: UserWarning: Ignoring formula terms variable 'ps' referenced by data variable 'b_bnds' via variable 'lev': Dimensions ('time', 'lat', 'lon') do not span ('lev', 'bnds')
  warnings.warn(msg)
/net/home/h05/itwl/projects/git/iris/lib/iris/fileformats/netcdf.py:685: UserWarning: Unable to find coordinate for variable 'ps'
  "{!r}".format(name)
/net/home/h05/itwl/projects/git/iris/lib/iris/fileformats/netcdf.py:685: UserWarning: Unable to find coordinate for variable 'ps'
  "{!r}".format(name)
>>> print(cubes)
0: ap_bnds / (unknown)                 (atmosphere_hybrid_sigma_pressure_coordinate: 2; -- : 2)
1: b_bnds / (unknown)                  (atmosphere_hybrid_sigma_pressure_coordinate: 2; -- : 2)
2: cloud_area_fraction_in_atmosphere_layer / (%) (time: 1; atmosphere_hybrid_sigma_pressure_coordinate: 2; latitude: 3; longitude: 4)
3: surface_air_pressure / (Pa)         (time: 1; latitude: 3; longitude: 4)
print(cubes[2])
cloud_area_fraction_in_atmosphere_layer    (time: 1; atmosphere_hybrid_sigma_pressure_coordinate: 2; latitude: 3; longitude: 4)
     Dimension coordinates:
          time                                            x                                               -            -             -
          atmosphere_hybrid_sigma_pressure_coordinate     -                                               x            -             -
          latitude                                        -                                               -            x             -
          longitude                                       -                                               -            -             x
     Auxiliary coordinates:
          surface_air_pressure                            x                                               -            x             x
          ap                                              -                                               x            -             -
          b                                               -                                               x            -             -
     Derived coordinates:
          air_pressure                                    x                                               x            x             x

seems possible

[bjlittle]

@bouweandela and @valeriupredoi How much of a deal breaker is this for your guys? I'm assuming an iris 3.0.1 would be preferable?

[valeriupredoi]

@bjlittle cheers for the super quick and detailed answers, mate! We appreciate very much your very quick response times to issues we raise 🍺 About this, I am afraid it generates a chicken-and-egg problem - the chicken being faulty model data that is already on the ESGF and the egg is us needing to fix it via our CMOR fixes before we put it through the analysis. The file I linked above is a replication of such a faulty model data file that we use in a test of a fix: if we can't load the derived coordinates it means we can't fix them in ESMValCore and that means we will not be able to run ESMValTool on it, well, at least without not turning off all the CMOR checks that, in turn, will generate a whole lot of other headaches. So yeah, I think this is something we kind of really need 😁

[bjlittle]

@valeriupredoi Okay, no worries dude. Leave it with me and I'll get back to you asap... just need to speak to the team about this 👍

[valeriupredoi]

@valeriupredoi Okay, no worries dude. Leave it with me and I'll get back to you asap... just need to speak to the team about this +1

you're legend @bjlittle cheers muchly! 🍺

[bjlittle]

@valeriupredoi Okay, I'm dropping everything today to service this - I can totally appreciate the situation you're in, so I'm going to push an iris 3.0.1 release today.

Hopefully, I can turn this around quickly for you, as I know you're working to a tight deadline.

I'll ping you here @valeriupredoi with my progress on this to keep you in the loop 👍