ISEA3H/IVEA3H: `getZoneWGS84Vertices` winds zones clockwise, opposite of the other grids

[ajfriend]

When using ISEA3H or IVEA3H, getZoneWGS84Vertices returns a zone's corners in clockwise order (viewed from outside the sphere). The other grids I checked (ISEA7H, IVEA7H, rHEALPix) return counterclockwise. Alternatively, getZoneRefinedWGS84Vertices returns counterclockwise for all five grids.

import dggal as al
al.pydggal_setup(al.Application())


def wise(pts):
    """Compute planar signed area on (lon, lat) to determine orientation.
    positive = counterclockwise.

    Returns CW or CCW accordingly.
    """
    n = len(pts)
    A = sum(
          pts[(i + 1) % n][1] * pts[i][0]
        - pts[(i + 1) % n][0] * pts[i][1]
        for i in range(n)
    )

    if A > 0:
        return 'CCW'
    else:
        return 'CW'


p = al.GeoPoint(42.1, 12.0)

for name in ('ISEA7H', 'IVEA7H', 'rHEALPix', 'ISEA3H', 'IVEA3H'):
    g = getattr(al, name)()
    z = g.getZoneFromWGS84Centroid(8, p)

    corners = [
        (float(v.lon), float(v.lat))
        for v in g.getZoneWGS84Vertices(z)
    ]
    refined = [
        (float(v.lon), float(v.lat))
        for v in g.getZoneRefinedWGS84Vertices(z, 8)
    ]

    print(f'{name:10s} corners: {wise(corners):3s}   refined: {wise(refined)}')

Output:

ISEA7H     corners: CCW   refined: CCW
IVEA7H     corners: CCW   refined: CCW
rHEALPix   corners: CCW   refined: CCW
ISEA3H     corners: CW    refined: CCW
IVEA3H     corners: CW    refined: CCW

Checking exhaustively at levels 0 to 4 (for both 3H grids), every zone comes back CW except two pentagons per level. The level 0 pentagons that come back CCW: AA-0-A and AB-0-A.

[jerstlouis]

@ajfriend Thanks for testing and detailing this.

There is currently no guarantee as to which order getZoneWGS84Vertices() returns vertices in.

Arguably, we should align this so it's CCW for everything, but since this is somewhat of an internal method and it might currently be associated to the neighbors orders etc., changing this now for the 3H grids would have some implications.

I would likely leave this issue open to tackle this eventually, but it is not a high priority and probably not something that would be fixed in 0.0.7 unless this is somehow important for some use case.

AA-0-A and AB-0-A are the polar pentagons (the icosahedron poles at the first vertex off the coast and Sweden and its opposite vertex near Antarctica).

[ajfriend]

Agreed that this isn't a high priority. For my own purposes, I can just reverse the rings, or use getZoneRefinedWGS84Vertices.

As with the other issues, I'm just raising them when I find them in case it ends up being useful to you or others. :)
This one was just surprising since I had some scripts that were working fine with the other grids, but broke on these due to the area computation coming out larger than expected due to the loop orientation.

As always, thanks for the quick response and attention to these!

[jerstlouis]

@ajfriend

area computation coming out larger than expected due to the loop orientation.

For precise (and very fast) area computation, getZoneArea() should be used. It is the true theoretical area of the zone (more accurate than calculating the geodetic area of the refined vertices which have a limited precision).

getZoneWGS84Vertices() is mostly informational. Almost all use cases should use getZoneRefinedWGS84Vertices() instead, since the actual edge between two vertices is not a geodesic, and so intermediate points are necessary to have a more accurate zone boundary (or area), and the raw vertices also will not form a proper polygon around the anti-meridian etc.

[ajfriend]

I have a use case where I'm primarily interested in a zone's "real" vertices, plus the "distortion" vertices where an edge crosses a face boundary. Any chance there could be a function that returns those (real + distortion) directly? The closest I can currently get is getZoneRefinedWGS84Vertices(), but that also samples the smooth non-geodesic curvature along each edge, when what I'm after is just the discrete change points.

[jerstlouis]

@ajfriend

I'm not sure what you mean by the "distortion" vertices. There are no such things in ISEA / IVEA projections!

In H3 specificlaly, the zone edges are one or two great circle arc segments on the internal authalic sphere, 1 when remaining inside one spherical icosahedron spherical triangle, and two when crossing that edge. That is because H3 uses a gnomonic projection for each of the 20 icosahedral faces, which maps great circles to straight lines in the projection. I think this is what you might be referring to? But even in H3, those great circle arcs on the internal authalic sphere actually do not map to a WGS84 geodesic, and individual refinement points along those great circle arcs should ideally be converted from authalic to geodetic latitude (and more generally, a different aspect of this is that coordinates should be converted between WGS84 geodetic and authalic latitude when using the H3 library, which does not do it automatically).

