Project data between different datums using custom parameters in QGIS!

Question

I want to convert or reproject in the fly some shape files and raster between between different datums using custom transformation parameters (7 parameters) but don't know how to do that in QGIS.

I want to know the string of the final transformation, transformation matrice of rotation (my parameters are in seconds and scale is in ppm) that QGIS uses and if that string works only in one direction of the transformation or it automatically calculate the reverse transformation too?

In my scenario i want to go through this:

Epsg:3127 to Epsg:32634 (Transf. parameters are tx=0.0527 m ty=0.0509 m tz=-0.06636 m rx=-0.001456" ry=-0.008809" rz=0.014238" scale=0.000958 ppm)

Epsg:32634 to Epsg:28404 (Transf. parameters are tx=44.183 m ty=0.580 m tz=-38.489 m rx=2.3867" ry=2.7072" rz=-3.5196" scale=8.2703 ppm)

And how can i combine this 2 transformation paths in one so i can go directly from Epsg:3127 to Epsg:28404 and vice versa.

Here is the .gtf (32634 to 2642) contents created from the tool "create custom geographic transformation"

1DøÁÒ½  øu¼Î ]  GEOGTRAN["H_TRANS",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],GEOGCS["GCS_Pulkovo_1942",DATUM["D_Pulkovo_1942",SPHEROID["Krasovsky_1940",6378245.0,298.3]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],METHOD["Coordinate_Frame"],PARAMETER["X_Axis_Translation",44.183],PARAMETER["Y_Axis_Translation",0.58],PARAMETER["Z_Axis_Translation",38.489],PARAMETER["X_Axis_Rotation",2.3867],PARAMETER["Y_Axis_Rotation",2.7072],PARAMETER["Z_Axis_Rotation",-3.5196],PARAMETER["Scale_Difference",8.2703]]

and here is the .gtf (32634 to 3127 or 6870) contents created from the tool "create custom geographic transformation"

1DøÁÒ½  øu¼Î u  GEOGTRAN["UTM_KRGJSH_WGS_GK20FIN",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],GEOGCS["GCS_ETRS_1989",DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],METHOD["Coordinate_Frame"],PARAMETER["X_Axis_Translation",0.0527],PARAMETER["Y_Axis_Translation",0.0509],PARAMETER["Z_Axis_Translation",-0.06636],PARAMETER["X_Axis_Rotation",-0.001456],PARAMETER["Y_Axis_Rotation",-0.008809],PARAMETER["Z_Axis_Rotation",0.014238],PARAMETER["Scale_Difference",0.000958]]

Answer 1

If I understand how QGIS works, you would want to update the PROJ.4 strings in the 3127 and 28404 entries. Because, as AndreJ points out, you define the transformations within the ProjCRS definitions with the +wgs84 parameter, you don't need to pass through 32634 (WGS84 / UTM zone 34 North). Instead, the data will pass through 4326 (WGS84).

In 3127, add or update this parameter:

(updated based on ArcGIS gtf definition) +wgs84=-0.0527,-0.0509,0.06636,-0.001456,-0.008809,0.014238,-0.000958

This transformation converts from ETRS89 to WGS84. The method must be the same as EPSG's "position vector" method as that's the method that PROJ.4 uses. The ArcGIS gtf file is defined from WGS84 to ETRS89, so the translations and scale signed must be changed. The rotation signs need to be flipped for direction and then again to switch from the Coordinate Frame method to the Position Vector one, so they have the same signs as the gtf file.

For 28404, you again need to add or update the +wgs84 parameter to:

Note: Below updated based on the ArcGIS gtf file. Translation and scale signs were flipped to change the direction "to WGS84" because the ArcGIS version was "from WGS84". The rotation parameters signs remain the same because they needed to be flipped one for direction and again to change from the Coordinate Frame method to the Position Vector method.

+wgs84=-44.183,-0.580,-38.489,2.3867,2.7072,-3.5196,-8.2703

(original text based on question before the addition of the gtf file contents)

Note that I flipped all the signs to change the directionality to from Albanian 1987 to WGS84. HOWEVER, the signs of tz, rx, ry, and rz do not match the entry in the EPSG registry for 7834. The EPSG registry says that the transformation uses the coordinate frame method which means that the signs of the rotations should be changed to make it compatible with the position vector method that PROJ.4 uses. If I use the EPSG definition, here's the parameter:

