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

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.

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
``````
• Is it always ~2m? Or is it just off as you move away from the centre? Commented Sep 25, 2016 at 3:48
• Hey @BradHards, it's 2 m up the top left of the larger picture, but it's less and different all over. I realise the larger picture doesn't represent it well, but the yellow and green dots are in the two projections, displayed in UTM. Note that up the top left the difference is different than the lower right. THe single green dot at the bottom is at 0, 0 on the site datum. I can't see a systematic pattern to the distortions... Commented Sep 25, 2016 at 4:50
• Can you add how you incorporated the shift values into the custom CRS? And if possible include some coordinate pairs with minimal and maximal offset. I guess there is something wrong how you calculated the origin obeying the rotation. Commented Sep 25, 2016 at 5:35
• Hi @AndreJ, I've added the coordinates I used. As I say above, the lat/lon shift is calculated with the 0,0 origin of the site datum. The rotation is the rotation from the 4 parameter transform converted from radians to degrees. The scale I used is potentially not right, but I don't know how to derive it. Commented Sep 25, 2016 at 5:47
• I used some old Python code I had laying around that uses least squares to calculate helmert params (7 parameter). It's good code and I trust it. As I said, the 4 param transform works perfectly, and I tested that to and from MGA and SITE coords. Commented Sep 25, 2016 at 7:44

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.

• Thanks, @AndreJ, I can confirm that this is much better at <10 cm error and mostly much smaller. Maybe I'll brute force it to try to find a minimum... You've confirmed that there's no better Proj CRS to use, though, so I'll work from there. Commented Sep 25, 2016 at 23:06
• You could also try fixing lonc to 147° East (as UTM), and use x_0 and y_0 instead. Commented Sep 26, 2016 at 5:27
• I've had a go at exploring the set of potential parameters in lat, lon, scale and alpha space and you're parameters here are impressively good! I can't find better using python brute force... I'll have a look at x_0 and y_0. Commented Sep 26, 2016 at 5:59
• Well this gets fairly close using X and Y offsets, but it doesn't look like it will converge and I don't really like not nkowing what the parameters should be rather than just fiddling until they look good! +proj=omerc +lat_0=0 +lonc=147 +alpha=43.090 +k_0=0.99994 +x_0=-3215850 +y_0=3758550 +gamma=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs Commented Sep 26, 2016 at 7:04
• It seems to be better when the origin is close to the data points. That keeps any projection distortion small. But it was worth a try. Commented Sep 26, 2016 at 17:49