NEW FEATURE: Visualization of UGRID elements defined on the nodes

[Chilipp]

Summary

Within the UGRID conventions it is possible to defined variables to live on the nodes, rather than the faces. Within psyplot, we do visualize them at the the moment, but silently this provides wrong results.

Reason

It should be possible to visualize these elements, too.

Detailed explanation

Within the get_cell_node_coord method of the UGRID decoder, we generate triangles using a delauney triangulation.

https://github.com/psyplot/psyplot/blob/36a23ce964d58e41468fd1f97188b75454c14b7d/psyplot/data.py#L1807

I think this is fine for now as we do not have generic methods to generate grids from the nodes. What grid is the best is a very scientific question and should rather be answered by a custom decoder class for the data.

None the less, the current implementation in psy-maps (and psy-simple) is wrong. In the _polycolor method, we say transformed, array=arr.ravel(). For variables on a node however, this does not generate correct results because the length of the array (which is the same as the number of nodes) is not the same as for transformed, which has the length of the number of triangles.

So array should actually be the mean of the nodes for each generated face element (i.e. for each generated triangle).

Examples

matplotlibs tripcolor method is doing exactly this: https://github.com/matplotlib/matplotlib/blob/e097bf4baf8f275fda91f224b537076caf17dd91/lib/matplotlib/tri/tripcolor.py#L108

ping @platipodium

[Chilipp]

@platipodium, I'd very much like to get your feedback here. Consider the test file for psyplot: https://github.com/psyplot/psyplot/blob/master/tests/simple_triangular_grid_si0.nc

ncdump -h simple_triangular_grid_si0.nc

netcdf simple_triangular_grid_si0 {
dimensions:
	nMesh2_node = 4 ;
	nMesh2_face = 2 ;
	Two = 2 ;
	Three = 3 ;
	time = UNLIMITED ; // (1 currently)
variables:
	int Mesh2 ;
		Mesh2:cf_role = "mesh_topology" ;
		Mesh2:long_name = "Topology data of 2D unstructured mesh" ;
		Mesh2:topology_dimension = 2 ;
		Mesh2:node_coordinates = "Mesh2_node_x Mesh2_node_y" ;
		Mesh2:face_node_connectivity = "Mesh2_face_nodes" ;
	float Mesh2_node_x(nMesh2_node) ;
		Mesh2_node_x:standard_name = "longitude" ;
		Mesh2_node_x:long_name = "Longitude of 2D mesh nodes" ;
		Mesh2_node_x:units = "degrees_east" ;
	float Mesh2_node_y(nMesh2_node) ;
		Mesh2_node_y:standard_name = "latitude" ;
		Mesh2_node_y:long_name = "Latitude of 2D mesh nodes" ;
		Mesh2_node_y:units = "degrees_north" ;
	int Mesh2_face_nodes(nMesh2_face, Three) ;
		Mesh2_face_nodes:cf_role = "face_node_connectivity" ;
		Mesh2_face_nodes:long_name = "Maps every triangular face to its three corner nodes" ;
		Mesh2_face_nodes:start_index = 0 ;
	float time(time) ;
		time:standard_name = "time" ;
		time:units = "days since 1951-01-01 00:00:00" ;
	float Mesh2_ndvar(time, nMesh2_node) ;
		Mesh2_ndvar:long_name = "node variable" ;
		Mesh2_ndvar:standard_name = "node_variable" ;
		Mesh2_ndvar:units = "None" ;
		Mesh2_ndvar:mesh = "Mesh2" ;
		Mesh2_ndvar:location = "node" ;
	float Mesh2_fcvar(time, nMesh2_face) ;
		Mesh2_fcvar:long_name = "face variable" ;
		Mesh2_fcvar:standard_name = "face_variable" ;
		Mesh2_fcvar:units = "None" ;
		Mesh2_fcvar:mesh = "Mesh2" ;
		Mesh2_fcvar:location = "face" ;

// global attributes:
		:title = "test mesh" ;
		:institution = "Universitaet Hamburg" ;
		:contact = "None" ;
		:source = "None" ;
		:references = "None" ;
		:comment = "None" ;
		:Conventions = "UGRID-0.9" ;
		:creation_date = "2015-01-26 09:19:01  01:00" ;
		:modification_date = "2015-01-26 09:19:01  01:00" ;
}

the Mesh2_ndvar variable lives on the nodes. A call to

ds = psy.open_dataset("simple_triangular_grid_si0.nc")
plt.tripcolor(ds.Mesh2_node_x.values, ds.Mesh2_node_y.values, ds.Mesh2_ndvar.values[0])
plt.scatter(ds.Mesh2_node_x.values, ds.Mesh2_node_y.values, color="red")

gives something like this

i.e. the nodes (red dots) are considered as the nodes of the triangles. Is this an approach you'd consider as a valid visualization of the node elements?

