Preserve `crs_wkt` when loading and saving NetCDF files

[neilCrosswaite]

✨ Feature Request

Preserve crs_wkt when loading and saving NetCDF files

Motivation

To add to some of the open issues and PRs on Coordinate reference system well known text (crs_wkt) requirements.

• Adopt WKT for CRS #3796
• Save crs_wkt on netCDF save #4719

I would like to:

• load a cube from a netcdf file which defines a crs_wkt on a grid-mapping variable
• make changes to the cube that do not affect the coordinate system
• Save the cube back to a new file and verify that the crs_wkt remains intact and is identical to the input crs_wkt

An example of this would be load a cube, calculate the mean temperature with respect to time which on saving it should have the same coorinate metadata.

There is a current problem with this that the crs_wkt just doesn't get loaded nor re-saved.

But there is a future problem in #4719 implementation of it that it can change the crs_wkt

These can be seen with a file loaded and saved using the save from #4719.

The input file has the crs_wkt:

PROJCRS[
    "unknown",
    BASEGEOGCRS[
        "unknown",
        DATUM[
            "OSGB 1936",
            ELLIPSOID["Airy 1830", 6377563.396, 299.3249646, LENGTHUNIT["metre", 1]],
            ID["EPSG", 6277],
        ],
        PRIMEM[
            "Greenwich", 0, ANGLEUNIT["degree", 0.0174532925199433], ID["EPSG", 8901]
        ],
    ],
    CONVERSION[
        "unknown",
        METHOD["Mercator (variant B)", ID["EPSG", 9805]],
        PARAMETER[
            "Latitude of 1st standard parallel",
            0,
            ANGLEUNIT["degree", 0.0174532925199433],
            ID["EPSG", 8823],
        ],
        PARAMETER[
            "Longitude of natural origin",
            0,
            ANGLEUNIT["degree", 0.0174532925199433],
            ID["EPSG", 8802],
        ],
        PARAMETER["False easting", 0, LENGTHUNIT["metre", 1], ID["EPSG", 8806]],
        PARAMETER["False northing", 0, LENGTHUNIT["metre", 1], ID["EPSG", 8807]],
    ],
    CS[Cartesian, 2],
    AXIS["(E)", east, ORDER[1], LENGTHUNIT["metre", 1, ID["EPSG", 9001]]],
    AXIS["(N)", north, ORDER[2], LENGTHUNIT["metre", 1, ID["EPSG", 9001]]],
]

When this is saved it is output as:

GEOGCRS[
    "OSGB36",
    DATUM[
        "Ordnance Survey of Great Britain 1936",
        ELLIPSOID["Airy 1830", 6377563.396, 299.3249646, LENGTHUNIT["metre", 1]],
    ],
    PRIMEM["Greenwich", 0, ANGLEUNIT["degree", 0.0174532925199433]],
    CS[ellipsoidal, 2],
    AXIS[
        "geodetic latitude (Lat)",
        north,
        ORDER[1],
        ANGLEUNIT["degree", 0.0174532925199433],
    ],
    AXIS[
        "geodetic longitude (Lon)",
        east,
        ORDER[2],
        ANGLEUNIT["degree", 0.0174532925199433],
    ],
    USAGE[
        SCOPE["Geodesy."],
        AREA[
            "United Kingdom (UK) - offshore to boundary of UKCS within 49°45'N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore."
        ],
        BBOX[49.75, -9.01, 61.01, 2.01],
    ],
    ID["EPSG", 4277],
]

Test Examples

I have written example tests for netCDF intergration, which can be viewed on my iris fork: View example tests.

The test loads a file which has a crs_wkt, tests that the crs_wkt is as expected, then it saves the cube and uses netCDF4 to verify that the correct metadata is still on the cube.

[trexfeathers]

Hi @neilCrosswaite,

Firstly, we're sorry for the frustrating slow progress on this topic, it's just the fastest we can go given the wide scope of the development team. It's good to know there is interest. You can see our 'plan of attack' in the task list on #3796.

But there is a future problem in #4719 implementation of it that it can change the crs_wkt

Iris has no contract regarding file-Cube-file round-tripping. We feel it's an important part of Iris' role to enforce correct file formatting - you can read more here #5165 - so Iris cannot guarantee that files will always be output identically to how they are input. Obviously this does NOT mean that Iris changes things just for the sake of it, we do make efforts to preserve what we can.

#4719 is not a finished solution, and we will need to review your example in more detail to see if Iris is being overzealous versus what we intended.

Unfortunately I will be unavailable for further discussion for several weeks after today.

Note for developers working on #3796

• Check whether Iris' modification of crs_wkt in this example is expected. Correct Save crs_wkt on netCDF save #4719 if necessary

[mo-marqh]

i'd like to add some support for this issue as raised by @neilCrosswaite, from another user community's perspective.

