Open Chilipp opened 3 years ago
@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
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?
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.
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)
I am not aware of plotting libraries that do this, unfortunately.
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
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/
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:
and for the edges, it's pretty much the same. Only step 4 can be dropped:
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).
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 esmf_support@ucar.edu.
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).
Alright! Thanks for the hints @platipodium ! I'll contact the esmf support tomorrow
oh, this issue has been accidently closed!
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 saytransformed, array=arr.ravel()
. For variables on a node however, this does not generate correct results because the length of thearray
(which is the same as the number of nodes) is not the same as fortransformed
, 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#L108ping @platipodium