[Chilipp]

actually @platipodium I understand now what you meant by the dual mesh. Do you know an efficient library for python to calculate this from the UGRID conventions? I think this would be the best in terms of visualization.

[platipodium]

The above visualization represents as face color the mean of the surrounding nodes. So it is valid, in a away. I don't think this is what people would like to see, though. There should be a difference between the nodes on the end of the diagonal (they have different values)

• This visualization should be possible
• A different visualization with faces on the dual mesh is probably better

I am not aware of plotting libraries that do this, unfortunately.

[Chilipp]

Alright, thanks for the feedback @platipodium ! I also could not find something that does this. But having thought about it, it's not too difficult to come up with an algorithm to calculate the dual mesh, so I'll probably do this and keep you posted

[platipodium]

An industry-standard dual-mesh algorithm should be found in the Earth System Modeling Framework https://github.com/esmf-org/esmf. Probably implemented in C++, but they do have a python Interface, too, which might show how to access the C++ backend.

And then, there is the qhull library with its Python wrapper https://pypi.org/project/pyhull/

[Chilipp]

An industry-standard dual-mesh algorithm should be found in the Earth System Modeling Framework https://github.com/esmf-org/esmf

yes, I saw these two, too. The ESMF library would probably be the best, but the documentation is extremely sparse and there are lots and lots of broken links inside. I could not manage to use it.

with the simple_triangular_grid_si0.nc for instance, it should be pretty straight-forward to load it via

import ESMF
mesh = ESMF.Mesh(filename="simple_triangular_grid_si0.nc", filetype=ESMF.FileFormat.UGRID)

but I just get a not very helpful error log with

20201121 132418.278 ERROR            PET0 ESMF_Mesh_C.F90:176 f_esmf_meshcreatefromfile Wrong argument specified  - - incorrect args for UGRID
20201121 132418.353 ERROR            PET0 ESMF_Mesh_C.F90:178 f_esmf_meshcreatefromfile Wrong argument specified  - Internal subroutine call returned Error

do you know anyone who is a bit more experienced with this library by chance?

---

And then, there is the qhull library with its Python wrapper https://pypi.org/project/pyhull/

I never used qhull so far, but it looks to me like a more basic software that you can use in case you don't have a mesh yet. This does not apply for our case, because we do have the primal mesh, just not the dual. But please correct me if I am wrong.

---

Is it really that difficult to come up with something? To me, I'd use an algorithm like this for nodes:

1. take the triangles from the face_node_connectivity
2. calculate the midpoints for each triangle
3. for each node:
   i. lookup the triangles that touch this node
   ii. assign their midpoints (i.e. center of mass) as a node of the dual mesh for this node
4. there will now be several nodes that do not have as many points as the others. These are the nodes at the boundary of the simulation domain. For these nodes,
   i. look up the edges that touch the node,
   ii. if the edge is contained in only one triangle, add the midpoint of the edge to the list of dual mesh nodes in 3ii.

and for the edges, it's pretty much the same. Only step 4 can be dropped:

1. take the triangles from face_node_connectivity
2. calculate the midpoints for each triangle
3. for each edge:
   i. lookup the triangles that have this edge (which will be 2 at most)
   ii. assign their midpoints (i.e. center of mass) as a node for the dual mesh

Am I oversimplyfing this? It's of course a bit of a challenge if you want to work with large meshes consisting of millions of triangles, but I have pretty good experience with using cython in this case, particularly for step 3.

Also I am not yet sure whether I can generalize this to use mixed meshes, i.e. meshes whose faces can be of any shape, i.e. triangular and quadratic. But I could live for the moment with triangular grids, only. I don't know of any use case where these mixed meshes occur.

---

In the end I'd say, the best strategy would be to use ESMF, but only if we have someone who has some experience with it and UGRID and can tell me how to generate the dual mesh from an existing UGRID-conform primal mesh.

If we cannot come up with someone, I think it would take me about two full days to implement the suggested algorithm, including tests. The latter should be feasible within the next three weeks (my schedule is pretty full and I am on vacation for one week).

[platipodium]

1. I have very good connections to ESMF people. I think their python interface is mostly managed by Ryan O'Kuingtthons, who had been doing some subcontracting work for us. Please write an email to [email protected].

2. I agree with your simple calculation. This should be easily extensible to quads (we actually operate on mixed meshes at HZG), simply take the center of mass of the quad. At the edge, you might have a triangle special situation.

3. Yeah, qhull is good for creating a new mesh. But they're more lightweight than ESMF and they might have (I don't know) also a call for getting the dual from an existing mesh. I would probably prefer using that library (but note their python interface hasn't been maintained for a while).

[Chilipp]

Alright! Thanks for the hints @platipodium ! I'll contact the esmf support tomorrow

[Chilipp]

oh, this issue has been accidently closed!