How to define a local projection with correct altitude for correct area-metrics

Question

Using PROJ, and defining custom projection, by its PROJ-strings.

How to correct WGS84 altitude, to my "fit to local" reference?

---

For example this Albers projection for Brazil:

+proj=aea +lat_0=-12 +lon_0=-54 +lat_1=-2 +lat_2=-22 +x_0=5000000 +y_0=10000000 +ellps=WGS84 +units=m

Can I add a z parameter, or something else, to correct the altitude?

It is using the WGS84 as Z-reference, but I need to apply a "correction factor" to increase the area measurements — because the WGS84 ellipsoid not fits the "median Earth's surface altitude of Brazilian territory" (illustrated).

... The projection (and coordinate) system must to represent conditions “at ground”, in the local zone.

Suppose, for instance, the correct median altitude is 600 meters above sea level.
PS: "Earth's surface Brazilian sea level", not false sea level of ellipsoid.

---

Notes

As commented before,

area increase when altitude increase (Earth radius).

Tests for answers

a	a19	a23
21688.8155010	21688.8154995	?

SELECT st_area(geom::geography) as a,
       st_area(g19) as a19, st_area(g23) as a23
FROM (
 SELECT *, ST_Transform(geom,952019) as g19, ST_Transform(geom,952023) as g23
 FROM (                    
   SELECT
          ST_GeomFromText(
           'POLYGON((-46.6342 -23.55057,-46.6342 -23.5492,-46.6328 -23.5492,-46.6328 -23.55057,-46.6342 -23.55057))',
           4326
          )
 ) t1(geom)
) t2;

To use with QGIS need a SRID, not string, so before:

INSERT INTO spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext) VALUES
   -- -- -- --
(  -- IBGE Albers, SRID number convention in my Project:
  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',
  -- GCS_SIRGAS2000 is, for practical applications, same as WGS84 
  $$PROJCS[
  "Conica_Equivalente_de_Albers_Brasil",
  GEOGCS[
    "GCS_SIRGAS2000",
    DATUM["D_SIRGAS2000",
       SPHEROID["Geodetic_Reference_System_of_1980",6378137,298.2572221009113]
    ],
    PRIMEM["Greenwich",0],
    UNIT["Degree",0.017453292519943295]
  ],
  PROJECTION["Albers"],
  PARAMETER["standard_parallel_1",-2],
  PARAMETER["standard_parallel_2",-22],
  PARAMETER["latitude_of_origin",-12],
  PARAMETER["central_meridian",-54],
  PARAMETER["false_easting",5000000],
  PARAMETER["false_northing",10000000],
  UNIT["Meter",1]
 ]$$
),
(  --  Albers2, SRID:
  952023, 'test', 52023,
  '...', NULL
);

Answer 1

The cartographic area is always the area reduced to an ellipsoid, and in this case (a cartographic projection), after reduction to the ellipsoid the area is transformed to a plane coordinate system.

The geometries may have Z coordinates, but as long as the reference systems are two-dimensional, that third dimension does not participate in any calculation.

---

A map projection to a plane grid is not exactly as you describe in the image figure. Although there may be projections with these geometric properties, it is always better to consider them in their general form:

A map projection is a mathematical conversion in two dimensions, defined by two equations: x, which is a function of latitude and longitude, and y, which is another function of latitude and longitude.

---

The definition of coordinate systems through proj4 strings is a bit obsolete, however it still works and will continue to do so due to the number of applications that depend on it. And it is possible to define a projection from an ellipsoid different from the one assumed by default, making it explicit in the definition of the reference system.

---

However, for latitude and longitude to change from one ellipsoid to another, a datum transformation is required.

In a datum transformation, latitude and longitude are converted to global, three-dimensional Cartesian coordinates, and an affine transformation in Cartesian space can be applied. This transformation could include translation, rotation and scaling of the global Cartesian coordinates, then transformed to longitude and latitude again.

---

Taking into account this basic introduction, and your practical problem of trying an engineering task (calculating the area on the ground) through a cartographic information system (where the areas are always reduced to the ellipsoid), I propose a practical solution:

Define a coordinate system through a proj4 string that performs a Cartesian coordinate scaling transformation of the datum, and a cartographic projection from an ellipsoid larger than the one used by default.

---

The datum transformation can be defined as seven Helmert parameters towards WGS-84. The seventh parameter is the scale difference in parts per million. Assuming a sphere with a radius of 6 million meters, 100 parts per million is 600 meters.

The reference ellipsoid for the cartographic projection can be defined by its semi-major axis parameter, measuring 600 meters longer than WGS-84, but the same reverse flattening.

+proj=aea +lat_0=-12 +lon_0=-54 +lat_1=-2 +lat_2=-22 +x_0=5000000 +y_0=10000000 +a=6378737 +rf=298.257223563 +towgs84=0,0,0,0,0,0,-100 +units=m +no_defs +type=crs

---

It's a good solution? Not at all. Not so much for the solution but for the problem. It is just one possible solution, but it does not come with guarantees about side effects.

