pyugrid: locate_faces doesn't account for wrapping

[djkirkham]

A face which straddles the meridian may be incorrectly returned by UGrid.locate_faces(), or not returned when it should be. E.g.:

>>> from gridded.pyugrid.ugrid import UGrid
>>> grid = UGrid(node_lon=[10, 10, 350, 350],                                                                               
...              node_lat=[40, 50, 50, 40],                                                    
...              faces=[[0, 1, 2, 3]],
...              edges=[[0, 1], [1, 2], [2, 3], [3, 0]])
>>> grid.locate_faces([5, 45]) # Should return 0
-1
>>> grid.locate_faces([15, 45]) # Should return -1
0

It should be possible to infer where the 'inside' of a face is by using the fact that the nodes of its vertices should be specified in counter-clockwise order.

[jay-hennen]

There's a lot of problematic things happening here.

Firstly, UGrid doesn't support 4-sided polygons. The fact that locate_faces is doing the correct thing in this case is more incidental than intentional.

Secondly, what is happening under the hood is that the polygon being tested is 'wrapping the other way around the world'. The library that powers the locate_faces function is not globe-aware or even grid-aware in any capacity; it is a pure math library and is treating the nodes provided as if they were points on an infinite 2D mathematical plane, and the faces describing polygons formed from the points. It's also a hard requirement that when laid out in this manner the polygons formed don't overlap.

[ChrisBarker-NOAA]

UGrid doesn't support 4-sided polygons.

Actually, I think it does -- at least partially, though I don't think we ever the updated cell_tree2D to work with them.

But yes, none of the rest of algorithms know that the earth is round. This is a classic problem / limitation of virtually all software that deal with this GIS analysis, mapping, etc.

And knowing the winding order doesn't actually help -- whether you go from point: (10, 50) to (350, 50) or from (350, 50) to (10, 350) -- you can either go the short or long way around the earth -- whichever way the quad is wound, it could be either a small quad around the prime meridian or a really wide quad wrapped mostly around the earth.

(which brings up a heuristic that could (and often is) be used -- always go the short way. This would certainly work for small things like grid cells, but would require some potentially expensive checks, and may not be compatible at all with the CellTree data structure.

For the most part, this is dealt with by choosing a coordinate system that doesn't wrap -180 to 360 or whichever. As long as you are consistent, code the is pretending it's all cartesian will work fine.

None of that helps if you have a grid that truly wraps around the earth -- global models (which may have issues at the poles, too. But in this case, you have a periodic grid which will have to be dealt with specially in other says as well.

I would think the "cheating" way would be to duplicate the cells that cross the dateline, and deal that duplication outside Gridded somehow.

Patches welcome :-)

[djkirkham]

interesting to learn that 4-sided faces are not fully supported. We (my team at the Met Office) are planning a piece of work to support easier analysis and plotting of cubedsphere grids, and as part of that we will be looking at improvements we can submit to the pyugrid code. As such it would be useful to get an idea of what it would take to support 4-sided faces.

knowing the winding order doesn't actually help.

I hadn't realised that knowing the winding order doesn't fully specify where a face - or even an edge - specified by nodes is located. Does that mean that the UGrid conventions themselves are ambiguous in this regard?

[ChrisBarker-NOAA]

As for quad grids -- someone added support in pyugrid a a while back, I'm not sure f it made it inot gridded at this point, but it would not take much to do so.

A bit more work is required to support unstructured quad grids in cellTree2D -- it's in C++ and has triangles hard-coded in. That's need iff you need fast cell location. And celltree2D does support structured quad grids, so it shouldn't be too hard.

As for ambiguity in the topology, it's not ambiguous in 2D -- ti's only ambiguous on a sphere. geometry on a sphere is so much more complicated (and expensive) that geometry libraries generally don't deal with it. And any non-global model can simply use a coordinate system that does wrap -- it gets us pretty far.

But, of course, for global models, that won't cut it. If you're working with global models, then you've figured out ways to deal with it -- so we'd be glad to have help updateing the code to handle it.

As for the UGrid convention -- it is currently silent on how to deal with periodic grids. I recall a conversation about it a few months ago -- but can't find it right now.

Maybe you could start that on the ugrid conventions repo:

https://github.com/ugrid-conventions/ugrid-conventions

[jay-hennen]

A bit more work is required to support unstructured quad grids in cellTree2D -- it's in C++ and has triangles hard-coded in. That's need iff you need fast cell location. And celltree2D does support structured quad grids, so it shouldn't be too hard.

Actually, I think cell_tree2D is handling the situation fine in this case, and dealing with the quad correctly. What I'm less confident of is pyugrid's support for quad unstructured grids. I remember specifically adding code for and testing cell_tree with unstructured quads and more, but I've never done so with pyugrid.

[ChrisBarker-NOAA]

Thanks @jay.

So it looks like quad grids are supported but not tested.

So that “just” leaves you with the world is round problem :-)