Finding applied transformation method and accuracy information when reprojecting layer in QGIS
Geographic Information Systems Asked by Aron on January 16, 2021
When I load two layers with different CRS I’m prompted to choose a transformation method and there is a column indicating the accuracy (see image 1). I would like to know where the information as part of the on-the-fly (OTF) transformation is coming from.
As far as I can tell, there is no such information displayed when I try to reproject a layer into a new CRS which is based on a different ellipsoid, using the ‘Reproject Layer’ tool. All I seem to be able to do is to choose the target CRS but there is no option to choose a particular method, neither there is the info about which method is applied by QGIS, and what the accuracy implications are (see image 2).
How can I address that particular component shown as part of the OTF with Python?
My research indicates I have to make use of the proj4 library but I don’t know how to do it yet.
Is the transformation code (tfm code – see image 1) also part of the library and how to address it?
Background of my question: As part of my coursework I intend to write a plugin which lists possible transformation methods and the accuracy when using the Reproject tool, similar to the OTF transformation. If this isn’t possible (I can imagine there is a reason why it is not already there), my alternative goal is to display the information about which method is applied and the accuracy at least.
I have carefully read the related posts on this plattform, however I get the impression they don’t help with what I’m after.
I am using QGIS version 3.10.5.
3 Answers
Proj library has a copy of the contents of the EPSG database in the SQLite database "proj.db". That's the immediate source that QGIS is using. The ultimate source is EPSG database, for example
http://www.epsg-registry.org/export.htm?gml=urn:ogc:def:coordinateOperation:EPSG::15948
reports the accuracy in section
<gml:coordinateOperationAccuracy>
<gmd:DQ_RelativeInternalPositionalAccuracy xmlns:gmd="http://www.isotc211.org/2005/gmd">
<gmd:result>
<gmd:DQ_QuantitativeResult>
<gmd:valueUnit xlink:href="urn:ogc:def:uom:EPSG::9001" xmlns:xlink="http://www.w3.org/1999/xlink"/>
<gmd:value>
<gco:Record xmlns:gco="http://www.isotc211.org/2005/gco">
<gco:Decimal>0.9</gco:Decimal>
</gco:Record>
</gmd:value>
</gmd:DQ_QuantitativeResult>
</gmd:result>
</gmd:DQ_RelativeInternalPositionalAccuracy>
</gml:coordinateOperationAccuracy>
From the proj.db you can find the info with SQL query
select * from grid_transformation
where code=15948;
Answered by user30184 on January 16, 2021
If you upgrade to QGIS 3.12, then you'll see the transformation details directly in the reproject algorithms and be able to manually pick a transform from there.
Answered by ndawson on January 16, 2021
I don't think you can do it with the proj4 library. I don't know if you can do it with pyproj.
What you must do is run projinfo command, that once installed is accessible from the command console. You will see how to implement it in Python.
Also take this answer just as a hint. The transformation is also needed in the export tool, and must be checked if there is not already defined a default datum transformation in the settings or the project properties.
From my environment:
C:>projinfo
Rel. 6.3.1, February 10th, 2020
usage: projinfo [-o formats] [-k crs|operation|ellipsoid] [--summary] [-q]
([--area name_or_code] | [--bbox west_long,south_lat,east_long,north_lat])
[--spatial-test contains|intersects]
[--crs-extent-use none|both|intersection|smallest]
[--grid-check none|discard_missing|sort] [--show-superseded]
[--pivot-crs always|if_no_direct_transformation|never|{auth:code[,auth:code]*}]
[--boundcrs-to-wgs84]
[--main-db-path path] [--aux-db-path path]*
[--identify] [--3d]
[--c-ify] [--single-line]
{object_definition} | (-s {srs_def} -t {srs_def})
-o: formats is a comma separated combination of: all,default,PROJ,WKT_ALL,WKT2:2015,WKT2:2019,WKT1:GDAL,WKT1:ESRI,PROJJSON
Except 'all' and 'default', other format can be preceded by '-' to disable them
{object_definition} might be a PROJ string, a WKT string, a AUTHORITY:CODE, or urn:ogc:def:OBJECT_TYPE:AUTHORITY::CODE
Get a summary of the available transformations from EPSG:3068 to EPSG:25833:
C:>projinfo --summary -s EPSG:3068 -t EPSG:25833
Candidate operations found: 2
unknown id, Inverse of Soldner Berlin + DHDN to ETRS89 (2) + UTM zone 33N, 3 m, Germany - West Germany all states
unknown id, Inverse of Soldner Berlin + DHDN to ETRS89 (8) + UTM zone 33N, 0.9 m, Germany - onshore, at least one grid missing
Get the PROJ pipelines format of the available transformations:
C:>projinfo -o PROJ -s EPSG:3068 -t EPSG:25833
Candidate operations found: 2
-------------------------------------
Operation No. 1:
unknown id, Inverse of Soldner Berlin + DHDN to ETRS89 (2) + UTM zone 33N, 3 m, Germany - West Germany all states
PROJ string:
+proj=pipeline +step +proj=axisswap +order=2,1 +step +inv +proj=cass +lat_0=52.4186482777778 +lon_0=13.6272036666667 +x_0=40000 +y_0=10000 +ellps=bessel +step +proj=push +v_3 +step +proj=cart +ellps=bessel +step +proj=helmert +x=598.1 +y=73.7 +z=418.2 +rx=0.202 +ry=0.045 +rz=-2.455 +s=6.7 +convention=position_vector +step +inv +proj=cart +ellps=GRS80 +step +proj=pop +v_3 +step +proj=utm +zone=33 +ellps=GRS80
-------------------------------------
Operation No. 2:
unknown id, Inverse of Soldner Berlin + DHDN to ETRS89 (8) + UTM zone 33N, 0.9 m, Germany - onshore, at least one grid missing
PROJ string:
+proj=pipeline +step +proj=axisswap +order=2,1 +step +inv +proj=cass +lat_0=52.4186482777778 +lon_0=13.6272036666667 +x_0=40000 +y_0=10000 +ellps=bessel +step +proj=hgridshift +grids=BETA2007.gsb +step +proj=utm +zone=33 +ellps=GRS80
Grid BETA2007.gsb needed but not found on the system. Can be obtained from the proj-datumgrid package at https://download.osgeo.org/proj/proj-datumgrid-1.8.zip
Get the default information (PROJ pipelines and WKT2:2019 definitions) of the available transformations:
C:>projinfo -s EPSG:3068 -t EPSG:25833
Candidate operations found: 2
-------------------------------------
Operation No. 1:
unknown id, Inverse of Soldner Berlin + DHDN to ETRS89 (2) + UTM zone 33N, 3 m, Germany - West Germany all states
PROJ string:
+proj=pipeline +step +proj=axisswap +order=2,1 +step +inv +proj=cass +lat_0=52.4186482777778 +lon_0=13.6272036666667 +x_0=40000 +y_0=10000 +ellps=bessel +step +proj=push +v_3 +step +proj=cart +ellps=bessel +step +proj=helmert +x=598.1 +y=73.7 +z=418.2 +rx=0.202 +ry=0.045 +rz=-2.455 +s=6.7 +convention=position_vector +step +inv +proj=cart +ellps=GRS80 +step +proj=pop +v_3 +step +proj=utm +zone=33 +ellps=GRS80
WKT2:2019 string:
CONCATENATEDOPERATION["Inverse of Soldner Berlin + DHDN to ETRS89 (2) + UTM zone 33N",
SOURCECRS[
PROJCRS["DHDN / Soldner Berlin",
BASEGEOGCRS["DHDN",
DATUM["Deutsches Hauptdreiecksnetz",
ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4314]],
CONVERSION["Soldner Berlin",
METHOD["Cassini-Soldner",
ID["EPSG",9806]],
PARAMETER["Latitude of natural origin",52.4186482777778,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",13.6272036666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["False easting",40000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",10000,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (x)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (y)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
ID["EPSG",3068]]],
TARGETCRS[
PROJCRS["ETRS89 / UTM zone 33N",
BASEGEOGCRS["ETRS89",
DATUM["European Terrestrial Reference System 1989",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4258]],
CONVERSION["UTM zone 33N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",15,
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]],
ID["EPSG",25833]]],
STEP[
CONVERSION["Inverse of Soldner Berlin",
METHOD["Inverse of Cassini-Soldner",
ID["INVERSE(EPSG)",9806]],
PARAMETER["Latitude of natural origin",52.4186482777778,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",13.6272036666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["False easting",40000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",10000,
LENGTHUNIT["metre",1],
ID["EPSG",8807]],
ID["INVERSE(EPSG)",19996]]],
STEP[
COORDINATEOPERATION["DHDN to ETRS89 (2)",
VERSION["IfAG-Deu W 3m"],
SOURCECRS[
GEOGCRS["DHDN",
DATUM["Deutsches Hauptdreiecksnetz",
ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4314]]],
TARGETCRS[
GEOGCRS["ETRS89",
DATUM["European Terrestrial Reference System 1989",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4258]]],
METHOD["Position Vector transformation (geog2D domain)",
ID["EPSG",9606]],
PARAMETER["X-axis translation",598.1,
LENGTHUNIT["metre",1],
ID["EPSG",8605]],
PARAMETER["Y-axis translation",73.7,
LENGTHUNIT["metre",1],
ID["EPSG",8606]],
PARAMETER["Z-axis translation",418.2,
LENGTHUNIT["metre",1],
ID["EPSG",8607]],
PARAMETER["X-axis rotation",0.202,
ANGLEUNIT["arc-second",4.84813681109536E-06],
ID["EPSG",8608]],
PARAMETER["Y-axis rotation",0.045,
ANGLEUNIT["arc-second",4.84813681109536E-06],
ID["EPSG",8609]],
PARAMETER["Z-axis rotation",-2.455,
ANGLEUNIT["arc-second",4.84813681109536E-06],
ID["EPSG",8610]],
PARAMETER["Scale difference",6.7,
SCALEUNIT["parts per million",1E-06],
ID["EPSG",8611]],
OPERATIONACCURACY[3.0],
ID["EPSG",1776],
REMARK["Mean of 109 stations. Replaces DHDN to ETRS89 (1) (tfm code 1309). May be taken as approximate transformation DHDN to WGS 84 - see code 1777."]]],
STEP[
CONVERSION["UTM zone 33N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",15,
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]],
ID["EPSG",16033]]],
USAGE[
SCOPE["unknown"],
AREA["Germany - West Germany all states"],
BBOX[47.27,5.87,55.09,13.84]]]
-------------------------------------
Operation No. 2:
unknown id, Inverse of Soldner Berlin + DHDN to ETRS89 (8) + UTM zone 33N, 0.9 m, Germany - onshore, at least one grid missing
PROJ string:
+proj=pipeline +step +proj=axisswap +order=2,1 +step +inv +proj=cass +lat_0=52.4186482777778 +lon_0=13.6272036666667 +x_0=40000 +y_0=10000 +ellps=bessel +step +proj=hgridshift +grids=BETA2007.gsb +step +proj=utm +zone=33 +ellps=GRS80
WKT2:2019 string:
CONCATENATEDOPERATION["Inverse of Soldner Berlin + DHDN to ETRS89 (8) + UTM zone 33N",
SOURCECRS[
PROJCRS["DHDN / Soldner Berlin",
BASEGEOGCRS["DHDN",
DATUM["Deutsches Hauptdreiecksnetz",
ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4314]],
CONVERSION["Soldner Berlin",
METHOD["Cassini-Soldner",
ID["EPSG",9806]],
PARAMETER["Latitude of natural origin",52.4186482777778,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",13.6272036666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["False easting",40000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",10000,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (x)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (y)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
ID["EPSG",3068]]],
TARGETCRS[
PROJCRS["ETRS89 / UTM zone 33N",
BASEGEOGCRS["ETRS89",
DATUM["European Terrestrial Reference System 1989",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4258]],
CONVERSION["UTM zone 33N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",15,
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]],
ID["EPSG",25833]]],
STEP[
CONVERSION["Inverse of Soldner Berlin",
METHOD["Inverse of Cassini-Soldner",
ID["INVERSE(EPSG)",9806]],
PARAMETER["Latitude of natural origin",52.4186482777778,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",13.6272036666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["False easting",40000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",10000,
LENGTHUNIT["metre",1],
ID["EPSG",8807]],
ID["INVERSE(EPSG)",19996]]],
STEP[
COORDINATEOPERATION["DHDN to ETRS89 (8)",
VERSION["BKG-Deu BeTA2007"],
SOURCECRS[
GEOGCRS["DHDN",
DATUM["Deutsches Hauptdreiecksnetz",
ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4314]]],
TARGETCRS[
GEOGCRS["ETRS89",
DATUM["European Terrestrial Reference System 1989",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4258]]],
METHOD["NTv2",
ID["EPSG",9615]],
PARAMETERFILE["Latitude and longitude difference file","BETA2007.gsb"],
OPERATIONACCURACY[0.9],
ID["EPSG",15948],
REMARK["Developed for ATKIS (Amtliches Topographisch-Kartographisches Informationssystem [Official Topographic and Cartographic Information System]). Provides a uniform transformation across the whole country. May be used as tfm to WGS84 - see tfm code 15949."]]],
STEP[
CONVERSION["UTM zone 33N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",15,
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]],
ID["EPSG",16033]]],
USAGE[
SCOPE["unknown"],
AREA["Germany - onshore"],
BBOX[47.27,5.86,55.09,15.04]]]
Grid BETA2007.gsb needed but not found on the system. Can be obtained from the proj-datumgrid package at https://download.osgeo.org/proj/proj-datumgrid-1.8.zip
That is intended to know the information, then you must apply the selected transformation with GDAL (for raster formats) or OGR (for vector formats) in the reprojection process.
Answered by Gabriel De Luca on January 16, 2021
Add your own answers!
Ask a Question
Get help from others!
Recent Answers
- Peter Machado on Why fry rice before boiling?
- haakon.io on Why fry rice before boiling?
- Jon Church on Why fry rice before boiling?
- Joshua Engel on Why fry rice before boiling?
- Lex on Does Google Analytics track 404 page responses as valid page views?
Recent Questions
- How can I transform graph image into a tikzpicture LaTeX code?
- How Do I Get The Ifruit App Off Of Gta 5 / Grand Theft Auto 5
- Iv’e designed a space elevator using a series of lasers. do you know anybody i could submit the designs too that could manufacture the concept and put it to use
- Need help finding a book. Female OP protagonist, magic
- Why is the WWF pending games (“Your turn”) area replaced w/ a column of “Bonus & Reward”gift boxes?