4

The projection named ISEA - Icosahedral Snyder Equal Area has a (supposedly) good implementation in the PROJ Library, proj.org/isea, that is used by PostGIS, QGIS, etc.
The example is agnostic - was expressed in standard SQL with no software dependence.

Problem: bad results, seems that I am not using correct resolution or other parameter in the ISEA definition.

I tried it first with +proj=isea, but the result seems wrong (PROJ Library bug?). See "use case" below. I changed to +proj=isea +resolution=1.0 and other values (0.1, ..., 40.0) and no effect, no changes in the result.

Supposing that PROJ Library has no bug, the question is "How to use ISEA?"

Use case

SELECT ST_Area(box_geom,true) as area, -- in meters
       ST_Area( ST_Transform(box_geom,isea_proj) ) as isea_area 
FROM (
   SELECT '+proj=isea', -- +resolution=1.0
   ST_SetSRID(ST_GeomFromGeoHash('6gyf4bf'), 4326)
) t1(isea_proj,box_geom);

Result: area=21326 and isea_area=21424, so a significative difference of 0.46%.

Note: calculationg percentual difference by abs(area-proj_area)/area.

A change in projection can't affect area, only when PostGIS do some internal calculation error. The usual in projection-change is something as ~0.0001% or less. The ~0.5% result is big!

PS: the used Geohash, 6gyf4bf, is a box located at Brazil/São Paulo. Using other box, at Equador/Quito (6rbntru) or Chile/Terra del Fuego (4qqhz1e) the error is the same.

The question ended here... Below a complement to remember me the other checks that I've done.

A complete use case to compare

Bellow I will compare with another equal-area projection, valid for Brazilian territory... It shows that PostGIS is working fine, and a good projection not affects area calculations.

Check a landmark point of GeoURI geo:-23.550385,-46.633956 and a small box containing it, the geometry of the Geohash 6gyf4bf... Comparing results of ISEA and Albers equal-area projection.

  • test1 The area of the transformed boxes (rectangles of the Geohashes 6gyf4bf and 6gyf4bg). Expected near same box areas.

  • test2 The distance between box centers.

In the test1 of areas, the expected is ~0.0001% or less in area difference when comparing with PostGIS spheroid calculation with WGS84. Test2 is to check distance behavior in both projections.

SELECT geohash,
       area1, 
       st_area(box1_geom_albers) area1_albers,
       st_area(box1_geom_isea) area1_isea,

   area2,
   st_area(box2_geom_albers) area2_albers,
   st_area(box2_geom_isea) area2_isea,

   d,
   st_distance(pt1_geom_albers, pt2_geom_albers) as d_albers,
   st_distance(pt1_geom_isea, pt2_geom_isea) as d_isea 



FROM (
 SELECT
   st_geohash(pt1_geom,8) as geohash, -- 6gyf4bf1
   st_area(box1_geom,true) as area1,
   st_area(box2_geom,true) as area2,


st_transform(pt1_geom,952019) pt1_geom_albers,
   st_transform(pt2_geom,952019) pt2_geom_albers,
   st_transform(box1_geom,952019) box1_geom_albers,
   st_transform(box2_geom,952019) box2_geom_albers,


st_transform(pt1_geom,isea_proj) pt1_geom_isea,
   st_transform(pt2_geom,isea_proj) pt2_geom_isea,
   st_transform(box1_geom,isea_proj) box1_geom_isea,
   st_transform(box2_geom,isea_proj) box2_geom_isea,
   st_distance(pt1_geom, pt2_geom, true) d
 FROM (
   SELECT '+proj=isea', -- +resolution=1.0
   ST_SetSRID(ST_GeomFromGeoHash('6gyf4bf'), 4326),
   ST_SetSRID(ST_GeomFromGeoHash('6gyf4bg'), 4326),
   ST_SetSRID( ST_Point( -46.633889, -23.550278), 4326), -- landmark
   ST_SetSRID( ST_Point( -46.632156, -23.549881), 4326)  -- centroid of 6gyf4bg
 ) t1(isea_proj,box1_geom, box2_geom, pt1_geom, pt2_geom)
) t2;

Results:

area1 area1_albers area1_isea area2 area2_albers area2_isea
21326.147181987762 21326.14718184221 21424.010447963992 21326.147181987762 21326.14718169578 21424.01044831646
d d_albers d_isea
182.32446429 183.18969403271427 191.02939743271384

  • Area Albers difference: 0.000000145552 = abs(21326.1472 - 21326.14719); so 0% of 21326.

  • Area ISEA difference: 97.863265976230 = abs(21326.1472 - 21424.01045); so ~0.5% of 21326.

  • Distance Albers difference: ~0.47% = 100*abs(182.32446429-183.18969403271427)/182.32446429

  • Distance ISEA difference: ~4.78% = 100*abs(182.32446429-191.02939743271384)/182.32446429

