search

SciTools / iris-esmf-regrid

A collection of structured and unstructured ESMF regridding schemes for Iris.
https://iris-esmf-regrid.readthedocs.io/en/latest
BSD 3-Clause "New" or "Revised" License
19 stars 17 forks source link

Cannot load a minimalistic ugrid file #188

Closed pletzer closed 2 years ago

pletzer commented 2 years ago

🐛 Bug Report

Cannot load a mesh from a ugrid file that contains:

variables: int Mesh2d ; Mesh2d:cf_role = "mesh_topology" ; Mesh2d:edge_node_connectivity = "Mesh2d_edge_nodes" ; Mesh2d:face_edge_connectivity = "Mesh2d_face_edges" ; Mesh2d:face_node_connectivity = "Mesh2d_face_nodes" ; Mesh2d:long_name = "Topology data of 2D unstructured mesh" ; Mesh2d:node_coordinates = "Mesh2d_node_x Mesh2d_node_y" ; Mesh2d:topology_dimension = 2 ; int Mesh2d_face_edges(nMesh2d_face, nMesh2d_vertex) ; ...

How To Reproduce

Steps to reproduce the behaviour:

  1. Download https://github.com/pletzer/mint/blob/master/data/lfric_diag_wind.nc
  2. Run the script below

import iris from iris.experimental.ugrid import PARSE_UGRID_ON_LOAD with PARSE_UGRID_ON_LOAD.context(): mesh = iris.load_cube("lfric_diag_wind.nc")

You'll get