+wgs84=-44.183,-0.580,-38.489,2.3867,2.7072,-3.5196,-8.2703

Answer 2

In QGIS, datum transformation parameters are defined from a geographic (latlong) CRS to WGS84 using the +towgs84parameter. The inverse datum shift is done internally, as well as transforming from one datum to another, and any projected CRS (like transverse mercator, UTM etc) based on the geographic CRS.

From the comments, it seems you want to transform from EPSG:6870 ETRS89 / Albania TM 2010 and EPSG:32634WGS 84 / UTM zone 34N to EPSG:2462 Albanian 1987 / Gauss-Kruger zone 4. WGS84 and ETRS89 are treated in QGIS as identical, with an accuracy of one meter.

Official parameters can be found in the EPSG database, and http://asig.gov.al/images/Relacioni%20permbledhes.pdf .

NOTE: the following calculations are done with QGIS 2.18.4 / GDAL 2.1.3 / PROJ 4.9.3. See Update 2 below for recent changes to the datum shift parameters

If you have QGIS running on Windows with proj.4 version 4.9.3, you have both included, and can run cs2cs on the command line. Putting the 6870 coordinates in a file named 6870.txt:

513384.34158 4471379.1036
513202.340 4471093.764

and running

cs2cs +init=epsg:6870 +to +init=epsg:2462 6870.txt >>out.txt

delivers

4428298.62  4471734.88 0.05
4428113.39  4471451.59 0.05

which is about 400m off from your Arcgis values.

Using the ETRS89 to WGS84 datum shift needs a full parameter set:

cs2cs +proj=tmerc +lat_0=0 +lon_0=20 +k=1 +x_0=500000 +y_0=0 +ellps=GRS80 +towgs84=0.0527,0.0509,-0.06636,-0.001456,-0.008809,0.014238,0.000958 +units=m +no_defs +to +init=epsg:2462 6870.txt >>out.txt

The result is:

4428299.07  4471735.02 0.07
4428113.84  4471451.73 0.07

So less than one meter away. Both points are near the village of Suke.

Note that epsg:2462 uses build-in datum shift values. You can get them with

gdalsrsinfo epsg:2462

PROJ.4 : '+proj=tmerc +lat_0=0 +lon_0=21 +k=1 +x_0=4500000 +y_0=0 +ellps=krass +towgs84=-44.183,-0.58,-38.489,-2.3867,-2.7072,3.5196,-8.2703 +units=m +no_defs '

---

UPDATE 1

With the datum shift string provided by mkennedy, the correct projection string for Albanian 1987 / Gauss-Kruger zone 4 would be

+proj=tmerc +lat_0=0 +lon_0=21 +k=1 +x_0=4500000 +y_0=0 +ellps=krass +towgs84=-44.183,-0.580,-38.489,2.3867,2.7072,-3.5196,-8.2703 +units=m +no_defs

With this, I get

4428592.16  4471838.25 0.40
4428406.93  4471554.97 0.40

and the points align in QGIS. Recalculating sample coordinates from http://asig.gov.al/images/Relacioni%20permbledhes.pdf works too, and it fits to the data in the EPSG database.

UPDATE 2

The issue was fixed in the EPSG database with changeset EPSG::2016.026 in September 2016, and in GDAL with Changeset 37081 in January 2017. It is now in the current OSGEO4W64 build of GDAL 2.2.0:

EPSG:2426 : '+proj=tmerc +lat_0=0 +lon_0=21 +k=1 +x_0=4500000 +y_0=0 +ellps=krass +towgs84=-44.183,-0.58,-38.489,2.3867,2.7072,-3.5196,-8.2703 +units=m +no_defs '

If you download the latest QGIS 2.18.10 for Windows, you have the correct datum shift values included. Unfortunately, Proj.4 has not yet been released after the commit, and Ubuntu packages don't have it either. So cs2cs will output wrong values, while on-the-fly-reprojection in the latest QGIS only on Windows will do it right.

You still have to add the ETRS89 datum shift manually, if you are unhappy with the shift of less than one meter against WGS84.