libspatialite7: erroneous geometry conversion using ST_Transform

Hello,

I have discovered a bug in libspatialite7’s ST_Transform() function where transforms from WGS84 geographic coordinates (EPSG:4326) to NAD83 UTM 10N (EPSG:26910) returns incorrect results on Ubuntu 20.04/Debian 11 (and perhaps other versions). The behavior on macOS and Alpine Linux appears to be correct. Note that transforms to WGS84 UTM 10N (EPSG:32610) seems to work fine on Ubuntu/Debian.

Here are the steps to reproduce:

# DB Setup

CREATE TABLE IF NOT EXISTS test(id INTEGER PRIMARY KEY);
SELECT AddGeometryColumn ('test', 'point', 4326, 'POINT', 'XY', 1);

INSERT INTO test (point) VALUES (GeomFromText('POINT(-122.338928 47.629273)', 4326));
INSERT INTO test (point) VALUES (GeomFromText('POINT(-122.338923 47.62927)', 4326));
INSERT INTO test (point) VALUES (GeomFromText('POINT(-122.33895 47.629253)', 4326));
INSERT INTO test (point) VALUES (GeomFromText('POINT(-122.392051 47.62888)', 4326));
INSERT INTO test (point) VALUES (GeomFromText('POINT(-122.392038 47.62888)', 4326));

# Results from Ubuntu 22.04 / Debian 11 (ARM64 and x86_64):

spatialite> SELECT AsText(point), AsText(ST_Transform(point, 26910)) FROM test;
POINT(-122.338928 47.629273)|POINT(549665.036792 5275308.276363)
POINT(-122.338923 47.62927)|POINT(549665.415273 5275307.94615)
POINT(-122.33895 47.629253)|POINT(549663.402919 5275306.039517)
POINT(-122.392051 47.62888)|POINT(545674.35684 5275231.954471)
POINT(-122.392038 47.62888)|POINT(545675.333511 5275231.962125)

spatialite> SELECT AsText(point), AsText(ST_Transform(point, 32610)) FROM test;
POINT(-122.338928 47.629273)|POINT(549665.158878 5275308.476981)
POINT(-122.338923 47.62927)|POINT(549665.53736 5275308.14677)
POINT(-122.33895 47.629253)|POINT(549663.525012 5275306.240133)
POINT(-122.392051 47.62888)|POINT(545674.479273 5275232.145082)
POINT(-122.392038 47.62888)|POINT(545675.455944 5275232.152738)
spatialite>

# Results from macOS 12 (Homebrew or MacPorts; ARM64 and x86_64):

spatialite> SELECT AsText(point), AsText(ST_Transform(point, 26910)) FROM test;
AsText(point) AsText(ST_Transform(point, 26910))
---------------------------- ----------------------------------
POINT(-122.338928 47.629273) POINT(549665.158879 5275308.47686)
POINT(-122.338923 47.62927) POINT(549665.537361 5275308.146648)
POINT(-122.33895 47.629253) POINT(549663.525012 5275306.240011)
POINT(-122.392051 47.62888) POINT(545674.479273 5275232.14496)
POINT(-122.392038 47.62888) POINT(545675.455944 5275232.152617)

spatialite> SELECT AsText(point), AsText(ST_Transform(point, 32610)) FROM test;
AsText(point) AsText(ST_Transform(point, 32610))
---------------------------- -----------------------------------
POINT(-122.338928 47.629273) POINT(549665.158878 5275308.476981)
POINT(-122.338923 47.62927) POINT(549665.53736 5275308.14677)
POINT(-122.33895 47.629253) POINT(549663.525012 5275306.240133)
POINT(-122.392051 47.62888) POINT(545674.479273 5275232.145082)
POINT(-122.392038 47.62888) POINT(545675.455944 5275232.152738)
spatialite>

# Results from Alpine Linux 3.15 (ARM64 and x86_64):

sqlite> SELECT AsText(point), AsText(ST_Transform(point, 26910)) FROM test;
POINT(-122.338928 47.629273)|POINT(549665.158879 5275308.47686)
POINT(-122.338923 47.62927)|POINT(549665.537361 5275308.146648)
POINT(-122.33895 47.629253)|POINT(549663.525012 5275306.240011)
POINT(-122.392051 47.62888)|POINT(545674.479273 5275232.14496)
POINT(-122.392038 47.62888)|POINT(545675.455944 5275232.152617)

sqlite> SELECT AsText(point), AsText(ST_Transform(point, 32610)) FROM test;
POINT(-122.338928 47.629273)|POINT(549665.158878 5275308.476981)
POINT(-122.338923 47.62927)|POINT(549665.53736 5275308.14677)
POINT(-122.33895 47.629253)|POINT(549663.525012 5275306.240133)
POINT(-122.392051 47.62888)|POINT(545674.479273 5275232.145082)
POINT(-122.392038 47.62888)|POINT(545675.455944 5275232.152738)

To make sure this wasn’t PROJ that was broken, I did the following PROJ queries:

# Results from Ubuntu 22.04 / Debian 11 (ARM64 and x86_64):

echo -122.338928 47.629273 | cct -c1,2 -d6 -z0 -t0 +proj=utm +zone=10 +ellps=GRS80 +datum=NAD83 +units=m +no_defs
549665.158878 5275308.476964 0.000057 0.0000

echo -122.338928 47.629273 | cct -c1,2 -d6 -z0 -t0 +proj=utm +zone=10 +ellps=WGS84 +datum=WGS84 +units=m +no_defs
549665.158878 5275308.476981 0.000000 0.0000

# Results from macOS 12 (Homebrew or MacPorts; ARM64 and x86_64):

echo -122.338928 47.629273 | cct -c1,2 -d6 -z0 -t0 +proj=utm +zone=10 +ellps=WGS84 +datum=WGS84 +units=m +no_defs
549665.158878 5275308.476981 0.000000 0.0000

echo -122.338928 47.629273 | cct -c1,2 -d6 -z0 -t0 +proj=utm +zone=10 +ellps=GRS80 +datum=NAD83 +units=m +no_defs
549665.158878 5275308.476964 0.000057 0.0000

# Results from Alpine Linux 3.15 (ARM64 and x86_64):

echo -122.338928 47.629273 | cct -c1,2 -d6 -z0 -t0 +proj=utm +zone=10 +ellps=GRS80 +datum=NAD83 +units=m +no_defs
549665.158878 5275308.476964 0.000057 0.0000

echo -122.338928 47.629273 | cct -c1,2 -d6 -z0 -t0 +proj=utm +zone=10 +ellps=WGS84 +datum=WGS84 +units=m +no_defs
549665.158878 5275308.476964 0.000057 0.0000

So it seems that the PROJ output is the same for all three (or six if you include architecture), so the problem likely lies in libspatialite7 on Debian/Ubuntu.

Please let me know if you need more information.

Thanks,

Brian