Ancils for small test domain fail checking during reconfiguration due to longitude rounding errors

[Paul.Gregory]

More fun and games with a small test domain.

Outer domain is specified in the regional nesting suite as

rg01_centre=-28.5,153.5
rg01_igbp_offset=0,0
rg01_name="Lismore"
rg01_nreslns=2
rg01_rs01_delta=0.10,0.10
rg01_rs01_npts=100,100

Inner domain is specified with

rg01_rs02_name="d1100"
rg01_rs02_delta=0.10,0.10
rg01_rs02_npts=50,50

Reconfiguration fails. First there are multiple warnings with the qrclim.sea ancillary file.

????????????????????????????????????????????????????????????????????????????????
??????????????????????????????      WARNING       ??????????????????????????????
?  Warning code: -15000
?  Warning from routine: ANCIL_CHECK_HORIZONTAL_GRID
?  Warning message: Mismatch between model and ancil field x grid spacing
?          Expected x grid spacing : .1000000000
?          Ancil    x grid spacing : .0999998748
?          Ancil file : /scratch/gb02/pag548/cylc-run/rCM3-test-UM-ancil/share/data/ancils/Lismore/d1100/qrclim.sea
?          Lookup num : 1
?          Stashcode  : 96
?  Warning from processor: 0
?  Warning number: 2
????????????????????????????????????????????????????????????????????????????????

There are seven of these warnings, followed by the error

???????????????????????????????????????????????????????????????????????????????
???!!!???!!!???!!!???!!!???!!!       ERROR        ???!!!???!!!???!!!???!!!???!!!
?  Error code: 1001
?  Error from routine: ANCIL_CHECK_MOD::REPORT_ANCIL_ERRORS
?  Error message:
?        ERRORS: 1 ancil files have failed ancil checking and resulted in this abort:
?        -- /scratch/gb02/pag548/cylc-run/rCM3-test-UM-ancil/share/data/ancils/Lismore/d1100/qrclim.sea
?
?        WARNINGS: 1 ancil files have failed ancil checking. These files did not cause an abort
?        due to the setting of l_ignore_ancil_grid_check=.true. in the items namelist:
?        -- /scratch/gb02/pag548/cylc-run/rCM3-test-UM-ancil/share/data/ancils/Lismore/d1100/qrparm.orog.mn
?
?        Any fields causing an error will have produced an ereport
?        warning earlier in the run, please search the log files for
?        each ancil filename for more details on each failure.
?  Error from processor: 0
?  Error number: 14
????????????????????????????????????????????????????????????????????????????????

If you load qrclim.sea into xconv the x coords are not defined at 0.1 degree increments at sufficient precision to please the UM.

Note the y coords of the same file are all displayed in 0.1 degree increments.

For the outer domain (100 by 100 points), both x and y coordinates of qrclim.sea suffer from precision / rounding issues.

For fun, I replaced all entries in app/um/rose-app.conf of
l_ignore_ancil_grid_check=.true.
to
l_ignore_ancil_grid_check=.false.
This made no change, presumably because each entry in rose-app.conf for ancil checking is tied to an ancil filename and STASH entry, and as yet there are none for qrclim.sea

So, given that
a) The same domain size works fine with 0.11 spacing
b) I’m running a large 0.1 degree domain right across most of Australia right now without dramas
c) So has everyone else since the dawn of time
d) I’ve never seen any issues with rounding errors created by the suites in-built ancillaries before
e) There is no extra information in the um_recon pe0 file

I’m a bit stumped. We can manually fix the ancillary with ants but I’m curious as to how this an happen in the first place.

Has anyone ever encountered this before?

[Scott (Scott Wales)]

As I understand it all of the ancillaries should be basing their grid off of the land sea mask. Looking at the RAS app app/ancil_clim_sea/rose-app.conf it’s being given a target land mask. See if you can check if this target mask matches the one you’re using for the run.

[Paul.Gregory]

Yeah, I’ve been digging around the RAS suite for the past few months to understand its workflow and think of ways to build consistent regional UM and MOM6 regional masks.