So, ISEA area error is 10^7 times bigger, and distance error ~10 ... Something wrong in ISEA transform.

SRIDs on use case:

INSERT INTO spatial_ref_sys (srid, auth_name, auth_srid, proj4text) VALUES
(  -- Brazilian Albers of 2019, SRID:
  952019, 'BR:IBGE', 52019,
  '+proj=aea +lat_0=-12 +lon_0=-54 +lat_1=-2 +lat_2=-22 +x_0=5000000 +y_0=10000000 +ellps=WGS84 +units=m +no_defs'),

(  -- ISEA:
  953000, 'INT:ISEA', 953000, '+proj=isea')


ON CONFLICT DO NOTHING;
alphabetasoup
  • 8,718
  • 4
  • 38
  • 78
Peter Krauss
  • 2,292
  • 23
  • 43
  • For the result to be comparable, you should calculate the spherical area (not ellipsoidal) in geographic coordinates. You could also define the Lambers conic equiareal projection from a spherical surface as well, and compare the three results. – Gabriel De Luca Dec 11 '21 at 15:27
  • Hi @GabrielDeLuca, good clues, it makes sense! Please help us expanding it as an answer-text, expressing it more mathematically, or changing proj4text strings. – Peter Krauss Dec 11 '21 at 21:12

2 Answers2

1

This is a Wiki (!), please collabore editing here... Or adding a correct and complete answer!

Trying to answer

The problem in "how to use" was the area calculation, so, the alternative answers (please check if wrong) can be:

  • Datum surface of ISEA is a perfect sphere, and WGS84's is aoblate spheroid: the differences in the model lead to differences in area calculation. Test1 below shows it, and how to correct it.

  • Check how to calculate ISEA area using steradians instead square meter. Perhaps is possible to change PROJ string to obtain coordinates in spherical coordinates (polar and azimuthal angles).

Evidence tests:

Test1

((A first guess, but it is not a complete answer - please comment or suggest other test.))

Using the area-altitude correction function, areakm_at_altitude(), offered by this other question-answer:

SELECT id, gname, altitude,
   ROUND(isea_area::numeric,3) AS isea_area,
   ROUND(area_at_altitude::numeric,3) AS area_at_altitude,
   ROUND((100.0*(area_at_altitude-isea_area)/isea_area)::numeric,2)::text || '%' AS iseaalt_diff_perc,
   ROUND(st_area_km::numeric,3) AS st_area_km,
   ROUND((100.0*(area_at_altitude-st_area_km)/isea_area)::numeric,2)::text || '%' AS isea_diff_perc
FROM (
  SELECT id, gname, altitude,
       ST_Area( ST_Transform(geom,'+proj=isea') )/1000 AS isea_area,
       areakm_at_altitude(geom, altitude) AS area_at_altitude,
       ST_Area(geom::geography)/1000 AS st_area_km
  FROM test_table
) t;

Results:

id gname altitude ISEA_area area_at_altitude ISEAalt_diff_perc st_area_km ISEA_diff_perc
1 (fromText) 0 21.788 21.689 -0.46% 21.689 0.00%
2 (fromText) 29.1704 21.788 21.788 0.00% 21.689 0.46%
3 (fromText) 800 21.788 24.418 12.07% 21.689 12.53%
10 6gyf4bf 0 21.424 21.326 -0.46% 21.326 0.00%
11 6gyf4bf 29.1704 21.424 21.424 0.00% 21.326 0.46%
20 6gzx9 0 22099.150 21994.555 -0.47% 21994.555 0.00%
21 6gzx9 29.1704 22099.150 22095.253 -0.02% 21994.555 0.46%
30 gbsuv 0 15804.329 15817.722 0.08% 15817.722 0.00%
31 gbsuv 29.1704 15804.329 15890.307 0.54% 15817.722 0.46%

In bold values are areas where altitude correction was applied, when ISEA_area and area_at_altitude area near equal (with ISEAalt_diff_perc near 0%), both are bold.

Test1 conclusions: seems that ISEA sphere is used as surface of area modeling, instead ellipsoid. For near locations (gname column with values "fromText", 6gyf4bf and 6gzx9) the altitude correction make sense, and all can use the same correction. For other locations (gbsuv) other altitude correction must be used.

Sample data


DROP TABLE test_table;
CREATE TABLE test_table AS
 SELECT 1 AS id, '(fromText)' as gname, 
        0.0::float AS altitude,
        'POLYGON'::text AS gtype,
        ST_GeomFromText(
         'POLYGON((-46.6342 -23.55057,-46.6342 -23.5492,-46.6328 -23.5492,-46.6328 -23.55057,-46.6342 -23.55057))',
         4326
       ) AS geom