There is no such 1 or 2 great circle arc segments when crossing the face or not with the ISEA / IVEA / RTEA projection, because the Snyder projection does not map geodesic to the straight lines used to construct the grids. There are no discrete change points, it's a continuous function along the entire edge (and as explained above, due to the geodetic/authalic conversion, this is actually still true with the real H3 zone boundaries, unless you're working with a spherical Earth model -- the internal and topologically correct H3 zone edges are the great circle arcs on the internal authalic sphere, never geodesics on the WGS84 ellipsoid, whether you convert the authalic latitudes returned by the library to geodetic or not).

[ajfriend]

Sorry, I was carrying over "distortion vertices" from H3. And I totally get the distinctions between spherical and ellipsoidal geodesic curves, the underlying ellipsoid models, and the issues around converting between a spherical model and WGS84. The distortion vertices actually end up being kind of orthogonal to all that, and they exist in all three systems: H3, ISEA-based, and IVEA-based.

To clarify, by "distortion vertex", I mean any point along a cell edge where the tangent changes discontinuously, excluding the "obvious" intended vertices. This allows for segment arcs to be spherical or ellipsoidal geodesics (or some other curve) as long as they're smooth. Visually, these are the "kinks" in cells that make hexagons not look like hexagons, or H3 pentagons with 10 vertices.

You capture it correctly for H3 in your paragraph above with "one or two great circle arc segments". The edges made of two great circle arcs have a distortion vertex joining them somewhere along the cell edge. By my definition of distortion vertex, I'm claiming that ISEA3H and IVEA3H have these, too.

Here's a plot (courtesy of Claude) that demonstrates:

The red Xs are the distortion vertices; the turn-angle panels show their size, about 28 degrees in ISEA3H and about 3 degrees in IVEA3H. In ISEA you can see the kinks directly in the boundary. The improved IVEA projection mitigates them to pretty much imperceptible, at least to my eye, but they're still there.

Those change points are what I'd love to get from the API directly: vertices plus the points where the boundary direction jumps.

The current getZoneRefinedWGS84Vertices() mixes two things: 1) sampling points along differentiable arcs that make up a zone edge, and 2) defining the endpoints of those smooth arcs. I think it would be nice to query those separately. One of them changes with increased refinement, and one of them doesn't.

[jerstlouis]

@ajfriend Yes the ISEA projection in particular has those cusps, and IVEA mitigates this to almost imperceptible visually (nice to see the actual stats here!!).

Figuring out where they land is not particularly simple or easy -- it's related to the particular choice of projection (RTEA / ISEA / IVEA) and where the zone boundaries fall within the icosahedron triangles. In the 3H grids they happen halfway, but in 7H they are not halfway within the edge (e.g., see here).

getZoneRefinedWGS84Vertices() is kind of a black box giving you reasonable geometry, and exactly how this is done is subject to improvements (like trying to find these higher distortion points to refine more where needed).

Currently we don't have a method specifically for what you're asking for -- I think that would require some thinking in terms of what would be most useful along the lines of what you're suggesting.

[ajfriend]

Makes sense. I definitely appreciate that this is easier to ask for than to implement :)

So no expectations on my end. As I'm using the library, I'm just trying to note things like this in the hope that it's some help. And happy to chat more about this in the future if it'd be helpful.

[jerstlouis]

Thanks @ajfriend, I definitely appreciate all the detailed testing and feedback.

I would be curious as to learn more about the actual higher level use case i.e., why you would like these specific distortion points as separate vertices? What is the problem that would solve or what improvements would that enable?

[ajfriend]

So glad you asked! :)

I'm working on a new evaluation metric for DGGSs that measures the angular distortion of grids via defining an aspect ratio for polygons (cells) on the sphere (and maybe eventually the ellipsoid): https://github.com/ajfriend/csar_py/

I'm using that to compare grids at https://ajfriend.com/dggs_compare/. The IVEA grids currently have the best worst-case aspect ratio across all the systems I'm comparing! I'm not yet explicitly handling the sphere vs. ellipsoid differences, but I don't expect that to change the AR numbers much (it has a much bigger effect on area distortion; @MrBenGriffin has looked at this already and I still need to catch up).

To build those stats, I'm randomly sampling points on the sphere and converting them to DGGS cells. I arbitrarily chose 1 million cells per resolution as my target, which is what ended up stress testing getZoneFromWGS84Centroid().

In some cases, the distortion vertices / kinks / cusps can actually bow out significantly enough to slightly change the AR value of a cell: ajfriend/dggs_compare#25 (comment) So I need those vertices to accurately measure the aspect ratio of all the cells. I can densify the edge, but that results in a lot more vertices. Ideally, I'd just be able to query the distortion vertices directly. (Theoretically, the smooth edge arcs themselves could bow out enough to change the AR, but it seems unlikely, and I haven't seen it in practice.)

I'm running these jobs on GitHub Actions, and they take a while, so I'm looking for any opportunity to optimize. I was running out of memory on one job, which prompted #22 (which I still need to respond to). And the ecrt.Instance.delete(arr) workaround seems to increase the runtime significantly, but I still need to look into it more.