The main initial task for ancillary generation in RAS is:
<Region>_<resolution>_ancil_lct

which calls the script ancil_lct.py , which in turn returns

1 lct_cube - the 9 level vegetation area fraction cube
2 [land_mask_cube , land_fraction_cube]
where the land_mask_cube is binary, land_fraction_cube is a float (i.e values b/w 0.0 to 1.0)

The land mask cube is then written to disk, which then defines the grid which all other ancillaries are built from.

In our case qrparm.mask (which links to qrparm.mask_cci) displays in xconv as

So I’m going to run the ancil regridding task through a debugger and see why it doesn’t want to respect these grid co-ordinates.

The actual task command is

ants-launch ancil_general_regrid.py --ants-config ${ANTS_CONFIG} \
       ${source} --target-lsm ${target_lsm} --invert-mask -o ${output}

where source=<path to>/GlobColour/v2/qrparm.sea.nc
and target_lsm=<path to>qrparm.mask_cci

[Scott (Scott Wales)]

Yeah, I’ve been digging around the RAS suite for the past few months to understand its workflow and think of ways to build consistent regional UM and MOM6 regional masks.

Start with a single source of truth (for global runs we used the ocean mask as it was a higher resolution) and interpolate the land mask and land fraction from that - land fraction was a conservative regrid of the ocean mask to the land grid with a bit of cleanup something like >=99% land grid points were just set to land and similar for ocean.

There’s an opt in the ancil_lct app vegfrac_only that will get it to use a target mask rather than generate the mask itsself

[Paul.Gregory]

Thanks Scott. I’ve been working with @kieranricardo to replicate their CM3 workflow, which is
a) Define a land mask in MOM6
b) Interpolate that land mask (using the ESMF algorithms used by NUOPC) onto the target UM grid definition.
c) Use that interpolated land mask as the ‘target mask’ for subsequent UM ancil tasks.

I think that reflects your suggestion above. The idea is this retains consistency b/w the NUOPC ‘CAPS’ created inside the NUOPC executable when it passes UM/MOM6 gridded information.

I’ve successfully run ancil_lct in ‘target mask’ mode.

BUT I found the subsequent downstream tasks still required a ASCII grid.nl file (i.e. the namelist definition of the grid) - having a target mask alone is not sufficient.

So now I’m just figuring out which parameters you need to define a grid.nl file.

[Scott (Scott Wales)]

Sounds good - if you’re not using a rotated pole or stretched grid the grid.nl should be pretty basic, pretty much just origin and spacing as I recall - the RAS should generate it

[Paul.Gregory]

So I decided to run through the regridding script ancil_general_regrid.py for this particular case.

After the load_data function has completed we have
a) The target cube (which is the land-sea mask)

target_cube
<iris 'Cube' of land_binary_mask / (1) (latitude: 50; longitude: 50)>

b) The source cubes (which is the monthly chlorophyll in sea water climatology)

source_cubes
[<iris 'Cube' of mass_concentration_of_chlorophyll_a_in_sea_water / (kg m-3) (time: 12; latitude: 481; longitude: 640)>]

The target cube values are read with the correct precision, e.g.

print(target_cube.dim_coords[1][0:10])
DimCoord :  longitude / (degrees)
    points: [
        151. , 151.1, 151.2, 151.3, 151.4, 151.5, 151.6, 151.7, 151.8,
        151.9]
    shape: (10,)
    dtype: float32
    standard_name: 'longitude'
    coord_system: GeogCS(6371229.0)

print(target_cube.dim_coords[0][0:10])
DimCoord :  latitude / (degrees)
    points: [
        -31. , -30.9, -30.8, -30.7, -30.6, -30.5, -30.4, -30.3, -30.2,
        -30.1]
    shape: (10,)
    dtype: float32
    standard_name: 'latitude'
    coord_system: GeogCS(6371229.0)

The source cubes lat/lons are as follows:

print(source_cubes[0].dim_coords[1][:10])
DimCoord :  latitude / (degrees)
    points: [
        -90.   , -89.625, -89.25 , -88.875, -88.5  , -88.125, -87.75 ,
        -87.375, -87.   , -86.625]
    shape: (10,)
    dtype: float32
    standard_name: 'latitude'
    var_name: 'latitude'
    coord_system: GeogCS(6371229.0)


