Scale of (r)HEALPix coordinate space

[jerstlouis]

AI usage

• Google AI supported the creation of this ticket.
• This ticket has been (mostly) written by myself (human), and is aimed at being read by humans.

Example of problem

See AddressForAll/suporte#86 for PostGIS output illustrating the ~15.12% discrepancy from ST_Area()

Problem description

The current implementation of rHEALPix scales the radians coordinates (x ranging from -pi .. pi) by the semi-major axis. In order for the total area of the coordinate space to correspond to the true ellipsoid area, I believe this should be semi-major axis X sqrt(8÷(3×pi)). Alternatively, the x and y could be pre-multiplied by that factor in e.g. healpix_sphere.

Note that this likely also applies to HEALPix.

I'm not sure whether this breaking change would be desired by the community, but at minimum this aspect should be properly documented, so that users know that the coordinates returned by PROJ are essentially scaled by a ~0.9213 factor compared to true meters (8÷(3×pi) = ~0.8488 factor for areas).

cc. @rggibb @ppKrauss @tinaok @tik65536 @allixender

[jerstlouis]

@rouault @kbevers @jjimenezshaw

Do you know if in general the units should be true meters so that ST_Area() on equal-area projections returns real areas?

https://proj.org/en/stable/usage/projections.html#projection-units mentions:

Scaling of output units can be done by applying the +k_0 argument (unitless).

So perhaps (for HEALPix and rHEALPix) adding support for the +k_0 parameter here could help with backward compatibility with the current scaled units.
(i.e., fix the units to be true meters, and document that +k_0=0.92131773... -- sqrt(8÷(3×pi)) needs to be used to reproduce the old behavior).

[kbevers]

Do you know if in general the units should be true meters so that ST_Area() on equal-area projections returns real areas?

I can't answer this with much confidence but it makes sense to me that a true meters would result in correct area determinations. Wether you can achieve true meters across all areas of a projection is another question. If reliable area determinations are needed it's better to do them using geodetic coordinates on the ellipsoid using @cffk's algorithms.

[jerstlouis]

@kbevers Thanks. Right, my question was not whether ST_Area() would result in correct areas or not -- this should be the case for an equal-area projection, provided the polygon in question does not cross an interruption, and that the polygon contours are segmented fine enough.

My question was about the units of the projection. In the case of this issue here, both distances and areas would be off by a constant factor (0.9213 for distances, 0.9213^2 for areas). So I am trying to understand whether this is normal behaviors for some projections, or something that should really be fixed.

[jjimenezshaw]

My question was about the units of the projection.

For an equal area projection, it makes sense to me to have it scaled to make them equal (but it is just an opinion).

You are suggesting to use +k_0 parameter. Isn't it always available? (so you can use it without changing the code)

[jerstlouis]

@jjimenezshaw

+k_0 parameter. Isn't it always available?

That might be the case if it is handled centrally outside the projection, but I was suggesting it as a backward compatibility option after "fixing" the scale (changing the current behavior).

[ppKrauss]

Hi, I understand that PROJ's rHEALPix implementation today is not useful and has a bug. It need a bug fix... Is it?

The official explanation says "To run a rHEALPix projection (on the default GRS80 ellipsoidal model) use ... +proj=rhealpix +north_square=2". It is false, there are no "ellipsoidal model".

---

Only a WGS84-compatible projection is really useful. The parameter in +proj=rhealpix +ellps=WGS84 has no appreciable effect... Even if it did have an effect, as @jerstlouis demonstrated, the actual difference is approximately 15%.

The actual surface area of the WGS84 ellipsoid is ~510,065,621.72 km² , and 510,065,621.72÷602,252,610 = 0.846929699, 1 - 0.846929699 = ~15.3% which likely explains that discrepancy.

Script showing the bug

The bug of the discrepancy 15% in a supposed WGS84 projection.

INSERT INTO spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext) VALUES
( -- rHEALPix variant 2, PROJ v4.8+
  955002,  'DGGS:002:rHEALPix',  2,
  '+proj=rhealpix +ellps=WGS84 +south_square=0 +north_square=2',
  NULL -- no srtext
) ON CONFLICT DO NOTHING;

CREATE or replace FUNCTION  round(float,int) RETURNS numeric
  IMMUTABLE PARALLEL SAFE BEGIN ATOMIC; SELECT ROUND($1::numeric,$2);
