Custom Projection in Proj4 has errors, what can I do to implement a 4 parameter transform?

Question

I'm working on a project where I'd like to define a projection in QGIS, using a Proj4 string, so that a user can import and export site coordinates from the local UTM or geographic coordinates.

I've had a good search around the place, and the best resource I found is this gis.stackexchange post, which gets me close, but I have errors still, which you can see in the images below. They are not systematic in a way that I can see, but might be related to scale.

This is my Proj4 string:

+proj=omerc +lat_0=-42.83254553 +lonc=147.31794080 +alpha=43.297962 +k_0=0.999592 +x_0=0 +y_0=0 +gamma=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

And I have used a bit of software to generate transformation parameters from known coordinate relationships, which results in these:

• X shift: 525985.101348
• Y shift: 5257731.440819
• Z rotation: 0.755692
• Scale: 0.999592

These were what I used to get the parameters in the Proj4 string, including coming up with the origin Lon and Lat (via transforming the 0,0 point back to MGA and thence forward to WGS84).

If I use the above parameters to do a manual transformation, the coordinates are transformed perfectly. If I use the above projection in QGIS I'm still getting errors of up to 2 m.

I wonder if anyone knows of a better projection string to use in Proj4, or if there is another way that can handle a site datum or custom coordinate system with a transform into a known system.

EDIT: more info

Here's some coordinate pairs. Note that the lat/lon are derived from the MGA coordinates. E and N mga is GDA94, MGA Zone 55. The Zs are irrelevant it's a simple vertical shift.

Note that I worked out the Lat/Lon shift for the Proj4 string by converting a 0,0 site coordinate to MGA, and used QGIS to convert this to WGS Lat/Lon. This give the origin Lat, Lon as per the other stackexchange post.

Emga,Nmga,Eplant,Nplant,Lon,Lat,Zmga,Zplant
525631.34,5258107.55,-515.61,31.13,147.3135953,-42.82917059,10.7,10.02
525634.87,5258111.13,-515.49,36.17,147.3136383,-42.82913824,10.33,9.65
525633.12,5258109.03,-515.33,33.44,147.313617,-42.82915721,12.43,11.74
525633.98,5258109.93,-515.32,34.68,147.3136275,-42.82914907,9.76,9.08
525692.84,5258084.4,-454.95,56.47,147.3143487,-42.829377,3.09,2.4
525692.41,5258083.45,-454.61,55.49,147.3143435,-42.82938557,3.14,2.46
525856.42,5258016.48,-289.24,119.25,147.3163532,-42.82998312,3.2,2.52
525876.22,5257998.88,-262.762,120.021,147.3165963,-42.83014095,9.59,8.91
525907.25,5258001.29,-241.82,143.07,147.3169758,-42.83011819,3.99,3.3
525906.52,5257997.55,-239.79,139.84,147.3169671,-42.8301519,3.23,2.54
525912.12,5257997.14,-235.43,143.38,147.3170356,-42.8301554,3.17,2.49
525925.34,5257969.74,-207,132.5,147.3171986,-42.8304017,3.41,2.73
525939,5257962.35,-191.98,136.49,147.3173661,-42.83046778,2.89,2.21
525934.89,5257911.72,-160.24,96.81,147.3173181,-42.83092385,2.91,2.23
525933.72,5257861.72,-126.8,59.6,147.3173061,-42.83137416,2.91,2.23
525981.79,5257875.04,-100.93,102.28,147.3178936,-42.83125257,41.97,41.29
526006.26,5257860.79,-73.34,108.7,147.3181937,-42.83138007,8.59,7.9
526086.76,5257842.88,-2.44,150.89,147.3191794,-42.83153861,8.79,8.1
525788.67,5257525.66,-1.84,-284.59,147.3155468,-42.83440535,10.52,9.83
525788.87,5257456.4,45.82,-334.89,147.3155524,-42.83502905,10.47,9.79
526149.54,5257795.97,75.46,159.8,147.3199497,-42.8319589,10.48,9.8
525789.24,5257402.47,83.09,-373.9,147.3155594,-42.83551469,42.49,41.81
526184.03,5257761.77,124.03,158.57,147.3203733,-42.8322657,37.11,36.43

Answer 1

To solve the problem half-manually, I put your plant coordinates in a text file named local.txt, and let cs2cs convert them to EPSG:28355:

cs2cs +proj=omerc +lat_0=-42.8325458 +lonc=147.3179408 +alpha=43.084 +k_0=0.99994 +x_0=0 +y_0=0 +gamma=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +to +init=epsg:28355 -f "%%.3f" local.txt >out.txt

Furthermore, I took your complete file into Libre Office, paste the result from cs2cs in three new columns, let LibreOffice calculate the offset per column, and the squares and sum of those columns:

The result (in red) is not very good, as you experienced.

Then I started playing with the variables lat_0, lonc, alpha, k_0. Obviously the first two can be changed to minimize the mean value of deltE and deltN respectively (the yellow value marks the changed value, and the green the minimized result), Similarly, alpha and k_0 can be changed to reduce the square sum. You have to go through all variables several times to get the minimum of the square sums.

The resulting proj string

+proj=omerc +lat_0=-42.8325458 +lonc=147.3179408 +alpha=43.084 +k_0=0.99994 +x_0=0 +y_0=0 +gamma=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

has offsets of less than 10 cm. You might get better values with inputs covering millimeters.

Comparing to your parameters, the origin is quite good, but the rotation is different. This might be because UTM Zone 55 has the central meridian (true North) at 147° East, while the local system has the axis of true North at 147.3179408° East. So the MGA Y axis at your origin is rotated against true North a bit. I guess this is not obeyed by the python code you used.