Comments (9)
@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?
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:
- take the triangles from the face_node_connectivity
- calculate the midpoints for each triangle
- 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 - 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:
- take the triangles from face_node_connectivity
- calculate the midpoints for each triangle
- 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.
-
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].
-
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.
-
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.
Related Issues (18)
- NEW FEATURE: Implement formatoption widgets
- Problem in displaying dateline crossing for rotated coordinates HOT 2
- Variables without time dimensions greyed out HOT 10
- Curvilinear Datagrid is not displayed correctly HOT 5
- Projection update does not work with cartopy 0.18
- ICON grid cells at poles are not displayed correctly
- Plotting operational ICON data and benchmarking HOT 6
- Deploy psy-maps docs with CircleCI and GH-Pages
- Add scipy as dependency
- Rotated coordinates are not shown correctly with plot="poly"
- Move psy-maps to psyplot organization
- display bug with icon-data HOT 2
- Mapplot covers whole figure when using map_extent and saved as .png HOT 1
- mapvector: lonlatbox and density don't work together HOT 5
- Setting a `matplotlib.colors` normalization as `bounds` causes error in `mapcombined`
- ortho projection does not focus on central longitude and latitude
- Strange Figures for rotated coordinates HOT 3
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from psy-maps.