END;

CREATE VIEW comparing_projs AS
 SELECT ghs, city_name,  t6.*
 FROM generate_series(2,4) t0(i), 
    LATERAL ( VALUES
      (ST_GeoHash( ST_Point(-46.633, -23.550, 4326), i ), 'Sao Paulo' ), 
      (ST_GeoHash( ST_Point(-60.671, 2.8194, 4326), i ), 'Boa Vista' ), 
      (ST_GeoHash( ST_Point(-34.90, -8.050, 4326), i ), 'Recife' ), 
      (ST_GeoHash( ST_Point(-51.23, 30.0331, 4326), i ), 'Porto Alegre' ), 
      (ST_GeoHash( ST_Point(-67.81, -9.9747, 4326), i ), 'Rio Branco' )
    ) t1 (ghs,city_name),
    LATERAL (SELECT ST_SetSRID( ST_GeomFromGeoHash(t1.ghs), 4326) AS geom) t2,
    LATERAL (SELECT ST_Area(t2.geom,true) AS wgs84_area_m2) t3,
    LATERAL (SELECT SQRT(t3.wgs84_area_m2) AS side_m) t4,
    LATERAL (SELECT ST_Area(ST_Transform(t2.geom,955002)) AS rhpx_m2) t5,
    LATERAL (SELECT
      round( wgs84_area_m2/1000^2, 2) as wgs84_area_km2,
      round( rhpx_m2/1000^2, 2) as rhpx_km2,  
      wgs84_area_m2/rhpx_m2 - 1.0    AS rhpx_diff  -- m2 change
    ) t6
 ORDER BY 1,2;
----
SELECT ghs, city_name,
               wgs84_area_km2, rhpx_km2,
               round(100.0*rhpx_diff, 2)::text||'%' AS rhpx_chg
FROM comparing_projs;

Results:

ghs	city_name	wgs84_area_km2	rhpx_km2	rhpx_chg
6g	Sao Paulo	706624.09	831256.86	-14.99%
6gy	Sao Paulo	22419.39	26411.50	-15.12%
6gyf	Sao Paulo	699.25	823.78	-15.12%
6q	Rio Branco	772756.87	907633.13	-14.86%
6qp	Rio Branco	23941.67	28204.33	-15.11%
6qpz	Rio Branco	749.58	883.07	-15.12%
7n	Recife	772756.87	907633.13	-14.86%
7nx	Recife	24126.64	28422.17	-15.11%
7nx4	Recife	753.48	887.66	-15.12%
d8	Boa Vista	780217.06	916225.26	-14.84%
d8s	Boa Vista	24297.71	28623.63	-15.11%
d8sb	Boa Vista	759.73	895.03	-15.12%
dv	Porto Alegre	670664.67	789577.81	-15.06%
dv7	Porto Alegre	21101.85	24859.72	-15.12%
dv7d	Porto Alegre	661.17	778.92	-15.12%

Other simple test

In an equatorial or tropical zones (like Brazilian cities of the test excep Porto Alegre) compare rHealPix with the "Lambert rHEALPix face", that is, the Lambert cylindrical equal-area projection — in PROJ string is +proj=cea +ellps=WGS84 (see cea).

The result is the same ~15% of discrepancy.

[jerstlouis]

Thanks for chiming in @ppKrauss .

is not useful

I wouldn't go that far in saying it's not useful.
It produces a useful truly equal-area projection, but the units of the result are scaled relative to true meters.

It is false, there are no "ellipsoidal model".

I am quite sure it does use the GRS80 ellipsoidal model by default. See https://github.com/OSGeo/PROJ/blob/master/src/projections/healpix.cpp#L618 . It checks the meridional eccentricity (es) and if it's non-zero it calls e_healpix_forward/inverse.

The parameter in +proj=rhealpix +ellps=WGS84 has no appreciable effect...

That is because WGS84 is practically identical to GRS80. The only difference is the semi-minor axis is 6356752.314140 instead of 6356752.314245, so we're talking about a 0.1 mm difference.

the actual difference is approximately 15%.

The issue here is unrelated to the ellipsoid, but an overall scaling issue when mapping radians to meters.
The same 15% difference applies to a spherical Earth model as well.
The spherical vs. ellipsoid difference would have been much less than 15%.