Bug#1031392: postgis: Polar stereographic axis order wrong when transforming coordinates
Paul Cochrane
paul.cochrane at posteo.de
Wed Feb 15 20:48:28 GMT 2023
Package: postgis
Version: 3.1.1+dfsg-1
Severity: important
Tags: patch upstream
Dear Maintainer,
at my employer, we primarily work with GIS data from the polar regions,
hence we are often converting coordinates between polar stereographic
projections (such as epsg:3976) and Mercator projection (epsg:4326). We
have also been very conservative when upgrading to new Debian versions,
hence some of our systems are still running "buster". When upgrading
one system to "bullseye", we found that the test suite for our software
failed, giving incorrect output for coordinates transformed from
epsg:3976 to epsg:4326.
The steps to reproduce this issue are:
- create a test "bullseye" system (e.g. in vagrant)
- ensure the system is up to date:
`sudo apt update && sudo apt upgrade`
- install Postgres with the PostGIS extensions:
`sudo apt install -y postgresql postgis`
- create a test database:
`sudo -u postgres createdb -E UTF8 -T template0 test_db`
- create the PostGIS extension on the test database:
`sudo -u postgres psql -c "CREATE EXTENSION IF NOT EXISTS postgis;" test_db`
- start `psql` and transform a test point from epsg:3976 to epsg:4326
```
sudo -u postgres psql -d test_db
psql (13.9 (Debian 13.9-0+deb11u1))
Type "help" for help.
test_db=# select ST_AsEWKT(ST_Transform(ST_GeomFromEWKT('SRID=3976;POINT(-25000 75000)'), 4326));
st_asewkt
--------------------------------------------------------
SRID=4326;POINT(108.43494882292202 -89.27021258118658)
(1 row)
```
The expected output is:
```
test_db=# select ST_AsEWKT(ST_Transform(ST_GeomFromEWKT('SRID=3976;POINT(-25000 75000)'), 4326));
st_asewkt
--------------------------------------------------------
SRID=4326;POINT(-18.43494882292201 -89.27021258118658)
(1 row)
```
As one can see, the first element of the coordinate pair is completely
wrong.
We managed to isolate the problem to the issue [Polar Stereographic Axis
order wrong](https://trac.osgeo.org/postgis/ticket/4842) in the upstream
PostGIS project. The issue has thus already been fixed upstream as of
PostGIS version 3.1.2. We also checked the behaviour on "bookworm"
(which uses PostGIS 3.3.2) and could confirm that the issue does not
appear there.
We were also able to patch the Debian source package
(postgis-3.1.1+dfsg) with the changeset mentioned in the PostGIS ticket
mentioned above, build a new package, install it, and can confirm that
the correct behaviour is observed with the patch included.
To create the patch, we downloaded the source package with `apt source
postgis`, making a copy of the source package tree to save a reference
copy of the source package:
```
cp -a postgis-3.1.1+dfsg/ postgis-3.1.1+dfsg.orig/
```
We then downloaded [the diff of the changeset fixing the
bug](https://trac.osgeo.org/postgis/changeset/4f6d81f8843914666b00683dcc93b81f6c13635d/git?format=diff&new=4f6d81f8843914666b00683dcc93b81f6c13635d)
calling the file fix-axis-flips-for-lat-first-crs-postgis-orig.patch and
applied it to the source package:
```
cd postgis-3.1.1+dfsg/
patch -p0 < ~/fix-axis-flips-for-lat-first-crs-postgis-orig.patch
```
The patch applied successfully, however with offsets, hence the file
./liblwgeom/lwgeom_transform.c.orig was created in the source package
tree by the `patch` command. We removed this file and then diffed the
original and patched source packages to create a patch which applies
cleanly to the original Debian package:
```
diff -Nru postgis-3.1.1+dfsg.orig/ postgis-3.1.1+dfsg/ > fix-axis-flips-for-lat-first-crs.patch
```
This patch is attached to this bug report.
Rebuilding the source package with this patch (debuild -b -uc -us) and
installing the patched package fixed the issue and we receive correct
output when converting coordinates from polar stereographic projections
to Mercator.
I'm aware that this issue probably only affects a small number of users,
however for my employer this is a show-stopper because we would be
delivering incorrect information to customers, therefore we can't
upgrade our "buster" systems to "bullseye" at this stage.
I realise that one option for my employer is to simply wait until
"bookworm" is released and to skip upgrading to "bullseye" completely.
However, it's possible that this issue affects other PostGIS users who
are running "bullseye" systems and who might not be aware that they have
the problem, thus we decided it worthwhile to report the issue.
If you need further information from me, please don't hesitate to ask!
I hope this has helped in some way.
Kind regards,
Paul Cochrane
-- System Information:
Debian Release: 11.6
APT prefers stable-updates
APT policy: (500, 'stable-updates'), (500, 'stable-security'), (500, 'stable')
Architecture: amd64 (x86_64)
Kernel: Linux 5.10.0-21-amd64 (SMP w/2 CPU threads)
Locale: LANG=C.UTF-8, LC_CTYPE=C.UTF-8 (charmap=UTF-8), LANGUAGE not set
Shell: /bin/sh linked to /usr/bin/dash
Init: systemd (via /run/systemd/system)
LSM: AppArmor: enabled
Versions of packages postgis depends on:
ii libc6 2.31-13+deb11u5
ii libgdal28 3.2.2+dfsg-2+deb11u2
ii libgeos-c1v5 3.9.0-1
ii libpq5 13.9-0+deb11u1
ii libproj19 7.2.1-1
Versions of packages postgis recommends:
ii postgis-doc 3.1.1+dfsg-1
ii postgresql-13-postgis-3 [postgresql-postgis] 3.1.1+dfsg-1
postgis suggests no packages.
-- no debconf information
-------------- next part --------------
diff -Nru postgis-3.1.1+dfsg.orig/liblwgeom/lwgeom_transform.c postgis-3.1.1+dfsg/liblwgeom/lwgeom_transform.c
--- postgis-3.1.1+dfsg.orig/liblwgeom/lwgeom_transform.c 2021-01-28 19:10:32.000000000 +0000
+++ postgis-3.1.1+dfsg/liblwgeom/lwgeom_transform.c 2023-02-15 14:26:05.392913740 +0000
@@ -245,7 +245,7 @@
const char *out_name, *out_abbrev, *out_direction;
double out_unit_conv_factor;
const char *out_unit_name, *out_unit_auth_name, *out_unit_code;
- /* Read only first axis, see if it is degrees / north */
+ /* Read only first axis */
proj_cs_get_axis_info(NULL,
pj_cs,
0,
@@ -256,13 +256,10 @@
&out_unit_name,
&out_unit_auth_name,
&out_unit_code);
- if (strcasecmp(out_abbrev, "E") == 0)
- rv = LW_FALSE;
- else if (strcasecmp(out_abbrev, "Lat") == 0)
- rv = LW_TRUE;
- else if (strcasecmp(out_direction, "north") == 0)
- rv = LW_TRUE;
- else if (strcasecmp(out_direction, "south") == 0)
+
+ /* Only swap Lat/Lon systems */
+ /* Use whatever ordering planar systems default to */
+ if (strcasecmp(out_abbrev, "Lat") == 0)
rv = LW_TRUE;
else
rv = LW_FALSE;
diff -Nru postgis-3.1.1+dfsg.orig/regress/core/tickets.sql postgis-3.1.1+dfsg/regress/core/tickets.sql
--- postgis-3.1.1+dfsg.orig/regress/core/tickets.sql 2021-01-28 19:10:32.000000000 +0000
+++ postgis-3.1.1+dfsg/regress/core/tickets.sql 2023-02-15 14:26:05.408913778 +0000
@@ -1298,7 +1298,9 @@
SELECT '#4689', _ST_DistanceTree('POLYGON ((30 10, 40 40, 20 40, 30 10))'::geography, 'POLYGON((81 6,140 35,-70 18,-51 0,-60 -46,106 -6,81 6))');
+-- Polar Stereographic Axis order issues
SELECT '#4748', ST_AsEWKT(ST_Transform(ST_SetSRID(ST_Point(-36.75, -54.25), 4326), 3031),1);
+SELECT '#4842', ST_AsEWKT(ST_Transform(ST_SetSRID(ST_Point(12.572, 66.081), 4326), 3413),1);
SELECT '#4718',
round(degrees(
diff -Nru postgis-3.1.1+dfsg.orig/regress/core/tickets_expected postgis-3.1.1+dfsg/regress/core/tickets_expected
--- postgis-3.1.1+dfsg.orig/regress/core/tickets_expected 2021-01-28 19:10:32.000000000 +0000
+++ postgis-3.1.1+dfsg/regress/core/tickets_expected 2023-02-15 14:26:05.424913817 +0000
@@ -435,6 +435,7 @@
ERROR: LWGEOM_addpoint: Invalid offset
#4689|0
#4748|SRID=3031;POINT(-2399498.7 3213318.5)
+#4842|SRID=3413;POINT(2218086.2 -1409161.3)
#4718|270.000|90.000
#4727|0
#4796|POLYGON((178632 397744,178738 397769,178745 397741,178749 397742,178752 397727,178748 397726,178759 397684,178699 397669,178694 397691,178646 397680,178632 397744))
More information about the Pkg-grass-devel
mailing list