2

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]]
  • Transformation parameters in QGIS are usually to wgs84 EPSG:4326. Your second CRS is based on WGS84, so I wonder why you want to use EPSG:32634 as intermediate step. BTW your final CRS looks like EPSG:2462 Albanian 1987/GK zone 4, but all values with reversed signs. – AndreJ Jun 27 '17 at 16:11
  • Because i don't have parameters for going from 3127 to 2462, is there any way i can combine this, like arcgis does automatically ( composite transformation) – Shefqet LULJA Jun 27 '17 at 20:50
  • In QGIS, you always provide parameters to WGS84. Qgis does the inverse automatically. EPSG:2462 has the right parameters built-in, so you only have to define a custom CRS similar to EPSG:3127, but with your wgs84 parameters. Then QGIS can transform from custom to 2462 and backwards. – AndreJ Jun 28 '17 at 5:24
  • Official source of parameters: asig.gov.al/images/Relacioni%20permbledhes.pdf and asig.gov.al/images/DokumentaPDF/Udhezimi_4381.pdf (unless the server is down) – AndreJ Jun 28 '17 at 9:47
  • yes from there i got those – Shefqet LULJA Jun 28 '17 at 10:26
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

  • Qgis uses +towgs84=-44.183,-0.58,-38.489,-2.3867,-2.7072,3.5196,-8.2703 for EPSG:2462. EPSG:28404 uses a more general transformation set for the Pulkovo 1942 datum. – AndreJ Jun 27 '17 at 18:42
  • The signs of the parameters i posted are flipped because they are ready to be used directly in Arcgis which his matrice flip the signs, those parameters are not the original postes by ITRS – Shefqet LULJA Jun 27 '17 at 20:34
  • @ShefqetLULJA I'm sorry but I don't understand. When you gave the parameter values, did you give them in the direction matching the CRS ordering you gave? ArcGIS will take transformations in any order (doesn't have to be older-to-newer for instance) and supports both position vector and coordinate frame so you shouldn't have had to flip them for ArcGIS. – mkennedy Jun 27 '17 at 20:40
  • To be be more clear, if you put those as i posted you are ok in ArcGis, the problem is in Qgis now. – Shefqet LULJA Jun 27 '17 at 20:54
  • What transformation methods did you use in ArcGIS? – mkennedy Jun 27 '17 at 23:53
1

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.

  • For converting i can use arcgis , but i want to use the wmts service with this projection in qgis and on the fly projecting. The built in 2462 datum shift doesnt overlap very well – Shefqet LULJA Jun 28 '17 at 18:49
  • You are right. Can you try with the custom CRS from my updated answer? – AndreJ Jun 28 '17 at 19:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.