---

Tests

For those who do not find the answer so long, here are some examples to help understand these transformations:

Scale the datum but not the ellipsoid:

SELECT ST_AsText(ST_Transform( geom, 
  '+proj=latlon +ellps=WGS84 +towgs84=0,0,0,0,0,0,-100 +units=m +no_defs +type=crs'), 5) AS geom
FROM (
  SELECT 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) AS t1;
                                                   geom                                                    
-----------------------------------------------------------------------------------------------------------
 POLYGON((-46.6342 -23.55056,-46.6342 -23.54919,-46.6328 -23.54919,-46.6328 -23.55056,-46.6342 -23.55056))
(1 row)

A small change can be seen in the fifth decimal place in the latitudes. This is because the points were scaled, and from their new position the normals are drawn towards the small ellipsoid.

Do not scale the datum but the ellipsoid:

SELECT ST_AsText(ST_Transform( geom, 
  '+proj=latlon +a=6378737 +rf=298.257223563 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +type=crs'), 5) AS geom
FROM (
  SELECT 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) AS t1;
                                                   geom                                                    
-----------------------------------------------------------------------------------------------------------
 POLYGON((-46.6342 -23.55058,-46.6342 -23.54921,-46.6328 -23.54921,-46.6328 -23.55058,-46.6342 -23.55058))
(1 row)

Now, the latitudes represent the normals to the large ellipsoid that pass through the points on the small ellipsoid. So far the only thing we have proven is that latitudes are not geocentric.

Scale the datum and don't define the ellipsoid:

SELECT ST_AsText(ST_Transform( geom, 
  '+proj=latlon +towgs84=0,0,0,0,0,0,-100 +units=m +no_defs +type=crs'), 5) AS geom
FROM (
  SELECT 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) AS t1;
ERROR:  could not parse proj string '+proj=latlon +towgs84=0,0,0,0,0,0,-100 +units=m +no_defs +type=crs'
CONTEXT:  SQL function "st_transform" statement 1

For better or worse, my postgis doesn't like that. Scaling the datum is not the problem, the problem is having defined a datum transformation without an ellipsoid.

Scale the datum and the ellipsoid:

SELECT ST_AsText(ST_Transform( geom, 
  '+proj=latlon +a=6378737 +rf=298.257223563 +towgs84=0,0,0,0,0,0,-100 +units=m +no_defs +type=crs'), 5) AS geom
FROM (
  SELECT 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) AS t1;
                                                  geom                                                   
---------------------------------------------------------------------------------------------------------
 POLYGON((-46.6342 -23.55057,-46.6342 -23.5492,-46.6328 -23.5492,-46.6328 -23.55057,-46.6342 -23.55057))
(1 row)

Wait, the geographic coordinates didn't change. Thinking about it it makes sense, we're scaling the datum and the ellipsoid by about the same amount, but then... Why are we transforming the datum? Couldn't we directly project the same geometry but from the large ellipsoid? Probably yes. In that case, sorry for the noise.

The areas:

SELECT 
  ST_Area(ST_Transform(geom, '+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 +type=crs')) AS area_0,
  ST_Area(ST_Transform(geom, '+proj=aea +lat_0=-12 +lon_0=-54 +lat_1=-2 +lat_2=-22 +x_0=5000000 +y_0=10000000 +ellps=WGS84 +towgs84=0,0,0,0,0,0,-100 +units=m +no_defs +type=crs')) AS area_1,
  ST_Area(ST_Transform(geom, '+proj=aea +lat_0=-12 +lon_0=-54 +lat_1=-2 +lat_2=-22 +x_0=5000000 +y_0=10000000 +a=6378737 +rf=298.257223563 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +type=crs')) AS area_2,
  ST_Area(ST_Transform(geom, '+proj=aea +lat_0=-12 +lon_0=-54 +lat_1=-2 +lat_2=-22 +x_0=5000000 +y_0=10000000 +a=6378737 +rf=298.257223563 +towgs84=0,0,0,0,0,0,-100 +units=m +no_defs +type=crs')) AS area_3,
  ST_Area(ST_Transform(geom, '+proj=aea +lat_0=-12 +lon_0=-54 +lat_1=-2 +lat_2=-22 +x_0=5000000 +y_0=10000000 +a=6378737 +rf=298.257223563 +units=m +no_defs +type=crs')) AS area_4
FROM (
  SELECT 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) AS t1;
      area_0       |      area_1       |       area_2       |      area_3       |       area_4       
-------------------+-------------------+--------------------+-------------------+--------------------
 21688.81549992827 | 21688.80786481774 | 21692.903468480195 | 21692.89583164296 | 21692.896284356077
(1 row)

Great exercise. It seems like a datum transformation was not necessary :-)

Anyway, don't use this for anything other than calculating many areas, estimating a multiplicative factor and from there simply using it when you want to calculate a new area, always keeping the geometries in their original spatial reference system without creating a new, custom one.