Iris regridding spherical vs. ellipsoid

regrid.py

import numpy as np
import iris

spherical_crs = iris.coord_systems.GeogCS(semi_major_axis=6367470)
ellipsoid_crs = iris.coord_systems.GeogCS(semi_major_axis=6378137.0, semi_minor_axis=6356752.314140356)
lambert_azimuthal_crs = iris.coord_systems.LambertAzimuthalEqualArea(
    latitude_of_projection_origin=54.9,
    longitude_of_projection_origin=-2.5,
    false_easting=0.0, false_northing=0.0,
    ellipsoid=ellipsoid_crs)

# defining original data
original_x = np.linspace(-1158000, 924000, 200)
original_y = np.linspace(1036000, 902000, 200)

original_x_coord = iris.coords.DimCoord(
    original_x,
    standard_name='projection_x_coordinate',
    units='m',
    coord_system=lambert_azimuthal_crs)

original_y_coord = iris.coords.DimCoord(
    original_y,
    standard_name='projection_y_coordinate',
    units='m',
    coord_system=lambert_azimuthal_crs)

original_data = np.random.random((200, 200))

original_cube = iris.cube.Cube(
    original_data,
    dim_coords_and_dims=((original_x_coord, 0), (original_y_coord, 1)))

# defining target grid, both on sphere and ellipse.
# (I believe that latitudes on a sphere and latitudes on an ellipse reference slightly different real locations)
new_x = np.linspace(-5, -2, 30)
new_y = new_y = np.linspace(52, 55, 30)

spherical_x = iris.coords.DimCoord(
    new_x,
    standard_name="longitude",
    units='degrees',
    coord_system=spherical_crs)

spherical_y = iris.coords.DimCoord(
    new_y,
    standard_name="latitude",
    units='degrees',
    coord_system=spherical_crs)

new_spherical_cube = iris.cube.Cube(
    np.zeros((30,30)),
    dim_coords_and_dims=((spherical_x, 0), (spherical_y, 1)))

ellipsoid_x = iris.coords.DimCoord(
    new_x,
    standard_name="longitude",
    units='degrees',
    coord_system=ellipsoid_crs)

ellipsoid_y = iris.coords.DimCoord(
    new_y,
    standard_name="latitude",
    units='degrees',
    coord_system=ellipsoid_crs)

new_ellipsoid_cube = iris.cube.Cube(
    np.zeros((30,30)),
    dim_coords_and_dims=((ellipsoid_x, 0), (ellipsoid_y, 1)))


# doing the regrid
spherical_regrid = original_cube.regrid(grid=new_spherical_cube, scheme=iris.analysis.Linear())
ellipsoid_regrid = original_cube.regrid(grid=new_ellipsoid_cube, scheme=iris.analysis.Linear())

# result has different coord_systems
assert spherical_regrid.coord_system() != ellipsoid_regrid.coord_system()

# but same data
assert np.array_equal(spherical_regrid.data, ellipsoid_regrid.data)