we are evaluating adding crs-wkt metadata to model outputs. if such metadata is provided for files loaded into an environment, and data from those files is then output, perhaps after cube filtering, aggregation, etc. but where no changes to the spatial referencing have occurred, then i'd expect as a user that the provided crs_wkt definitions would be preserved in the output file.

I can understand newly created coordinate reference system instances not having a crs_wkt, but existing ones loaded from files should assume that the data provider knows what they are doing and why, I feel. I'd expect this to be preserved (as an attribute?) and available for export to netCDF files.
i'd also like to be able to set this attribute on a coordinate reference system instance within some user controlled Python code, but that's another related ticket, perhaps.

[trexfeathers]

From @SciTools/peloton: a note for developers. Once we have an implementation, it will be worth investigating whether ncdata can be used to augment cases where Iris is too opinionated.

[trexfeathers]

We have now finished our planned work on the extended grid mapping syntax and WKT (#6357), to be released imminently. This very much does not deliver what you are asking for, sorry.

The situation is closely analogous to how we handle the axis attribute. axis and crs_wkt are ascribed meaning by CF, but both are considered 'supplementary' to other information, so the information they convey is already available elsewhere in the Iris data model ¹. Because they do not add value, and in some cases risk divergence - e.g. where grid_mapping says one thing but crs_wkt says another - on balance we choose not to track this information within Iris. That means we don't parse it on load, and we regenerate it on save (since all the necessary information is available elsewhere).

This obviously is no good for your use case; ways forward:

Track WKT within the Iris data model

A reasonable amount of added complexity for something that makes no difference to Iris behaviour, only to file-file round-tripping, but can still be discussed if important to enough users

• A space for tracking crs_wkt in the Iris data model. Currently the CoordSystem class does not carry attributes at all, although maybe this could be a normal class member if it's a one-off
• An ability to prevent divergence between the CoordSystem instance and the WKT
  • Minimum: any changes to the CoordSystem mean the WKT gets deleted
  • Or: an actual understanding of the WKT meaning so that alignment can be properly validated

Workaround

We know Iris cannot be all things to all people, and that is one of the reasons we created Ncdata - to help people customise Iris NetCDF I/O to suit their specific needs. The below code demonstrates:

from ncdata.iris import from_iris, to_iris
from ncdata.netcdf4 import from_nc4, to_nc4
from ncdata.threadlock_sharing import enable_lockshare

enable_lockshare(iris=True)

pre_load = from_nc4("wkt.nc")
grid_mapping = pre_load.variables["latitude_longitude"]
wkt_original = grid_mapping.get_attrval("crs_wkt")

cubes = to_iris(pre_load)
#
# Insert Iris operations here
#
post_save = from_iris(cubes)

grid_mapping = post_save.variables["latitude_longitude"]
grid_mapping.set_attrval("crs_wkt", wkt_original)
to_nc4(post_save, "saved.nc")

Footnotes

1. I am aware there is some advanced information ONLY available in crs_wkt, but these are not things that Iris handles or has plans to handle ↩

[trexfeathers]

(Removing this from #3796 as I imagine the conversation might run for a while longer)

[mo-marqh]

The situation is closely analogous to how we handle the axis attribute. axis and crs_wkt are ascribed meaning by CF, but both are considered 'supplementary' to other information, so the information they convey is already available elsewhere in the Iris data model 1. Because they do not add value, and in some cases risk divergence - e.g. where grid_mapping says one thing but crs_wkt says another - on balance we choose not to track this information within Iris. That means we don't parse it on load, and we regenerate it on save (since all the necessary information is available elsewhere).

1. I am aware there is some advanced information ONLY available in crs_wkt, but these are not things that Iris handles or has plans to handle

I think that this analysis does not fully recognise the fact that well known text coordinate reference systems are far more rich and nuanced than the simple coordinate reference system analogues that CF refers to as grid mappings, despite the awareness note.
This was a primary justification for the inclusion of crs_wkt within the CF data model.
There are numerous cases where the information to properly decide on how to choose elements of a crs_wkt string cannot be encoded within a CF grid mapping element set. I presume that this applies to the Iris crs representation model as well.

But, i think i'm just restating what you already note here.

thank you for sharing the ways forward

I agree with the assessment that a full implementation of the ISO19111 data model and well known text encoding is far more than is really useful within the scoping of Iris.

I can see some value in storing any provided crs_wkt string from a netCDF file on load, an operating a simple trust mechanism, where either:

1. nothing has changed, keep the crs_wkt string and output it on write.
2. any changes to the CoordSystem mean the WKT gets deleted

this would seem like a useful feature.

I'll try to experiment with the ncdata package and validate it's usage as a work around option; I think i can see the intent from your code snippet.

[trexfeathers]

From @SciTools/peloton: thanks for the detail, this does sound like more work would be valuable, and you've provided a sensible way to either preserve or discard crs_wkt. I will make sure this doesn't go off the radar