Failure to compute cell boundary in geographic coordinates

[aaime]

The following simple test script fails with an error:

from rhealpixdggs import dggs, ellipsoids
from rhealpixdggs.ellipsoids import Ellipsoid
from rhealpixdggs.dggs import RHEALPixDGGS, Cell

WGS84_TB16 = Ellipsoid(a=6378137.0, b=6356752.314140356, e=0.0578063088401, f=0.003352810681182, lon_0=-131.25)
dggs = RHEALPixDGGS(ellipsoid=WGS84_TB16, north_square=0, south_square=0, N_side=3)
cell = dggs.cell(('Q', 6, 6, 6, 6, 6))
print(*cell.vertices(False))

resulting in:

Error: input coordinates (-1.570796, 0.791862) are out of bounds
Traceback (most recent call last):
  File "/tmp/test.py", line 8, in <module>
    print(*cell.vertices(False))
  File "/home/aaime/.local/lib/python3.6/site-packages/rhealpixdggs/dggs.py", line 2245, in vertices
    result = [self.rdggs.rhealpix(*p, inverse=True) for p in result]
  File "/home/aaime/.local/lib/python3.6/site-packages/rhealpixdggs/dggs.py", line 2245, in <listcomp>
    result = [self.rdggs.rhealpix(*p, inverse=True) for p in result]
  File "/home/aaime/.local/lib/python3.6/site-packages/rhealpixdggs/dggs.py", line 453, in rhealpix
    return f(u, v, inverse=inverse)
  File "/home/aaime/.local/lib/python3.6/site-packages/rhealpixdggs/projection_wrapper.py", line 121, in __call__
    lam, phi = f(u, v, radians=radians, inverse=True)
  File "/home/aaime/.local/lib/python3.6/site-packages/rhealpixdggs/pj_rhealpix.py", line 509, in f
    x, y, e=e, north_square=north_square, south_square=south_square
TypeError: iteration over a 0-d array

Hope I have the right version of rHEALPix, installed it using:

sudo pip3 install --upgrade rhealpixdggs

and then run the test script using python3.

[aaime]

The same issue happens with other cells, e.g., "R88888"

[rggibb]

I'll chase this, it seems to work at higher resolutions and then bails when the fifth level is reached. This may be related to something that Matt mentioned that his team fixed in the GA's version, so I'll check what they did in their repo.

[geofizzydrink]

Yes. This is the same issue identified in AusPIX. It appears to be due to a projection issue at the boundary between the equatorial and polar cells (where the equal area projection switches). I've tested this at multiple resolutions (down to level 15) and the error occurs only in the northern most vertices (for the south polar cell) and the southern most vertices (for the north polar cell).

The work around for this bug that was implemented by GA is as follows:

1. wrap the "cell.vertices" function in a "try/except" statement to catch the error.
2. if the exception is triggered a) for the South Polar Cell: a.1) get the immediate north neighbour cell and get its southernmost vertices. a.2) replace the coordinates of the northernmost vertices of the cell that triggered the error with the vertex coordinates obtained from a.1). b) for the North Polar Cell: b.1) get the immediate south neighbour cell and get its northernmost vertices. b.2) replace the coordinates of the southernmost vertices of the cell that triggered the error with the vertex coordinates obtained from b.1).

NB: This is not a proper fix for the issue - merely a work around.

A proper fix will require an examination of the polar projection algorithm used by rHealPIX and determine why it is getting an "out of bounds" error when converting from planar (cube face) coordinates to global (lon/lat) coordinates. It could be a boundary computation issue, which might be solved by extending the projection boundary slightly. But the work hasn't been done yet.

Cheers, Matt.

[aaime]

I can confirm those two cells compute properly now. But only after the source code is updated as described in https://github.com/manaakiwhenua/rhealpixdggs-py/issues/10

[alpha-beta-soup]

Sample given currently results in:

>>> from rhealpixdggs import dggs, ellipsoids
>>> from rhealpixdggs.ellipsoids import Ellipsoid
>>> from rhealpixdggs.dggs import RHEALPixDGGS, Cell
>>> WGS84_TB16 = Ellipsoid(a=6378137.0, b=6356752.314140356, e=0.0578063088401, f=0.003352810681182, lon_0=-131.25)
>>> dggs = RHEALPixDGGS(ellipsoid=WGS84_TB16, north_square=0, south_square=0, N_side=3)
>>> cell = dggs.cell(('Q', 6, 6, 6, 6, 6))
>>> print(*cell.vertices(False))
(np.float64(-131.25), np.float64(-41.51722626942552)) (np.float64(-130.87962962962962), np.float64(-41.51722626942552)) (np.float64(-130.87962962962965), np.float64(-41.937853910786615)) (np.float64(-131.25), np.float64(-41.937853910786615))