UNION ALL


SELECT 2, '(fromText)', 29.1704, 'POLYGON Z', ST_GeomFromText(
   'POLYGON Z((-46.6342 -23.55057 29.1704,-46.6342 -23.5492 29.1704,-46.6328 -23.5492 29.1704,-46.6328 -23.55057 29.1704,-46.6342 -23.55057 29.1704))',
   4326
 )


UNION ALL


SELECT 3, '(fromText)', 800.0, 'POLYGON Z', ST_GeomFromText(
   'POLYGON Z((-46.6342 -23.55057 800.0,-46.6342 -23.5492 800.0,-46.6328 -23.5492 800.0,-46.6328 -23.55057 800.0,-46.6342 -23.55057 800.0))',
   4326
 )


UNION ALL


SELECT 10, '6gyf4bf', 0.0,     'POLYGON', ST_SetSRID(ST_GeomFromGeoHash('6gyf4bf'), 4326)
UNION ALL
 SELECT 11, '6gyf4bf', 29.1704, 'POLYGON', ST_SetSRID(ST_GeomFromGeoHash('6gyf4bf'), 4326)


UNION ALL


SELECT 20, '6gzx9', 0.0,     'POLYGON', ST_SetSRID(ST_GeomFromGeoHash('6gzx9'), 4326)
UNION ALL
 SELECT 21, '6gzx9', 29.1704, 'POLYGON', ST_SetSRID(ST_GeomFromGeoHash('6gzx9'), 4326)


UNION ALL


SELECT 30, 'gbsuv', 0.0,     'POLYGON', ST_SetSRID(ST_GeomFromGeoHash('gbsuv'), 4326)
UNION ALL
 SELECT 31, 'gbsuv', 29.1704, 'POLYGON', ST_SetSRID(ST_GeomFromGeoHash('gbsuv'), 4326)
;
Peter Krauss
  • 2,292
  • 23
  • 43
1

A map projection is a 2D transformation (in mathematical terms; in cartographic terms it is a conversion), from one reference surface to a plane (or to a developable surface towards a plane).

For different projections to be comparable, the reference surface must be the same.

Supposing that PROJ Library has no bug, the question is "How to use ISEA?"

The behavior you observe is not a bug in the PROJ library.

The ISEA projection, by design, only supports a spherical reference surface.

To demonstrate this, you just have to compare it with other projections using a spherical reference surface, for example the Normal Sphere (+ellps=sphere as PROJ built-in ellipsoid definition):

INSERT INTO spatial_ref_sys (srid, auth_name, auth_srid, proj4text) VALUES
(  -- Normal Sphere Geographic Coordinates:
  953001, 'USER', 953001, '+proj=lonlat +ellps=sphere'),
(  -- Normal Sphere ISEA:
  953002, 'USER', 953002, '+proj=isea +ellps=sphere'),
(  -- Normal Sphere Brazilian Albers:
  953003, 'USER', 953003,
  '+proj=aea +lat_0=-12 +lon_0=-54 +lat_1=-2 +lat_2=-22 +x_0=5000000 +y_0=10000000 +ellps=sphere +units=m +no_defs')
ON CONFLICT DO NOTHING;

WITH data AS (
  SELECT 
    ST_SetSRID(ST_GeomFromGeoHash('6gyf4bf'), 953001) AS box_geom
)
SELECT
  ST_Area(box_geom, true) AS area,
  ST_Area(ST_Transform(box_geom, 953002)) AS isea_area,
  ST_Area(ST_Transform(box_geom, 953003)) AS albers_area
FROM
  data;
area isea_area albers_area
21376.0711152554 21376.0711155145 21376.0711144206

If you want the area projected from the sphere to be the same as projected from an ellipsoid, that is not possible.

There may be some approaches. But keep in mind that, except at specific points, the ellipsoid's radius of curvature varies in all directions.

It occurs to me that you could calculate the mean between Meridional and Prime vertical sections radii of curvature for a central point of the polygon to be projected, and use that radius.

But if the question is whether that is the way to use the ISEA projection: I don't think so.

Gabriel De Luca
  • 14,289
  • 3
  • 20
  • 51
  • Thanks Gabriel for good explanations (!) and to show how to use +ellps=sphere... It seems that it is the "Rosetta Stone" to interoperate and to compare with the WGS84 standard (SRID 4326) and others. – Peter Krauss Dec 24 '21 at 12:36
  • You are welcome. Can use +R=xxx parameter instead of +ellps=xxx too. As you saw in your question, they can be interoperable as long as you assume the error product of using a different spheroid. – Gabriel De Luca Dec 24 '21 at 12:41