Reprojecting local CRS using PROJ

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 '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. Feb 20 '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 '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? Feb 20 '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.

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:

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.