Investigate: subsetting a MeshCoord should return a MeshCoord?

[trexfeathers]

I might have some code that could help with this:

iris/lib/iris/tests/stock/__init__.py

Lines 840 to 886 in 6adac1b

	def remove_duplicate_nodes(mesh: ugrid.Mesh):
	"""Remove duplicate nodes from a mesh.
	Mesh.from_coords() does not do this due to complications like lazy
	arrays. Not a problem here.
	"""
	# TODO:
	# Contained in a function because this _could_ be generalised into a
	# public function. Would need to make it handle Dask arrays and masked
	# indices.
	# Example nodes: [a, b, c, a, c, b, d]
	# (Real nodes are X-Y pairs so a 2d array).
	# Example faces: [[0, 1, 2, 6], [3, 4, 5, 6]]
	# I.e. faces made by connecting: a-b-c-d , a-c-b-d
	# Processed nodes: [a, b, c, d]
	# Processed faces: [[0, 1, 2, 3], [0, 2, 1, 3]]
	nodes = np.stack([c.points for c in mesh.node_coords])
	face_node = mesh.face_node_connectivity
	# first_instances = a full length array but always with the index of
	# the first instance of each node e.g.: [0, 1, 2, 0, 2, 1, 3]
	nodes_unique, first_instances = np.unique(
	nodes,
	return_inverse=True,
	axis=1,
	)
	# E.g. indexing [0, 1, 2, 0, 2, 1, 3] with [[0, 1, 2, 6], [3, 4, 5, 6]]
	# -> [[0, 1, 2, 3], [0, 2, 1, 3]]
	indices_unique = first_instances[face_node.indices]
	# Connectivity indices expected to be a masked array.
	indices_unique = np.ma.masked_array(indices_unique, mask=face_node.indices.mask)
	# Replace the original node coords and face-node connectivity with the
	# unique-d versions.
	node_x, node_y = [
	AuxCoord(nodes_unique[i], **c.metadata._asdict())
	for i, c in enumerate(mesh.node_coords)
	]
	mesh.add_coords(node_x=node_x, node_y=node_y)
	conn_kwargs = dict(indices=indices_unique, start_index=0)
	mesh.add_connectivities(
	ugrid.Connectivity(**(face_node.metadata._asdict() | conn_kwargs))
	)
	return mesh

[trexfeathers]

This is very closely linked to #4751 - it's probably a choice between 1 and the other

[stephenworsley]

Another potential route for subsetting a mesh is proposed here #6014.

---

After considering this I see essentially two alternatives to the current behaviour for subsetting a MeshCoord.

1. A MeshCoord is returned whose Mesh is a reduced version of the original.
2. A MeshCoord (or similar Coord) is returned whose mesh is a "view" onto the original Mesh.

This second alternative, implemented above as the MeshIndexSet class, has some precedence in the UGRID concet of a "location index set". It would be more sophisticated and difficult to implement, but could be useful for types of operations which require splitting up and recombining data from a mesh.

It should be noted that a MeshCoord superficially appears like an AuxCoord (indeed, this appearance is why the existing behaviour "falls out" as a natural consequence). Therefore, changing the behaviour to return MeshCoords yields a result which looks superficially similar to the existing behaviour, and only adds extra information. This means there is a good case to be made that changing this behaviour would not be a breaking change, especially since returning a MeshCoord as the result of a slice of a MeshCoord already seems like the more natural result.
There is a similar case to make that the MeshIndexSet only adds information, however since this does not seem like the more natural result, I would be hesitant about the idea of returning this as the result of a slice of a MeshCoord and would consider this more likely to be a breaking change. There is still a good case that this could be the result of some other method of subsetting a MeshCoord.

In the meantime, it is worth highlighting that the existing behaviour of returning an AuxCoord may be liable to change if we are unable to put in a change before taking mesh code out of experimental. A warning when slicing a MeshCoord may suffice.

[trexfeathers]

The below code demonstrates the currently typical way of recreating a mesh after subsetting (#4751 refers to this as 'remeshing').

If we make the suggestions here work then the below code could create some confusing situations - no MeshCoord will have been replaced with an AuxCoord. But it will be easy to avoid that happening if we just add a clause into Mesh.from_coords() that knows to perform a 'no-op' if the coordinates passed in are already MeshCoords - just return the existing Mesh.

If we made that work I would be comfortable implementing changes even after mesh support becomes non-experimental.

iris/lib/iris/experimental/geovista.py

Lines 324 to 333 in 409f92c

	coords_on_mesh_dim = region_cube.coords(dimensions=mesh_dim)
	new_mesh = Mesh.from_coords(
	*[c for c in coords_on_mesh_dim if c.has_bounds()]
	)
	new_mesh_coords = new_mesh.to_MeshCoords(cube.location)
	for coord in new_mesh_coords:
	region_cube.remove_coord(coord.name())
	region_cube.add_aux_coord(coord, mesh_dim)