Editing cube coordinate datum....

[nickschurch]

Hi,

I'm new to iris so this might be a stupid Q, but here goes. I'm truing to regrid some data from a grib2 file onto the same grid as a netCDF file. When I load the data from the netCDF (target) file I get:

<stdin>:1: FutureWarning: Ignoring a datum in netCDF load for consistency with existing behaviour. In a future version of Iris, this datum will be applied. To apply the datum when loading, use the iris.FUTURE.datum_support flag.

which seems like a sensible standard warning that can be dealt with by setting iris.FUTURE.datum_support = True.

However, if I do this and then run the regrid, I get an error:

pyproj.exceptions.CRSError: Invalid projection: +proj=lonlat +datum=Ordnance Survey of Great Britain 1936 +no_defs +type=crs: (Internal Proj Error: proj_create: Error 1027 (Invalid value for an argument): Unknown value for datum)

This actually maes sense because the datum Ordnance Survey of Great Britain 1936 doesn't appear to be the right way to specify the correct datum which should be EPSG:6277, but I can't find out where this datum is stored in either the source or target cubes, or how to edit that information and give it the right information.

Any help gratefully received...
Ubuntu 22.04 LTS, Python 3.11.4, iris 3.7.0, iris_grib 0.19.0

[nickschurch]

I think it might be somewhere in the crs of the target cube, since this has a datum attribute and produces the same error:

>>> target_cube.coord_system().as_cartopy_crs().datum
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/nick/mambaforge/envs/resas-22-27-d2-1/lib/python3.11/site-packages/iris/coord_systems.py", line 710, in as_cartopy_crs
    return ccrs.TransverseMercator(
           ^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/nick/mambaforge/envs/resas-22-27-d2-1/lib/python3.11/site-packages/cartopy/crs.py", line 1448, in __init__
    super().__init__(proj4_params, globe=globe)
  File "/home/nick/mambaforge/envs/resas-22-27-d2-1/lib/python3.11/site-packages/cartopy/crs.py", line 667, in __init__
    super().__init__(*args, **kwargs)
  File "/home/nick/mambaforge/envs/resas-22-27-d2-1/lib/python3.11/site-packages/cartopy/crs.py", line 198, in __init__
    super().__init__(self.proj4_init)
  File "/home/nick/mambaforge/envs/resas-22-27-d2-1/lib/python3.11/site-packages/pyproj/crs/crs.py", line 348, in __init__
    self._local.crs = _CRS(self.srs)
                      ^^^^^^^^^^^^^^
  File "pyproj/_crs.pyx", line 2378, in pyproj._crs._CRS.__init__
pyproj.exceptions.CRSError: Invalid projection: +proj=tmerc +datum=Ordnance Survey of Great Britain 1936 +lon_0=-2.0 +lat_0=49.0 +k=0.9996012717 +x_0=400000.0 +y_0=-100000.0 +units=m +no_defs +type=crs: (Internal Proj Error: proj_create: Error 1027 (Invalid value for an argument): Unknown value for datum)

[HGWright]

Thanks for opening this @nickschurch! It would really help us to see the file (or a subset of it).

[nickschurch]

Grid_OSGB_2km.zip
Here you go. Happy datum-fun...

[trexfeathers]

Definitely echo that a file example would be great, as we are still learning about datums.

In terms of fixing: while you can access target_cube.coord_system(), the coordinate system is actually a member of a coordinate originally, so that is where edits should take place. Coordinate do not share coordinate system objects, so some effort is needed to keep them aligned.

horizontal_coords = [target_cube.coord(axis=a, dim_coords=True) for a in ("x", "y")]
for coord in horizontal_coords:
    if hasattr(coord, "coord_system"):
        coord.coord_system.datum = "OSGB36"

OSGB36 actually refers to EPSG:27700, rather than 6277. I think they are equivalent, and I don't know the string that might correspond to 6277. (This is probably part of the flaw with Proj strings - you can only define a predetermined selection of things, as opposed to WKT which can define things from scratch).

[nickschurch]

Yep, 6277 is a part of 27700 so I think you are right, that should be the code to set. It's a bit weird to me that the iris structure keeps coordinated seperat for the coordinate system - surely a recipe for a mis-match mess.

I'll try this and see if it works though. Do I need to edit the vertical coordinates similarly and how do I write them back to the cube sensibly?

The below doesn't work, so I guess there's some setting function I need...

target_cube.coord(axis="x")  = horizontal_coords[1]
target_cube.coord(axis="y")  = horizontal_coords[2]

[trexfeathers]

Hi @nickschurch, sorry this got lost during the festive period.

I'll try this and see if it works though. Do I need to edit the vertical coordinates similarly and how do I write them back to the cube sensibly?

No, the vertical coordinates are not considered when Iris is working with projections (by the way the API for updating ANY coordinate is the same).

The below doesn't work, so I guess there's some setting function I need...

target_cube.coord(axis="x")  = horizontal_coords[1]
target_cube.coord(axis="y")  = horizontal_coords[2]

coord_to_remove = target_cube.coord(axis="x")
cube_dims = coord_to_remove.cube_dims(target_cube)
target_cube.remove_coord(coord_to_remove)
# Or add_aux_coord(), as appropriate.
target_cube.add_dim_coord(horizontal_coords[1], cube_dims)