Discrepancies between ST_Transform() and PROJ on Debian version of libspatialite7

[Brian Miles]

Hello,

I have discovered discrepancies in results from the ST_Transform() function where transforms from WGS84 geographic coordinates (EPSG:4326) to NAD83 UTM (EPSG:26910 specifically) returns different/unexpected/incorrect results (as in differing from output from PROJ directly) on Debian 11 (using libspatialite7:5.0.1-2).

The behavior on macOS and Alpine Linux appears to be the same as each other, but different from Debian. Further, macOS and Alpine results are consistent with PROJ output. Note that transforms to WGS84 UTM (EPSG:32610) seems to work fine on 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 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)


# Commentary
Interestingly, when I run spatialite with PROJ_NETWORK=ON on both macOS and Alpine, I get the same results as under Debian (both with and without PROJ_NETWORK=ON). This is the case even though project resource files are installed on all three (the proj-data file and bundled with the Homebrew proj package, as far as I can tell: https://github.com/Homebrew/homebrew-core/blob/HEAD/Formula/proj.rb#L32).

# Comparison to PROJ
In either case, the results from spatialite on Debian differ from what I am seeing from running PROJ directly (which matches the default output from spatialite in macOS and Alpine):

## PROJ output from Debian 11 (ARM64 and x86_64):

$ echo -122.338928 47.629273 | proj -d 6 +proj=utm +zone=10 +ellps=GRS80 +datum=NAD83 +units=m +no_defs
549665.158879    5275308.476860

$ echo -122.338928 47.629273 | proj -d 6 +proj=utm +zone=10 +ellps=WGS84 +datum=WGS84 +units=m +no_defs
549665.158878    5275308.476981


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

$ echo -122.338928 47.629273 | proj -d 6 +proj=utm +zone=10 +ellps=GRS80 +datum=NAD83 +units=m +no_defs
549665.158879    5275308.476860

$ echo -122.338928 47.629273 | proj -d 6 +proj=utm +zone=10 +ellps=WGS84 +datum=WGS84 +units=m +no_defs
549665.158878    5275308.476981


## PROJ output from Alpine Linux 3.15 (ARM64 and x86_64):

$ echo -122.338928 47.629273 | proj -d 6 +proj=utm +zone=10 +ellps=GRS80 +datum=NAD83 +units=m +no_defs
549665.158879   5275308.476860

$ echo -122.338928 47.629273 | proj -d 6 +proj=utm +zone=10 +ellps=WGS84 +datum=WGS84 +units=m +no_defs
549665.158878   5275308.476981


# Summary
So it seems that the PROJ output is the same for all three (or six if you include architecture), and there may be some problem in libspatialite7 on Debian, but I am not sure how to move forward with this. I have reported this to the Debian libspatialite7 maintainers, and they suggested I report this upstream, hence this e-mail.


Please let me know if you need more information.

Thanks,

Brian

[Even Rouault]

Maybe a partial response, but hopefully relevant points:

- to make a fair comparison, you need to ensure that all PROJ versions are the same. There are continuous enhancement/bug fixes (or sometimes regressions) that affect results if using different PROJ versions

- the homebrew recipee https://github.com/Homebrew/homebrew-core/blob/HEAD/Formula/proj.rb#L32 points to the old (and small) proj-datumgrid-1.8 package. When using PROJ_NETWORK=ON, you indirectly use the full big package of grids that can be locally downloaded from http://download.osgeo.org/proj/proj-data-1.9.tar.gz . And NAD83 to WGS84 transformations are impacted because proj-data-1.9 includes grids for such datum transformation, whereas proj-datumgrid-1.8 does not

Even

--
You received this message because you are subscribed to the Google Groups "SpatiaLite Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to spatialite-use...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/spatialite-users/81901984-df95-4d7d-a984-bcef420a0fe1n%40googlegroups.com.

-- 
http://www.spatialys.com
My software is free, but my time generally not.

[Brian Miles]

Hi Even,

This appears to be part of the issue, at least from libspatialite’s point of view. I tried installing Alpine:edge, which provides PROJ 9 and proj-data now uses https://download.osgeo.org/proj/proj-data-1.9.zip. The results now correspond to what I am seeing in Debian 11.

I am still confused as to why the `proj` binary is returning different results from SpatiaLite.

Thanks so much for your help.

Brian