Error when reading file having spherical ellipsoid and inverse_flattening=0. GeogCS case (1 0 0)

[larsbarring]

🐛 Bug Report

Ellipsoid parameters inverse_flattening = 0 should be accepted without warning.

How To Reproduce

Steps to reproduce the behaviour:

1. Grab a suitable file form esgf: e.g. a CORDEX run with the ICTP RegCM4-6 model:
   tas_EUR-11_MPI-M-MPI-ESM-LR_historical_r1i1p1_ICTP-RegCM4-6_v1_mon_197001-198012
2. Check relevant metadata:

ncdump -h tas_EUR-11_MPI-M-MPI-ESM-LR_historical_r1i1p1_ICTP-RegCM4-6_v1_mon_197001-198012.nc
netcdf tas_EUR-11_MPI-M-MPI-ESM-LR_historical_r1i1p1_ICTP-RegCM4-6_v1_mon_197001-198012 {
dimensions:
y = 527 ;
x = 527 ;
time = UNLIMITED ; // (132 currently)
bnds = 2 ;
variables:
char crs ;
crs:proj4_params = "+proj=lcc +lat_1=30.00 +lat_2=65.00 +lat_0=48.00 +lon_0=9.75 +x_0=-6000. +y_0=-6000. +ellps=sphere +a=6371229. +b=6371229. +units=m +no_defs" ;
crs:CoordinateTransformType = "Projection" ;
crs:CoordinateAxisTypes = "GeoX GeoY" ;
crs:semi_major_axis = 6371229. ;
crs:inverse_flattening = 0. ;
crs:false_easting = -6000. ;
crs:false_northing = -6000. ;
crs:grid_mapping_name = "lambert_conformal_conic" ;
crs:standard_parallel = 30., 65. ;
crs:longitude_of_central_meridian = 9.75 ;
crs:latitude_of_projection_origin = 48. ;

3. Load the file

 >python
Python 3.8.1 | packaged by conda-forge | (default, Jan 29 2020, 14:55:04) 
[GCC 7.3.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import iris
>>> print(iris.__version__)
2.4.0
>>> mydatadir = ".../cordex/output/EUR-11/ICTP/MPI-M-MPI-ESM-LR/historical/r1i1p1/RegCM4-6/v1/mon/tas/v20190502/"
>>> testfile = "tas_EUR-11_MPI-M-MPI-ESM-LR_historical_r1i1p1_ICTP-RegCM4-6_v1_mon_197001-198012.nc"
>>> c = iris.load(mydatadir + testfile)
.../.conda/envs/climix-devel-3/lib/python3.8/site-packages/iris/coord_systems.py:180: RuntimeWarning: divide by zero encountered in double_scalars
  semi_minor_axis = semi_major_axis - ((1.0 / inverse_flattening)

Expected behaviour

Iris should accept inverse_flattening = 0. as valid input for the case of a spherical ellipsoid. After all this is one way to define a spherical ellipsoid, the other one being semi_major_axis = semi_minor_axis, which iris handles.

Environment

• OS & Version: CentOS Linux 7 (Core)
• Iris Version: 2.4.0

Additional context

Looking at the GeogCS code in master branch it seems like only a small change is needed:
line 197-202)

    elif semi_major_axis is not None and (
        semi_minor_axis is None and (
              inverse_flattening is None or
              inverse_flattening == 0)
    ):
        semi_minor_axis = semi_major_axis
        inverse_flattening = 0.0  

[larsbarring]

In fact, although iris issues a warning, it should definitely be considered an error because the resulting ellipsoid parameters make no sense:

print(c.coord("projection_x_coordinate").coord_system.ellipsoid)                                                                                                                                                                      
GeogCS(semi_major_axis=6371229.0, semi_minor_axis=-inf)

[jamesp]

Hi @larsbarring, thanks for both the bug report, and the potential fix!
Have you committed to iris or another github project before? If not I can show you how to create a PR and submit the proposed fix.

I wonder if

    elif semi_major_axis is not None and (
        semi_minor_axis is None and not inverse_flattening
    ):
        semi_minor_axis = semi_major_axis
        inverse_flattening = 0.0  

might be simpler, but we should check if it breaks anything in the test suite

[larsbarring]

Hi @jamesp
Basically I am not a software guy/developer, so --
-- Only pretty basic commits and PRs, but not at the iris level of sophistication in terms of structured requirements and infrastructure support/automation.
-- Your simplification makes sense.
-- Yes, but I have not managed to wrap my head around what unit tests mean and how they are done (other than they are supposed to test the code ....)