Bug#1008612: libspatialite7: erroneous geometry conversion using ST_Transform
Sebastiaan Couwenberg
sebastic at xs4all.nl
Wed Mar 30 18:07:34 BST 2022
On 3/30/22 18:44, Brian Miles wrote:
> I am a little confused as to why we are seeing different results between spatialite and proj itself for what should be the same transformation.
SpatiaLite has its own CRS definitions in the spatial_ref_sys table like
PostGIS:
$ spatialite -line /tmp/test.db "SELECT * FROM spatial_ref_sys WHERE
srid IN (26910, 32610)"
SpatiaLite version ..: 5.0.1 Supported Extensions:
- 'VirtualShape' [direct Shapefile access]
- 'VirtualDbf' [direct DBF access]
- 'VirtualText' [direct CSV/TXT access]
- 'VirtualGeoJSON' [direct GeoJSON access]
- 'VirtualXL' [direct XLS access]
- 'VirtualNetwork' [Dijkstra shortest path - obsolete]
- 'RTree' [Spatial Index - R*Tree]
- 'MbrCache' [Spatial Index - MBR cache]
- 'VirtualFDO' [FDO-OGR interoperability]
- 'VirtualBBox' [BoundingBox tables]
- 'VirtualSpatialIndex' [R*Tree metahandler]
- 'VirtualElementary' [ElemGeoms metahandler]
- 'VirtualRouting' [Dijkstra shortest path - advanced]
- 'VirtualKNN' [K-Nearest Neighbors metahandler]
- 'VirtualGPKG' [OGC GeoPackage interoperability]
- 'VirtualXPath' [XML Path Language - XPath]
- 'SpatiaLite' [Spatial SQL - OGC]
PROJ version ........: Rel. 9.0.0, March 1st, 2022
GEOS version ........: 3.10.2-CAPI-1.16.0
RTTOPO version ......: 1.1.0
TARGET CPU ..........: x86_64-linux-gnu
srid = 26910
auth_name = epsg
auth_srid = 26910
ref_sys_name = NAD83 / UTM zone 10N
proj4text = +proj=utm +zone=10 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0
+units=m +no_defs
srtext = PROJCS["NAD83 / UTM zone
10N",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS
1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","26910"]]
srid = 32610
auth_name = epsg
auth_srid = 32610
ref_sys_name = WGS 84 / UTM zone 10N
proj4text = +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
srtext = PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS
84",DATUM["WGS_1984",SPHEROID["WGS
84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]
These differ from the definitions in PROJ:
$ projinfo -k crs EPSG:26910
PROJ.4 string:
+proj=utm +zone=10 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
+type=crs
WKT2:2019 string:
PROJCRS["NAD83 / UTM zone 10N",
BASEGEOGCRS["NAD83",
DATUM["North American Datum 1983",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4269]],
CONVERSION["UTM zone 10N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",-123,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Engineering survey, topographic mapping."],
AREA["North America - between 126°W and 120°W - onshore and
offshore. Canada - British Columbia; Northwest Territories; Yukon.
United States (USA) - California; Oregon; Washington."],
BBOX[30.54,-126,81.8,-119.99]],
ID["EPSG",26910]]
$ projinfo -k crs EPSG:32610
PROJ.4 string:
+proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +type=crs
WKT2:2019 string:
PROJCRS["WGS 84 / UTM zone 10N",
BASEGEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["UTM zone 10N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",-123,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Engineering survey, topographic mapping."],
AREA["Between 126°W and 120°W, northern hemisphere between
equator and 84°N, onshore and offshore. Canada - British Columbia (BC);
Northwest Territories (NWT); Nunavut; Yukon. United States (USA) -
Alaska (AK)."],
BBOX[0,-126,84,-120]],
ID["EPSG",32610]]
> For example:
>
> 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
>
> this matches the value you are getting below when PROJ_NETWORK is not set.
>
> Incidentally, I believe I do have proj-data installed:
>
> ```
> $ dpkg -s proj-data
> Package: proj-data
> [...]
> Source: proj
> Version: 7.2.1-1
> Description: Cartographic projection filter and library (datum package)
> [...]
> .
> This package contains auxiliary projection datum grids used by the
> library and tools.
> Homepage: https://proj.org/
> ```
>
> And I received the same result when running `cct` with `PROJ_NETWORK=ON`.
>
> Any ideas what’s going on here?
That's a different proj-data, it's the binary package built from the
proj source package containing datumgrids from the proj-datumgrid
upstream project (amongst others). Those grids are no longer updated,
all updates happen in the PROJ-data upstream project now:
https://proj.org/resource_files.html#proj-data
Due to the unresolved license issue the proj-data source package
building the proj-data-grids binary packages aren't in Debian, and since
no progress is being made getting the issue resolved it's not getting
into Debian any time soon.
If you need those grids, getting them from the CDN by using
PROJ_NETWORK=ON is the easiest option.
Kind Regards,
Bas
--
GPG Key ID: 4096R/6750F10AE88D4AF1
Fingerprint: 8182 DE41 7056 408D 6146 50D1 6750 F10A E88D 4AF1
More information about the Pkg-grass-devel
mailing list