NEW FEATURE: Visualization of UGRID elements defined on the nodes about psy-maps HOT 9 OPEN

@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

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?

from psy-maps.

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.

from psy-maps.

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.

from psy-maps.

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

from psy-maps.

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/

from psy-maps.

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).

from psy-maps.

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).

from psy-maps.

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

from psy-maps.

oh, this issue has been accidently closed!

from psy-maps.