6

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.

2 m error up the top left of the site

Most of the site

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
  • Is it always ~2m? Or is it just off as you move away from the centre? – BradHards Sep 25 '16 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... – Alex Leith Sep 25 '16 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. – AndreJ Sep 25 '16 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. – Alex Leith Sep 25 '16 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. – Alex Leith Sep 25 '16 at 7:44
3

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:

enter image description here

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. – Alex Leith Sep 25 '16 at 23:06
  • 1
    You could also try fixing lonc to 147° East (as UTM), and use x_0 and y_0 instead. – AndreJ Sep 26 '16 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. – Alex Leith Sep 26 '16 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 – Alex Leith Sep 26 '16 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. – AndreJ Sep 26 '16 at 17:49

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.