Regridding very similar grids may result in lat/lon coordinates that are not exactly identical

[bouweandela]

Issue reported by @tomaslovato in #485 (comment):

The grid coordinates comparison done in _horizontal_grid_is_close uses the function np.allclose which has default tolerance values set to rtol=1e-05, atol=1e-08 that may not be always enough.
In particular, when lat/lon coordinates are not the same for all regridded datasets any following manipulation made between cubes will produce an iris error from the cube coordinates equality check.

I found this issue while working with CMIP6 ocean data on regular grid from GFDL-CM4and GFDL-ESM4 datasets.

Here below a sample of the latitude field from dissic_Omon_GFDL-CM4_historical_r1i1p1f1_gr_201001-201412.nc

 lat = -89.5, -88.5, -87.5, -86.5, -85.5, -84.5, -83.5, -82.5, -81.5, -80.5,
    -79.5, -78.5, -77.5, -76.5, -75.5, -74.5, -73.5, -72.5, -71.5, -70.5,
    -69.5, -68.5, -67.5, -66.5, -65.5, -64.5, -63.500000000000007105,
    -62.500000000000007105, -61.500000000000007105, -60.5, -59.5, -58.5,
    ...
    63.500000000000007105, 64.5, 65.5, 66.5, 67.5, 68.5, 69.5, 70.5, 71.5,
    72.5, 73.5, 74.5, 75.5, 76.5, 77.5, 78.5, 79.5, 80.5, 81.5, 82.5, 83.5,
    84.5, 85.5, 86.5, 87.5, 88.5, 89.5 ;

I see a couple of possible fixes:

• force the lat/lon coordinates of the cube to be the same as the target_grid at the end of the _regrid function
• in _horizontal_grid_is_close replace the np.allclose with a more stringent iris equality check (as proposed by @valeriupredoi in the above discussion) and modify in here the lat/lon coordinates of the cube.

Originally posted by @tomaslovato in #485 (comment)

[bouweandela]

See #1177 for information on how this was solved for vertical coordinates.

[tomaslovato]

I had a look at the code changes in #1177 and as I understand the solution in there was to adapt the tolerance of the np.allclose function.
As the grid points are recognized to be close enough within the existing function _horizontal_grid_is_close,
could it be a rather simple but effective solution to force the lon/lat values of the source cube with those of the target_grid cube ?
Maybe something like in the following

diff --git a/esmvalcore/preprocessor/_regrid.py b/esmvalcore/preprocessor/_regrid.py
index de822a73a..a39592b37 100644
--- a/esmvalcore/preprocessor/_regrid.py
+++ b/esmvalcore/preprocessor/_regrid.py
@@ -639,6 +639,11 @@ def regrid(cube, target_grid, scheme, lat_offset=True, lon_offset=True):
                     "'%s': %s", original_dtype, cube.core_data().dtype,
                     str(exc))
             cube.data = da.ma.masked_equal(cube.core_data(), fill_value)
+    else:
+        # force target coordinates
+        for coord in ['latitude', 'longitude']:
+            cube.coord(coord).points = target_grid.coord(coord).points
+            cube.coord(coord).bounds = target_grid.coord(coord).bounds

     return cube

[bouweandela]

That looks fine to me. Note that @zklaus requested that the user has control over the tolances in #1177 (comment), so I guess that would also be needed here?

[tomaslovato]

I think that the test performed in _horizontal_grid_is_close is sufficient to determine that the grids a 'similar', so the main point for the coordinate substitution is just to ensure the grid match for the following manipulations ...

if @zklaus agrees on the above solution I'll open a PR to fix it ...

[zklaus]

I like the idea of forcing the target grid on the result. I do think that some user control on the required precision might be useful, particularly in light of projections where the coordinates may take very large values (km in the non-homogeneous projection plane or other fun stuff) in which case much more relaxed requirements may make sense.

But I'd say go ahead with the forcing in a PR and let's discuss fine-tweaking of tolerances there.