print(source_cubes[0].dim_coords[2][:10])
DimCoord :  longitude / (degrees)
    points: [
        0.    , 0.5625, 1.125 , 1.6875, 2.25  , 2.8125, 3.375 , 3.9375,
        4.5   , 5.0625]
    shape: (10,)
    dtype: float32
    standard_name: 'longitude'
    var_name: 'longitude'
    coord_system: GeogCS(6371229.0)

Let’s now step over the regrid step to see what the regridded_cubes dimensions are:

print(regridded_cubes[0])
mass_concentration_of_chlorophyll_a_in_sea_water / (kg m-3) (time: 12; latitude: 50; longitude: 50)
    Dimension coordinates:
        time                                                     x             -              -
        latitude                                                 -             x              -
        longitude                                                -             -              x
    Cell methods:
        0                                                   time: mean
    Attributes:
        Conventions                                         'CF-1.7'
        STASH                                               m01s00i096
        grid_staggering                                     6
        history                                             '2019-02-08T12:33:57: chlorophyll.py\nhttps://code.metoffice.gov.uk/tra ...'
        source                                              'GLOBCOLOR, http://www.globcolour.info/, see http://fcm1/projects/ANCIL ...'

regridded_cubes[0].dim_coords[1].points
array([-31. , -30.9, -30.8, -30.7, -30.6, -30.5, -30.4, -30.3, -30.2,
       -30.1, -30. , -29.9, -29.8, -29.7, -29.6, -29.5, -29.4, -29.3,
       -29.2, -29.1, -29. , -28.9, -28.8, -28.7, -28.6, -28.5, -28.4,
       -28.3, -28.2, -28.1, -28. , -27.9, -27.8, -27.7, -27.6, -27.5,
       -27.4, -27.3, -27.2, -27.1, -27. , -26.9, -26.8, -26.7, -26.6,
       -26.5, -26.4, -26.3, -26.2, -26.1], dtype=float32)

regridded_cubes[0].dim_coords[2].points
array([151. , 151.1, 151.2, 151.3, 151.4, 151.5, 151.6, 151.7, 151.8,
       151.9, 152. , 152.1, 152.2, 152.3, 152.4, 152.5, 152.6, 152.7,
       152.8, 152.9, 153. , 153.1, 153.2, 153.3, 153.4, 153.5, 153.6,
       153.7, 153.8, 153.9, 154. , 154.1, 154.2, 154.3, 154.4, 154.5,
       154.6, 154.7, 154.8, 154.9, 155. , 155.1, 155.2, 155.3, 155.4,
       155.5, 155.6, 155.7, 155.8, 155.9], dtype=float32)

So far so good - although we note that the bounds of latitude are showing some variation in precision. But not longitude