(iris-dev) niwa-35791~/iris-test$ python terror.py Traceback (most recent call last): File "terror.py", line 4, in mesh = iris.load_cube("lfric_diag_wind.nc") File "/Users/pletzera/iris/lib/iris/init.py", line 341, in load_cube cubes = _load_collection(uris, constraints, callback).cubes() File "/Users/pletzera/iris/lib/iris/init.py", line 276, in _load_collection result = _CubeFilterCollection.from_cubes(cubes, constraints) File "/Users/pletzera/iris/lib/iris/cube.py", line 104, in from_cubes for cube in cubes: File "/Users/pletzera/iris/lib/iris/init.py", line 261, in _generate_cubes for cube in iris.io.load_files(part_names, callback, constraints): File "/Users/pletzera/iris/lib/iris/io/init.py", line 208, in load_files for cube in handling_format_spec.handler( File "/Users/pletzera/iris/lib/iris/fileformats/netcdf.py", line 974, in load_cubes mesh_coords, mesh_dim = _build_mesh_coords(mesh, cf_var) File "/Users/pletzera/iris/lib/iris/experimental/ugrid/load.py", line 492, in _build_mesh_coords mesh_coords = mesh.to_MeshCoords(location=cf_var.location) File "/Users/pletzera/iris/lib/iris/experimental/ugrid/mesh.py", line 1972, in to_MeshCoords result = [ File "/Users/pletzera/iris/lib/iris/experimental/ugrid/mesh.py", line 1973, in self.to_MeshCoord(location=location, axis=ax) for ax in self.AXES File "/Users/pletzera/iris/lib/iris/experimental/ugrid/mesh.py", line 1947, in to_MeshCoord return MeshCoord(mesh=self, location=location, axis=axis) File "/Users/pletzera/iris/lib/iris/experimental/ugrid/mesh.py", line 2835, in init points, bounds = self._construct_access_arrays() File "/Users/pletzera/iris/lib/iris/experimental/ugrid/mesh.py", line 3046, in _construct_access_arrays points_coord = self.mesh.coord(include_edges=True, axis=axis) File "/Users/pletzera/iris/lib/iris/experimental/ugrid/mesh.py", line 1614, in coord result = self._coord_manager.filter( File "/Users/pletzera/iris/lib/iris/experimental/ugrid/mesh.py", line 2283, in filter raise CoordinateNotFoundError(emsg) iris.exceptions.CoordinateNotFoundError: 'Expected to find exactly 1 coordinate, but found none.'

Expected behaviour

One should be able to load a ugrid file with variables that don't have coordinates. This would be expected for "extensive" variables, which are defined over the entire cell, face or edge. (All nodal variables are expected, however, to have coordinates.)

Screenshots

N/A

Environment

Additional context

Click to expand this section... ``` Add additional verbose information in a collapsible section. ``` See [here](https://gist.github.com/pierrejoubert73/902cc94d79424356a8d20be2b382e1ab) for further details.
github-actions[bot] commented 2 years ago

✨ Congratulations! ✨ Thanks for submitting your first issue to iris-esmf-regrid. We really appreciate it and will get back to you as soon as possible. Awesome job 👍

stephenworsley commented 2 years ago

This looks like it's due to a known problem in Iris where you aren't able to construct a "pointless" coordinate (see item 3 in this discussion https://github.com/SciTools/iris/discussions/4438). Due to the way we're interpreting meshes through the existing coordinate infrastructure in iris, we can only handle the case where face coords and edge coords exist (when data is stored on faces or edges). This is due to existing assumptions that all corrdinates ought to have points. We should be able to get this sorted eventually as it is a priority for future Iris/UGRID development, but for the time being you may want to use a work around.

I think it should still be possible to load the mesh in this file with iris.experimental.ugrid.load_mesh, then it should be possible to add "dummy coordinates" to the mesh, then you can generate MeshCoord's to describe the mesh using the to_MeshCoord method, then you can add those coords to a cube loaded from the data without the PARSE_UGRID_ON_LOAD context. That would look something like:

import numpy as np

import iris
from iris.experimental.ugrid import PARSE_UGRID_ON_LOAD, load_mesh

uav_cube = iris.load_cube("lfric_diag_wind.nc", "upward_air_velocity")
ew_cube = iris.load_cube("lfric_diag_wind.nc", "eastward_wind")
with PARSE_UGRID_ON_LOAD.context():
  mesh = load_mesh("lfric_diag_wind.nc")

edge_length = mesh.edge_node_connectivity.indices_by_location().shape[0]
face_length = mesh.face_node_connectivity.indices_by_location().shape[0]

dummy_edge_lat = AuxCoord(np.zeros(edge_length), standard_name="latitude")
dummy_edge_lon = AuxCoord(np.zeros(edge_length), standard_name="longitude")
dummy_face_lat = AuxCoord(np.zeros(face_length), standard_name="latitude")
dummy_face_lon = AuxCoord(np.zeros(face_length), standard_name="longitude")

mesh.add_coords(
  edge_x=dummy_edge_lon,
  edge_y=dummy_edge_lat,
  face_x=dummy_face_lon,
  face_y=dummy_face_lat,
)

for coord in mesh.to_MeshCoords("edge"):
   ew_cube.add_aux_coord(coord, 2)

for coord in mesh.to_MeshCoords("face"):
   uav_cube.add_aux_coord(coord, 2)
pp-mo commented 2 years ago

This looks like it's due to a known problem in Iris where you aren't able to construct a "pointless" coordinate (see item 3 in this discussion SciTools/iris#4438).

I've been considering fixing that problem, so if now is a good time I'd like to explore that. The reason we deferred it is that allowing a Coord to have bounds but no points would be rather a big change to the Iris datamodel. However, I feel that we could implement a somewhat "stopgap" solution, whereby a suitable MeshCoord would be created, but will raise an error only if an attempt to access the .points is made.

Would it help to have a prototype version of Iris which implements that (e.g. a feature branch) ?

github-actions[bot] commented 2 years ago

@SciTools-incubator/esmf-regrid-devs This issue is stale due to a lack of activity in the last 90 days. Remove stale label or comment, otherwise this issue will close automatically in 7 days time.

github-actions[bot] commented 2 years ago

@SciTools-incubator/esmf-regrid-devs This stale issue has been automatically closed due to no community activity