rHEALPix projection: discrepancy between pyproj/PROJ and native Python implementation in rhealpixdggs

[ndemaio]

Hello! Manaaki Whenua Landcare Research maintains rhealpixdggs, a Python package that implements the HEALPix and rHEALPix projections for use in a DGGS (repository available at https://github.com/manaakiwhenua/rhealpixdggs-py ).

We're aiming to modernise that package and would like to move away from the internal versions of the HEALPix and rHEALPix projections. We'd like to use those provided by pyproj and the PROJ libraries instead.

We've noticed some discrepancies in cell coordinates when projecting a planar dart cell in a DGGS with 3 subcells per side onto the sphere. It's causing some of our tests to fail when using pyproj/PROJ, highlighting an edge case around the 180/-180 longitudinal boundary.

A code snippet that helps check for the discrepancy is:

from rhealpixdggs.dggs import WGS84_003

rdggs = WGS84_003
c = rdggs.cell(["N", 6])
vertices = c.vertices(plane=False)

Examining vertices when using the Python implementation of rHEALPix in rhealpixdggs gives:

[(-180.0, 74.39069094862855),
 (-150.0, 41.87385774189841),
 (-180.0, 41.87385774189841),
 (149.99999999999997), 41.873857741898405)]

Examining vertices when using the rHEALPix projection provided by pyproj in rhealpixdggs gives:

[(179.99999999999994, 74.39069094862855),
 (-150.0, 41.87385774189841),
 (-180.0, 41.87385774189841),
 (149.99999999999997, 41.873857741898405)]

Note the discrepancy with inverted sign in vertices[0][0] (longitude of northernmost vertex in dart).

Also note that there's no such discrepancy in vertices[0][2] (longitude of non-vertex in dart).

We're not sure how to handle this but it's currently blocking progress on the migration to pyproj/PROJ. Any ideas on what the best way forward might be?

Some in-depth info for delving further into the problem

What gets called internally when using pyproj/PROJ to project a planar vertex (u, v) of cell N6 onto the sphere is

rdggs.rhealpix(u, v, inverse=True, region='north_polar')

which then calls (with rhealpixdggs.projection_wrapper as pw)

f = pw.Projection(
            ellipsoid=self.ellipsoid,
            proj='rhealpix',
            north_square=0
            south_square=0
            region=region, # 'north_polar'
        )
f(u, v, inverse=inverse)

with pw.Projection being defined as

pyproj.Proj(proj=proj, a=a, e=e, **kwargs)

where a and e are parameters taken from the WGS84 ellipsoid.

Environment Information

pyproj info:
pyproj: 3.6.1
PROJ: 9.3.0
data dir: /home/demaion/dev/rhealpixdggs-py/.venv/lib/python3.11/site-packages/pyproj/proj_dir/share/proj
user_data_dir: /home/demaion/.local/share/proj
PROJ DATA (recommended version): 1.15
PROJ Database: 1.2
EPSG Database: v10.094 [2023-08-08]
ESRI Database: ArcGIS Pro 3.1 [2023-19-01]
IGNF Database: 3.1.0 [2019-05-24]

System:
python: 3.11.5 (main, Aug 29 2023, 09:48:00) [GCC 11.4.0]
executable: /home/demaion/dev/rhealpixdggs-py/.venv/bin/python
machine: Linux-5.15.153.1-microsoft-standard-WSL2-x86_64-with-glibc2.35

Python deps:
certifi: 2024.6.2
Cython: None
setuptools: 65.5.0
pip: 24.1.1

[jjimenezshaw]

Hi.

I guess this is a question for PROJ. (I never know what is better, the mailing list https://lists.osgeo.org/mailman/listinfo/proj or github). Pyproj is a python layer on top of PROJ, and apparently your problem is related to the transformation itself.

About the examples you are providing, please make an example that does not depend on any external software. Yes, you are using rhealpixdggs.dggs, but you should create a test only with PROJ commands. Probably cs2cs or cct.

[ndemaio]

Thanks for the quick reply! I'll take the question to PROJ then.

I think it's safe to close this topic.