Reprojecting local CRS using PROJ

Question

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 (x₁,y₁).

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?

Answer 1

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.