We have point data (x,y in m from a local 0,0) from an archaeological campaign, that we would like to display embedded into real-world cartography.

We know either the real world coordinates of 0,0 (e.g. 337861.49,5510278.69) and a rotation angle 48,5° or the coordinates of a second point (x1,y1).

Version A

One obvious idea would be to use an existing projection in PROJ like e.g. oblique mercator projection (like PROJ omerc) for this, however we end up with a large error (~200 m) and no rotation whatsoever.

Using PROJ tpqed we have a rotation but again a non-offsettable shift in E and N and also mismatched scaling.

Version B

Another option would be starting from a fitting cartesian, projected CRS like EPSG:25832 and just perform an affine transformation. This should be doable via SC_DerivedCRS in a WKT file, however I couldn't find any usable instructions how to do this with WKT2.

How do I wrangle the errors in A or the syntax for the WKT file?

  • Maybe working with the tools translategeometry and rotatefeatures in QGIS could suffice? – Erik Feb 20 at 13:01
  • I would like to have something that is using existing data and is not depending on arbitrary manual shifts. PROJ also has the added benefit of integrating this as a pipeline. – St. Hermes Feb 20 at 13:07
  • Yeah, but if you know one coordinate and a rotation angle, this should suffice. I doubt you're talking about areas several kilometers large? – Erik Feb 20 at 13:12
  • Your known origin point (337861.49,5510278.69) looks like it's in UTM, or State Plane (if this is in the USA)... is that right? If you know that point, you know the rotation angle between your local coordinate system and true north, and you know your local x,y points, and if they're in the same units (e.g., meters), why not just do the trigonometry to rotate each point yourself, and arrive at points in the same UTM zone or State Plane coordinate system? – Paulo Raposo Feb 20 at 14:22

Let me give you a Version C:

  • Find the Affine conversion parameters and process your coordinates with ogr2ogr and the -ct parameter to include a pipeline. Transform your original coordinates to a known CRS and set the WKT representation of that CRS to the transformed coordinates vector file.

The better way to go from a local survey Cartesian CRS to any geographic or projected CRS is through an intermediate geocentric CRS.

It is also a good option assume that your Cartesian CRS is a custom Transverse Mercator projection.

But since your CRS is not North oriented, in both cases you will need to rotate it with an Affine transformation pipeline.

So, lets divide the task. First, start with a North oriented CRS:

From https://upload.wikimedia.org/wikipedia/commons/2/2c/2D_affine_transformation_matrix.svg:


Therefore, for a teetha angle of 48.5 degrees, the pipeline to transform your data will be:

+proj=pipeline +step +proj=affine +s11=0.662620048 +s12=0.74895572 +s21=-0.74895572 +s22=0.662620048

With that pipeline you will be able to transform your original data to a North oriented CRS. You can define a custom Transverse Mercator CRS centered in your (0,0) point and set it for the transformed vector file.

Now, let me see the WGS84 geographic 3D coordinates:

C:\>cs2cs -f "%.8f" EPSG:25832 +to EPSG:4979
337861.49 5510278.69
49.72314746     6.75041204 0.00000000

They are latitude, longitude and ellipsoidal height assuming that the source coordinates were given in EPSG:25832 system and assuming the default datum transformation from ETRS89 and WGS84 (for this case it is without a datum transformation).

We will assume that the vertical of the place is coincident with the normal to the ellipsoid surface. But you will need to find the ellipsoidal height of your (0,0) point if you want to transform it with better accuracy. If you don't know the ellipsoidal height but you know that it is on the ground, take it from SRTM data.

Knowing the WGS84 latitude, longitude and ellipsoidal height of the center of your Cartesian North oriented (i.e., Topocentric) CRS, you can transform it from Topocentric to Geocentric coordinates:

From http://www.epsg.org/Portals/0/373-07-02.pdf:


Where lambda_0 is the longitude and phi_0 is the latitude of the topocentric point.

U, V ad W are easting, northing and elevation of the source coordinates to be converted.

To know the X_0, Y_0 and Z_0 parameters for the translation, convert the topocentric point geographic 3D coordinates to geocentric ones:

C:\>cs2cs -f "%.3f" EPSG:4979 +to EPSG:4978
49.72314746 6.75041204 0
4102764.717     485624.123 4842938.566

Same as before, calculate the rotation parameters and write the affine transformation pipeline. Transform your vector data and set to it EPSG:4978 (WGS84 Geocentric CRS). You will be able to project then from EPSG:4978 to your preferred CRS.

| improve this answer | |

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.