EPSG:5070 mismatch between spatial_ref_sys and PROJ

[jbkoch]

I discovered a coordinate transformation issue with projection EPSG:5070 when upgrading an Amazon RDS instance from/to:

POSTGIS="3.1.5 c60e4e3" [EXTENSION] PGSQL="130" GEOS="3.9.1-CAPI-1.14.2" PROJ="Rel. 5.2.0, September 15th, 2018" GDAL="GDAL 2.4.4, released 2020/01/08" LIBXML="2.9.1" LIBJSON="0.13.1" LIBPROTOBUF="1.3.2" WAGYU="0.5.0 (Internal)" RASTER

POSTGIS="3.1.7 aafe1ff" [EXTENSION] PGSQL="140" GEOS="3.9.1-CAPI-1.14.2" PROJ="8.0.1" GDAL="GDAL 3.4.3, released 2022/04/22" LIBXML="2.9.1" LIBJSON="0.15" LIBPROTOBUF="1.3.2" WAGYU="0.5.0 (Internal)" RASTER

Running the following query on both systems produces mismatched results:

--
select
  st_astext( geom_5070 )
from st_geomfromtext(
  'POLYGON((
    -94.23933949466884 41.46010784655825,
    -94.23902022173033 41.46011128632304,
    -94.23901390824655 41.45984441537449,
    -94.23933317987697 41.45984097564519,
    -94.23933949466884 41.46010784655825))', 4326 ) geom_4326,
    st_transform( geom_4326, 5070 ) geom_5070;

POLYGON((145953.8164326792 2051852.6976375678,145980.27310072441 2051853.5729657637,145981.34988068073 2051823.7175862629,145954.89322071325 2051822.842260323,145953.8164326792 2051852.6976375678))

POLYGON((145953.2894197142 2051852.3578379445,145979.74635016784 2051853.2334285255,145980.8230537888 2051823.3784043754,145954.36613142194 2051822.5028163139,145953.2894197142 2051852.3578379445))

I was able to fix the issue by using the PROJ.4 string defined in PROJ 8.0.1:

projinfo EPSG:5070

+proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +type=crs

I'm not sure if there are other projections that need updated as well, e.g., ESRI:102039 is the same as EPSG:5070?

[robe]

jbkoch for proj 6 and above, the spatial_ref_sys proj4text is not used unless there is no corresponding auth/srid found in the PROJ proj.db.

The proj I see currently defined in spatial_ref_sys is:

+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs 

Which agree is different from your proj 8

That said this should only be an issue for older PROJs and I'm hesitant to mess with those because

a) They are all essentially deprecated by proj project and I'm not sure the type=crs even works with older proj, so may cause damage. b) For PostGIS 3.1, we do not currently push updates, so even if we fixed for 3.1, you wouldn't get the updated change anyway without manual intervention.

We didn't start shipping of allowing update for spatial_ref_sys until PostGIS 3.3.

@pramsey you think this issue is worth fixing, or should we just let dead dogs rest.

[jbkoch]

@robe Thank you for the information as I was not aware of the precedence used. This is somewhat of a special case as this is happening in the official Amazon RDS Postgresql instances of which I have no control over packaging. It would appear that something is going wrong in their dependency stack to not have the correct projection info in proj.db for PROJ 8.0.1. I am not sure who manages that project or who would be the correct person to contact about this issue?

[pramsey]

I think this is a "close and walk away" problem. The differences are sub-meter, the versions are old, the solution is… I'm not even sure what it is, so it would be quite a bit of work.

[robe]

Replying to jbkoch:

@robe Thank you for the information as I was not aware of the precedence used. This is somewhat of a special case as this is happening in the official Amazon RDS Postgresql instances of which I have no control over packaging. It would appear that something is going wrong in their dependency stack to not have the correct projection info in proj.db for PROJ 8.0.1. I am not sure who manages that project or who would be the correct person to contact about this issue?

At jbkoch, are you saying the 8.0.1 on Amazon RDS is wrong? I thought you were complaining about the proj 5.2.0 on your 3.1.5 install which would be wrong, because we ship a different proj4, and your 3.1.5 is using proj 5.2.0, which does use the proj4text in spatial_ref_sys.

If there is an issue with what AMZ is shipping, I can talk with the AMZ packaging group about that.

When you said

I was able to fix the issue by using the PROJ.4 string defined in PROJ 8.0.1

Copy this where? If you updated your spatial_ref_sys on your AMZ RDS (running PROJ 8.0.1) with this proj4text, it shouldn't have made a difference as it would still be using the proj.db shipped with PROJ, and would only fall back on spatial_ref_sys, if that entry was entirely missing in the proj.db that Amazon ships.

On my Amazon RDS database running PROJ 8.0.1:

POSTGIS="3.3.2 4975da8" [EXTENSION] PGSQL="150" GEOS="3.9.1-CAPI-1.14.2" PROJ="8.0.1" GDAL="GDAL 3.4.3, released 2022/04/22" LIBXML="2.9.1" LIBJSON="0.15" LIBPROTOBUF="1.3.2" WAGYU="0.5.0 (Internal)" TOPOLOGY RASTER

I get:

select
  st_astext( geom_5070 )
from st_geomfromtext(
  'POLYGON((
    -94.23933949466884 41.46010784655825,
    -94.23902022173033 41.46011128632304,
    -94.23901390824655 41.45984441537449,
    -94.23933317987697 41.45984097564519,
    -94.23933949466884 41.46010784655825))', 4326 ) geom_4326,
    st_transform( geom_4326, 5070 ) geom_5070;

POLYGON((145953.2894197142 2051852.3578379445,145979.74635016784 2051853.2334285255,145980.8230537888 2051823.3784043754,145954.36613142194 2051822.5028163139,145953.2894197142 2051852.3578379445))