regridded_cubes[0].dim_coords[1].bounds
array([[-31.05    , -30.95    ],
       [-30.95    , -30.849998],
       [-30.849998, -30.75    ],
       [-30.75    , -30.650002],
       [-30.650002, -30.55    ],
       [-30.55    , -30.45    ],
       [-30.45    , -30.349998],
       [-30.349998, -30.25    ],
       [-30.25    , -30.150002],
       [-30.150002, -30.05    ],
       [-30.05    , -29.95    ],
       [-29.95    , -29.849998],
...

regridded_cubes[0].dim_coords[2].bounds
array([[150.95   , 151.05   ],
       [151.05   , 151.15   ],
       [151.15   , 151.25   ],
       [151.25   , 151.35   ],
       [151.35   , 151.45   ],
       [151.45   , 151.55   ],
       [151.55   , 151.65   ],
       [151.65   , 151.75   ],
       [151.75   , 151.85   ],
       [151.85   , 151.95   ],
       [151.95   , 152.05   ],
       [152.05   , 152.15   ],
       [152.15   , 152.25   ],
       [152.25   , 152.35   ],
       [152.35   , 152.45   ],
       [152.45   , 152.55   ],
       [152.55   , 152.65   ],
       [152.65   , 152.75   ],
       [152.75   , 152.85   ],
       [152.85   , 152.95   ],
....

Because target_lsm_path is defined, we now call ants.analysis.make_consistent_with_lsm.

What is the value of regridded_cubes dimension after this?

print(regridded_cubes[0].dim_coords[1][:10])
DimCoord :  latitude / (degrees)
    points: [
        -31. , -30.9, -30.8, -30.7, -30.6, -30.5, -30.4, -30.3, -30.2,
        -30.1]
    bounds: [
        [-31.05    , -30.95    ],
        [-30.95    , -30.849998],
        ...,
        [-30.25    , -30.150002],
        [-30.150002, -30.05    ]]
    shape: (10,)  bounds(10, 2)
    dtype: float32
    standard_name: 'latitude'
    var_name: 'latitude'
    coord_system: GeogCS(6371229.0)

print(regridded_cubes[0].dim_coords[2][:10])
DimCoord :  longitude / (degrees)
    points: [
        151. , 151.1, 151.2, 151.3, 151.4, 151.5, 151.6, 151.7, 151.8,
        151.9]
    bounds: [
        [150.95, 151.05],
        [151.05, 151.15],
        ...,
        [151.75, 151.85],
        [151.85, 151.95]]
    shape: (10,)  bounds(10, 2)
    dtype: float32
    standard_name: 'longitude'
    var_name: 'longitude'
    coord_system: GeogCS(6371229.0)

At this point the regridded cubes are then saved and the script exits.

Let’s read in the saved data from disk to check it is consistent with what was computed.

a = ants.load(output_path) 

print(a[0].dim_coords[1])
DimCoord :  latitude / (degrees)
    points: [-31. , -30.9, ..., -26.2, -26.1]
    shape: (50,)
    dtype: float32
    standard_name: 'latitude'
    coord_system: GeogCS(6371229.0)
print(a[0].dim_coords[2])
DimCoord :  longitude / (degrees)
    points: [151.     , 151.1    , ..., 155.79999, 155.9    ]
    shape: (50,)
    dtype: float32
    standard_name: 'longitude'
    coord_system: GeogCS(6371229.0)

We can convert this to a data array.

da = xr.DataArray.from_iris(a[0])

In [12]: da.latitude
Out[12]: 
<xarray.DataArray 'latitude' (latitude: 50)> Size: 200B
array([-31. , -30.9, -30.8, -30.7, -30.6, -30.5, -30.4, -30.3, -30.2, -30.1,
       -30. , -29.9, -29.8, -29.7, -29.6, -29.5, -29.4, -29.3, -29.2, -29.1,
       -29. , -28.9, -28.8, -28.7, -28.6, -28.5, -28.4, -28.3, -28.2, -28.1,
       -28. , -27.9, -27.8, -27.7, -27.6, -27.5, -27.4, -27.3, -27.2, -27.1,
       -27. , -26.9, -26.8, -26.7, -26.6, -26.5, -26.4, -26.3, -26.2, -26.1],
      dtype=float32)
Coordinates:
  * latitude  (latitude) float32 200B -31.0 -30.9 -30.8 ... -26.3 -26.2 -26.1
Attributes:
    standard_name:  latitude
    units:          degrees

In [13]: da.longitude
Out[13]: 
<xarray.DataArray 'longitude' (longitude: 50)> Size: 200B
array([151.     , 151.1    , 151.2    , 151.3    , 151.4    , 151.5    ,
       151.6    , 151.7    , 151.8    , 151.9    , 152.     , 152.1    ,
       152.2    , 152.3    , 152.4    , 152.5    , 152.59999, 152.7    ,
       152.8    , 152.9    , 153.     , 153.09999, 153.2    , 153.3    ,
       153.4    , 153.5    , 153.59999, 153.7    , 153.8    , 153.9    ,
       154.     , 154.09999, 154.2    , 154.3    , 154.4    , 154.5    ,
       154.59999, 154.7    , 154.79999, 154.9    , 155.     , 155.09999,
       155.2    , 155.29999, 155.4    , 155.5    , 155.59999, 155.7    ,
       155.79999, 155.9    ], dtype=float32)
Coordinates:
  * longitude  (longitude) float32 200B 151.0 151.1 151.2 ... 155.7 155.8 155.9
Attributes:
    standard_name:  longitude
    units:          degrees

So it seems the regridded latitude values lose precision after being written to disk, or being read from disk.

But can we be sure about the internal representation? These are the full longitude values read from disk

In [16]: print(a[0].dim_coords[2].points)
[151.      151.1     151.2     151.3     151.4     151.5     151.6
 151.7     151.8     151.9     152.      152.1     152.2     152.3
 152.4     152.5     152.59999 152.7     152.8     152.9     153.
 153.09999 153.2     153.3     153.4     153.5     153.59999 153.7
 153.8     153.9     154.      154.09999 154.2     154.3     154.4
 154.5     154.59999 154.7     154.79999 154.9     155.      155.09999
 155.2     155.29999 155.4     155.5     155.59999 155.7     155.79999
 155.9    ]

versus the regridded cube in memory

regridded_cubes[0].dim_coords[2].points
array([151. , 151.1, 151.2, 151.3, 151.4, 151.5, 151.6, 151.7, 151.8,
       151.9, 152. , 152.1, 152.2, 152.3, 152.4, 152.5, 152.6, 152.7,
       152.8, 152.9, 153. , 153.1, 153.2, 153.3, 153.4, 153.5, 153.6,
       153.7, 153.8, 153.9, 154. , 154.1, 154.2, 154.3, 154.4, 154.5,
       154.6, 154.7, 154.8, 154.9, 155. , 155.1, 155.2, 155.3, 155.4,
       155.5, 155.6, 155.7, 155.8, 155.9], dtype=float32)

I’m pretty stumped.

[Scott (Scott Wales)]

This is a um ancil file right? Coordinates are stored as lat[-1], deltalat, lon[-1], deltalon (i.e. origin point in the file is one grid space to the bottom left of the actual grid origin).

If you list the file with mule-pumf these will be lookup fields (UMDP F03)

59. BZY
60. BDY
61. BZX
62. BDX

How do these values compare across the different files?

You might need to check the fields in python directly with mule to check the full precision

import mule
f = mule.AncilFile.from_file("InputFile")
for field in f:
    print(f.bzy) # Etc.

[Scott (Scott Wales)]

Thinking on it though this does ring a bell - maybe coordinates were being computed by ANTS in float32 rather than float64? I might have gotten around this by hacking the ANTS code to change the precision - it’s been a while sorry

[Paul.Gregory]

Scott - no need to apologise

Yes this as a UM ancil file : qrclim.sea which is rejected by the um_recon task at the ANCIL_CHECK_HORIZONTAL_GRID stage due to a mis-match in grid spacing.

I’ll follow up your suggestions and get back to you.

[ashjbarnes (Ashley Barnes)]

@lachlanswhyborn you said you’d seen this issue before right?

[lachlanswhyborn (Lachlan Whyborn)]

I did see this in the past, but I didn’t come up with a good solution- I had to take the grid points from a file that I knew did work, that got the same effective resolution but the one bad point was “fixed” (meaning “what the UM tools expected it to be”).

[lachlanswhyborn (Lachlan Whyborn)]

Did you work out which point was causing the problem? If I remember correctly, for me it was the second last point on the grid, with the UM compatible grid being equally spaced for every point bar the last interval. My suspicion is the issue comes from a F64 to F32 conversion

[Paul.Gregory]

Ok here are some results from using mule to read the spacing values from the fields themselves.

In [46]: ff = mule.AncilFile.from_file('qrparm.mask')

In [47]: first_field = ff.fields[0]

In [48]: bdx = first_field.bdx
...: bdy = first_field.bdy

In [49]: print(f"Field-level DX: **{**bdx**}**")
...: print(f"Field-level DY: **{**bdy**}**")

Field-level DX: 0.10000000000000012
Field-level DY: 0.09999999999999998

In [50]: ff = mule.AncilFile.from_file('qrclim.sea')

In [51]: first_field = ff.fields[0]

In [52]: bdx = first_field.bdx
...: bdy = first_field.bdy

In [53]: print(f"Field-level DX: **{**bdx**}**")
...: print(f"Field-level DY: **{**bdy**}**")

Field-level DX: 0.09999987483024597
Field-level DY: 0.09999999403953552

Which matches the original error caught by ANCIL_CHECK_HORIZONTAL_GRID above

?  Warning message: Mismatch between model and ancil field x grid spacing
?          Expected x grid spacing : .1000000000
?          Ancil    x grid spacing : .0999998748

As an aside - this check is run from line 310 for this file : https://github.com/ACCESS-NRI/UM/blob/dfbb8ba8f583191d66f4cb1c31f081e2627f2073/src/control/ancillaries/ancil_check_mod.F90

      CALL check_lookup_values(lookup_real(bdx, i_lookup),                     &
           grid_space_x_local, i_lookup, i_stash, "x grid spacing", "",        &
           routinename, ancil_file)

which calls this check_lookup_values_real and flags the mis-match here from line 621 onwards

IF (near_equal_real(expected_value, ancil_lookup_value)) THEN
  IF (ancil_file % l_ignore_ancil_grid_check(i_stash) ) THEN
    ! Don't call ereport if logical set to downgrade errors to warnings
    local_error_status = ancil_warning
  ELSE
    ! Always call a warning here but set error status for processing later
    icode = -15000
    local_error_status = ancil_error
    WRITE(cmessageLong, '(A,F0.10,A,F0.10,A)')                                 &
       "Mismatch between model and ancil field " // name_of_check              &
       // TRIM(additional_error_info) // newline //                            &
       "Expected " // name_of_check // " : ", expected_value, newline//        &
       "Ancil    " // name_of_check // " : ", ancil_lookup_value, newline //   &
       TRIM(print_current_ancil_info(ancil_file%filename, i_lookup,            &
       ancil_file%stashcodes(i_stash)))
    CALL ereport(routinename, icode, cmessageLong)
  END IF
ELSE
  local_error_status = ancil_passed
END IF

So what bemuses me is that inside the regridding function itself, all the points look to be represented correctly.

Shall I dig into the ants save functionality? Should we check the values of the bdx field at write-time? And see how this is computed?

The regridding function calls

save.ancil(regridded_cubes, output_path)

The contents of save.ancil are (once you skip all the comments)

   cubes = ants.utils.cube.as_cubelist(cubes)
   ancilfile = _cubes_to_ancilfile(cubes)
   _mule_set_lbuser2(ancilfile)
   ancilfile.to_file(filename)

[Scott (Scott Wales)]

I guess it’s up to you if you fix the saving code or if you just fixup the output file coordinates with mule after it’s created. If you do work out a fix I’m sure the devs would be happy - ANTS is in the public metoffice github space

[Paul.Gregory]

I’ll keep digging into the save functionality and see what I find.

A simple 50x50 domain at 0.1 degree should run out of the box.

[Paul.Gregory]

Ok I think I found the same issue you did.

TLDR; the iris code has a simple switch to determine whether the points along any dimension should be 32 or 64 bit. This is determined in the function regular points in iris/util.py : iris/lib/iris/util.py at main · SciTools/iris · GitHub

def regular_points(zeroth, step, count):
    """Make an array of regular points.

    Create an array of `count` points from `zeroth` + `step`, adding `step` each
    time. In float32 if this gives a sufficiently regular array (tested with
    points_step) and float64 if not.

    Parameters
    ----------
    zeroth : number
        The value *prior* to the first point value.
    step : number
        The numeric difference between successive point values.
    count : number
        The number of point values.

    Notes
    -----
    This function does maintain laziness when called; it doesn't realise data.
    See more at :doc:`/userguide/real_and_lazy_data`.
    """

    def make_steps(dtype: np.dtype):
        start = np.add(zeroth, step, dtype=dtype)
        steps = np.multiply(step, np.arange(count), dtype=dtype)
        return np.add(start, steps, dtype=dtype)

    points = make_steps(np.float32)
    _, regular = iris.util.points_step(points)
    if not regular:
        points = make_steps(np.float64)
    return points

The test regular is computed in function points_step:

def points_step(points):
    """Determine whether `points` has a regular step.

    Parameters
    ----------
    points : numeric, array-like
        The sequence of values to check for a regular difference.

    Returns
    -------
    numeric, bool
        A tuple containing the average difference between values, and whether
        the difference is regular.

    Notes
    -----
    This function does not maintain laziness when called; it realises data.
    See more at :doc:`/userguide/real_and_lazy_data`.
    """
    # Calculations only make sense with multiple points
    points = np.asanyarray(points)
    if points.size >= 2:
        diffs = np.diff(points)
        avdiff = np.mean(diffs)
        # TODO: This value for `rtol` is set for test_analysis to pass...
        regular = np.allclose(diffs, avdiff, rtol=0.001)
    else:
        avdiff = np.nan
        regular = True
    return avdiff, regular

For our simple test domain (50x50 at 0.1 spacing) - the regular test passes and hence the points used to define downstream ancillary grids are cast as 32 bit floats. The grid spacing computed from these 32 bit floats doesn’t match the spacing from the original land sea mask, and hence you have errors in the UM reconfiguration step.

If you increase the relative tolerance to 0.00001, the regular test fails and the grid points appear to be `exact’ inside the python representation.

Full breakdown:

The mask coordinates are built in the converter function.

def _make_cube(field, converter):
    # Convert the field to a Cube.
    metadata = converter(field)
    cube_data = field.core_data()
    cube = iris.cube.Cube(
        cube_data,
        attributes=metadata.attributes,
        cell_methods=metadata.cell_methods,
        dim_coords_and_dims=metadata.dim_coords_and_dims,
        aux_coords_and_dims=metadata.aux_coords_and_dims,
    )

At this stage we have read the following fields from the binary mask file:

field.bdx
np.float64(0.10000000000000012)
field.bzx
np.float64(150.8)
field.lbnpt
np.int64(50)

The cube’s dim_coords_and_dims are built in this section of function iris.fileformats.pp_load_rules.convert:

def convert(f):
    """Convert a PP field into the corresponding items of Cube metadata.

    Parameters
    ----------
    f : :class:`iris.fileformats.pp.PPField`

    Returns
    -------
    :class:`iris.fileformats.rules.ConversionMetadata` object.

which calls the function _all_other_rules:

    # All the other rules.
    (
        references,
        standard_name,
        long_name,
        units,
        attributes,
        cell_methods,
        dim_coords_and_dims,
        other_aux_coords_and_dims,
    ) = _all_other_rules(f)

For this example, inside _all_other_rules we use the method DimCoord.from_regular():

    # "Normal" (i.e. not cross-sectional) lats+lons (--> vector coordinates)
    if f.bdx != 0.0 and f.bdx != f.bmdi and len(f.lbcode) != 5 and f.lbcode[0] == 1:
        dim_coords_and_dims.append(
            (
                DimCoord.from_regular(
                    f.bzx,
                    f.bdx,
                    f.lbnpt,
                    standard_name=f._x_coord_name(),
                    units="degrees",
                    circular=(f.lbhem in [0, 4]),
                    coord_system=f.coord_system(),
                ),
                1,
            )
        )

We call this method with the inputs

f.bzx
np.float64(150.8)
f.bdx
np.float64(0.10000000000000012)
f.lbnpt
np.int64(50)
f.coord_system()
GeogCS(6371229.0)
f.lbhem
np.int64(3)

The class method from_regular defined in iris.coords contains the function_regular_points: iris/lib/iris/util.py at main · SciTools/iris · GitHub

def regular_points(zeroth, step, count):
    """Make an array of regular points.

    Create an array of `count` points from `zeroth` + `step`, adding `step` each
    time. In float32 if this gives a sufficiently regular array (tested with
    points_step) and float64 if not.

    Parameters
    ----------
    zeroth : number
        The value *prior* to the first point value.
    step : number
        The numeric difference between successive point values.
    count : number
        The number of point values.

    Notes
    -----
    This function does maintain laziness when called; it doesn't realise data.
    See more at :doc:`/userguide/real_and_lazy_data`.
    """

    def make_steps(dtype: np.dtype):
        start = np.add(zeroth, step, dtype=dtype)
        steps = np.multiply(step, np.arange(count), dtype=dtype)
        return np.add(start, steps, dtype=dtype)

    points = make_steps(np.float32)
    _, regular = iris.util.points_step(points)
    if not regular:
        points = make_steps(np.float64)
    return points

For our test case, these values are:

zeroth
np.float64(150.8)
step
np.float64(0.10000000000000012)
count
np.int64(50)

The default values of points computed in this function is:

points
array([150.90001, 151.00002, 151.1 , 151.20001, 151.3 , 151.40001, 151.50002, 151.6 , 151.70001, 151.8 , 151.90001, 152.00002, 152.1 , 152.20001, 152.3 , 152.40001, 152.50002, 152.6 , 152.70001, 152.8 , 152.90001, 153.00002, 153.1 , 153.20001, 153.3 , 153.40001, 153.50002, 153.6 , 153.70001, 153.8 , 153.90001, 154.00002, 154.1 , 154.20001, 154.3 , 154.40001, 154.50002, 154.6 , 154.70001, 154.8 , 154.90001, 155.00002, 155.1 , 155.20001, 155.3 , 155.40001, 155.50002, 155.6 , 155.70001, 155.8 ], dtype=float32)

because the values of diffs are:

array([0.1000061 , 0.09999084, 0.1000061 , 0.09999084, 0.1000061 , 0.1000061 , 0.09999084, 0.1000061 , 0.09999084, 0.1000061 , 0.1000061 , 0.09999084, 0.1000061 , 0.09999084, 0.1000061 , 0.1000061 , 0.09999084, 0.1000061 , 0.09999084, 0.1000061 , 0.1000061 , 0.09999084, 0.1000061 , 0.09999084, 0.1000061 , 0.1000061 , 0.09999084, 0.1000061 , 0.09999084, 0.1000061 , 0.1000061 , 0.09999084, 0.1000061 , 0.09999084, 0.1000061 , 0.1000061 , 0.09999084, 0.1000061 , 0.09999084, 0.1000061 , 0.1000061 , 0.09999084, 0.1000061 , 0.09999084, 0.1000061 , 0.1000061 , 0.09999084, 0.1000061 , 0.09999084], dtype=float32)

avdiff is equal

avdiff
np.float32(0.099999875)

So the default test

regular = np.allclose(diffs, avdiff, rtol=0.001)

passes in this case, and so we have 32-bit representations of the points which loses a lot of precision for our value bdx.

If you change the above to

regular = np.allclose(diffs, avdiff, rtol=0.00001)

The points are returned as 64-bit floats:

In [3]: mask.dim_coords[1].points
Out[3]: 
array([150.9, 151. , 151.1, 151.2, 151.3, 151.4, 151.5, 151.6, 151.7,
       151.8, 151.9, 152. , 152.1, 152.2, 152.3, 152.4, 152.5, 152.6,
       152.7, 152.8, 152.9, 153. , 153.1, 153.2, 153.3, 153.4, 153.5,
       153.6, 153.7, 153.8, 153.9, 154. , 154.1, 154.2, 154.3, 154.4,
       154.5, 154.6, 154.7, 154.8, 154.9, 155. , 155.1, 155.2, 155.3,
       155.4, 155.5, 155.6, 155.7, 155.8])

In [4]: mask.dim_coords[1].dtype
Out[4]: dtype('float64')

So, given that ants.io.load.load is a wrapper for the iris load function, what’s the best way to increase the tolerance of the iris load function in ants?

I’m thinking about writing my own iris loader using the functions listed above, but enforcing 64-bit representation of the points.

Scott - if the above has jogged your memory, can you remember how you worked around this change the precision?

[Paul.Gregory]

Anyway, in the interim I wrote a couple of lines of python to replace the original 32-bit lat/lon points with new 64-bit precision points using the values ancillary header files.