Which corresponds to your second answer.

If I run:

    UPDATE spatial_ref_sys SET proj4text = '+proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +type=crs'
    WHERE srid=5070;

and check the results, as I would expect, I get the same answer:

POLYGON((145953.2894197142 2051852.3578379445,145979.74635016784 2051853.2334285255,145980.8230537888 2051823.3784043754,145954.36613142194 2051822.5028163139,145953.2894197142 2051852.3578379445))

This should only make a difference if you are updating your 3.1.5 instance which ships with PROJ 5.2.0 (predates the existence of the new proj api, so uses the spatial_ref_sys.proj4text)

[jbkoch]

Replying to robe:

Replying to jbkoch:

@robe Thank you for the information as I was not aware of the precedence used. This is somewhat of a special case as this is happening in the official Amazon RDS Postgresql instances of which I have no control over packaging. It would appear that something is going wrong in their dependency stack to not have the correct projection info in proj.db for PROJ 8.0.1. I am not sure who manages that project or who would be the correct person to contact about this issue?

At jbkoch, are you saying the 8.0.1 on Amazon RDS is wrong? I thought you were complaining about the proj 5.2.0 on your 3.1.5 install which would be wrong, because we ship a different proj4, and your 3.1.5 is using proj 5.2.0, which does use the proj4text in spatial_ref_sys.

Yes, this is what I believe to be the issue. The 3.1.5 install with the older proj 5.2.0 had the "correct/reference" values and when I upgraded to the 3.1.7 install with the newer proj 8.0.1, the values were off. I was able to match the old values by using the PROJ.4 string which I grabbed from a local proj 8.0.1 build (more details below).

Replying to robe:

If there is an issue with what AMZ is shipping, I can talk with the AMZ packaging group about that.

I would very much appreciate that if you too believe there is a potential issue here.

Replying to robe:

When you said

I was able to fix the issue by using the PROJ.4 string defined in PROJ 8.0.1

Copy this where? If you updated your spatial_ref_sys on your AMZ RDS (running PROJ 8.0.1) with this proj4text, it shouldn't have made a difference as it would still be using the proj.db shipped with PROJ, and would only fall back on spatial_ref_sys, if that entry was entirely missing in the proj.db that Amazon ships.

On my Amazon RDS database running PROJ 8.0.1:

POSTGIS="3.3.2 4975da8" [EXTENSION] PGSQL="150" GEOS="3.9.1-CAPI-1.14.2" PROJ="8.0.1" GDAL="GDAL 3.4.3, released 2022/04/22" LIBXML="2.9.1" LIBJSON="0.15" LIBPROTOBUF="1.3.2" WAGYU="0.5.0 (Internal)" TOPOLOGY RASTER

I get:

select
  st_astext( geom_5070 )
from st_geomfromtext(
  'POLYGON((
    -94.23933949466884 41.46010784655825,
    -94.23902022173033 41.46011128632304,
    -94.23901390824655 41.45984441537449,
    -94.23933317987697 41.45984097564519,
    -94.23933949466884 41.46010784655825))', 4326 ) geom_4326,
    st_transform( geom_4326, 5070 ) geom_5070;

POLYGON((145953.2894197142 2051852.3578379445,145979.74635016784 2051853.2334285255,145980.8230537888 2051823.3784043754,145954.36613142194 2051822.5028163139,145953.2894197142 2051852.3578379445))

Which corresponds to your second answer.

If I run:

    UPDATE spatial_ref_sys SET proj4text = '+proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +type=crs'
    WHERE srid=5070;

and check the results, as I would expect, I get the same answer:

POLYGON((145953.2894197142 2051852.3578379445,145979.74635016784 2051853.2334285255,145980.8230537888 2051823.3784043754,145954.36613142194 2051822.5028163139,145953.2894197142 2051852.3578379445))

This should only make a difference if you are updating your 3.1.5 instance which ships with PROJ 5.2.0 (predates the existence of the new proj api, so uses the spatial_ref_sys.proj4text)

If you put the 8.0.1 PROJ.4 string into the st_transform function on the 3.1.7 instance as so:

select
  st_astext( geom_5070 )
from st_geomfromtext(
  'POLYGON((
    -94.23933949466884 41.46010784655825,
    -94.23902022173033 41.46011128632304,
    -94.23901390824655 41.45984441537449,
    -94.23933317987697 41.45984097564519,
    -94.23933949466884 41.46010784655825))', 4326 ) geom_4326,
    st_transform( geom_4326, '+proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +type=crs' ) geom_5070;

You should get the old/matched results back:

POLYGON((145953.8164326792 2051852.6976375678,145980.27310072441 2051853.5729657637,145981.34988068073 2051823.7175862629,145954.89322071325 2051822.842260323,145953.8164326792 2051852.6976375678))

I was able to make this permanent by changing the auth_name column to something other than "EPSG" in the spatial_ref_sys table, i.e.:

UPDATE "public"."spatial_ref_sys" SET "auth_name" = 'JK', "proj4text" = '+proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +type=crs' WHERE "srid" = 5070;

[jbkoch]

So just to make sure there wasn't something going on with PROJ itself, I build the same dependency stack from source on a local machine:

PROJ v8.0.1
GEOS v3.9.1
GDAL v3.4.3

And I ran the following command (added the input and output files as attachments):

ogr2ogr -t_srs EPSG:5070 1884810.5070.ogr2ogr.gpkg 1884810.4326.gpkg

I got the same older matching results back:

POLYGON((145953.8164326792 2051852.6976375678,145980.27310072441 2051853.5729657637,145981.34988068073 2051823.7175862629,145954.89322071325 2051822.842260323,145953.8164326792 2051852.6976375678))