8

I'm trying to define a custom CRS using the WKT Syntax. However when I do the projection I'm off by about 2km.

Here is my rotation point.

Local X and Y:
X: 4635.396 Y: 2397.085

MGA94 Zone50:
x: 560255.527 y: 7427753.462

Control Points:

Mine X | Mine Y| MGA 94(50) X| MGA 94(50) Y
2453.122|3210.002|563053.406|7431461.771
-1735.225|-853.24|557798.872|7428929.256
5663.648|7386.58|567416.171|7434410.368
12607.859|-1438.839|571218.306|7423848.605
8502.84|2620.24|568605.287|7428993.832
-2500.032|3457.767|558433.331|7433259.449

These are the steps I'm following base on WKT for local mine grid:

  1. Convert the MGA94 Zone50 (EPSG:28350) x and y to longitude and Latitude ("EPSG:4326"). I've used the python package pyproj
from pyproj import Transformer
transformer = Transformer.from_crs("EPSG:28350", "EPSG:4326", always_xy=True)
print(transformer.transform(564420.896, 7430150.547))

This gives the points (117.62970383981178, -23.236582623614485)

  1. Base on the link above I've put used the WKT synta
PROJCS["Hotine_Oblique_Mercator_Azimuth_Center",
GEOGCS["GCS_GRS 1980(IUGG, 1980)", 
DATUM["D_unknown",
SPHEROID["GRS80",6378137,298.257222101]], 
PRIMEM["Greenwich",0],
UNIT["Degree",0.017453292519943295]],
PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],
PARAMETER["latitude_of_center",-23.25839260829391]
PARAMETER["longitude_of_center", 117.58908484003899],
PARAMETER["azimuth",-18.39841101],
PARAMETER["scale_factor",0.999585495], 
PARAMETER["false_easting",0],
PARAMETER["false_northing",0], 
UNIT["Meter",1]]
  1. I then paste the code into QGIS custom CRS.

When I apply this custom CRS to a polygon layer in QGIS the polygon appears about 2km away from the actual location.

Can anyone offer any advice on how to achieve more accuracy?

7
  • 1
    I've had issues with the Hotine Oblique projection before, after significant research I found that one of the the parameters wasn't supported (azimuth I think, it was a while ago), I don't know if this has been fixed. There's an old post trac.osgeo.org/grass/ticket/1 which might help. Feb 27 '20 at 3:55
  • Do you have any other suggestions for projections to use? Feb 27 '20 at 5:46
  • 1
    This might help gis.stackexchange.com/questions/63107/… seeing as you have local x,y from MGA94/Zone50 x,y transformation matrix. Feb 27 '20 at 5:56
  • 1
    In the linked article, false Easting and Northing are not zero for the "first try". That might give the error when using your parameters.
    – AndreJ
    Feb 27 '20 at 6:37
  • @MichaelStimson I gave it a try, but the issue is that I can't get the rotation with PROJ4 or WKT parameters. Feb 27 '20 at 8:26
5

Update - See python script below for an answer

Original String (Red)

+proj=omerc +lat_0=-23.2583926082939 +lonc=117.589084840039 +alpha=-0 +gamma=0 +k=0.999585495 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs

gamma string by -18 (Green)

+proj=omerc +lat_0=-23.2583926082939 +lonc=117.589084840039 +alpha=-0 +gamma=-18 +k=0.999585495 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs

This results in a tilt in some axis: Alpha by -18

alpha string by -18 (Green)

+proj=omerc +lat_0=-23.2583926082939 +lonc=117.589084840039 +alpha=-18 +gamma=0 +k=0.999585495 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs

This results in another tilt:

enter image description here

So somewhere between these 4 parameters by using trial and error (or a python script) i should be able to figure this out.


EDIT: If anyone is curious I developed a nasty python script that lets you put an initial guess of coordinates and it finds the lowest error with the control points.

import pyproj
import math
import numpy as np
from statistics import mean
import scipy.optimize as optimize

#This function converts the numbers into text
def text_2_CRS(params):
    # print(params)  # <-- you'll see that params is a NumPy array
    x_0, y_0, gamma, alpha, lat_0, lonc = params # <-- for readability you may wish to assign names to the component variables
    pm = '+proj=omerc +lat_0='+ str(lat_0) +' +lonc='+ str(lonc) +' +alpha=' + str(alpha) + ' +gamma=' + str(
        gamma) + ' +k=0.999585495 +x_0=' + str(x_0) + ' +y_0=' + str(y_0) + ' +ellps=GRS80 +units=m +no_defs'
    return pm

#Optimisation function
def convert(params):
    pm = text_2_CRS(params)
    trans_points = []
    #Put your control points in mine grid coordinates here
    points_local = [[5663.648, 7386.58],
                    [20265.326, 493.126],
                    [1000, -10000],
                    [-1000, -10000],
                    [1331.817, 2390.206],
                    [5794, -1033.6],
                    ]
    # Put your control points here mga here
    points_mga = [[567416.145863305, 7434410.3451835],
                  [579090.883705669, 7423265.25196681],
                  [557507.390559793, 7419390.6658927],
                  [555610.407664593, 7420021.64968145],
                  [561731.125709093, 7431037.98474379],
                  [564883.285081307, 7426382.75146683],
                  ]
    for i in range(len(points_local)):
        #note that EPSG:28350 is MGA94 Zone 50
        trans = pyproj.transform(pyproj.Proj(pm), pyproj.Proj("EPSG:28350"), points_local[i][0], points_local[i][1])
        trans_points.append(trans)
    error = []
    #this finds the difference between the control points
    for i in range(len(points_mga)):
        x1 = trans_points[i][0]
        y1 = trans_points[i][1]
        x2 = points_mga[i][0]
        y2 = points_mga[i][1]
        error.append(math.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2))

    print("Current Params are: ")
    with np.printoptions(precision=3, suppress=True):
        print(params)
    print("Current average error is: " + str(mean(error)) + " meters")
    print("String to use is: " + pm)
    print('')

    return mean(error)


#Add your inital guess
x_0 = 950
y_0 = -1200
gamma = -18.39841101
alpha=-0
lat_0 = -23.2583926082939
lonc = 117.589084840039


#define your control points
points_local = [[5663.648,7386.58],
          [20265.326,493.126],
          [1000,-10000],
          [-1000,-10000],
          [1331.817,2390.206],
          [5794,-1033.6],
          ]

points_mga = [[567416.145863305,7434410.3451835],
          [579090.883705669,7423265.25196681],
          [557507.390559793,7419390.6658927],
          [555610.407664593,7420021.64968145],
          [561731.125709093,7431037.98474379],
          [564883.285081307,7426382.75146683],
          ]


params = [x_0, y_0, gamma,alpha, lat_0, lonc]

error = convert(params)

print(error)

result = optimize.minimize(convert, params, method='Powell')
if result.success:
    fitted_params = result.x
    print(fitted_params)
else:
    raise ValueError(result.message)

This leaves me the final Proj4 code of:

+proj=omerc +lat_0=-23.258566991042546 +lonc=117.58903931496924 +alpha=-0.00092995750016844 +gamma=-18.167694329590468 +k=0.999585495 +x_0=972.059643024533 +y_0=-1213.4486096382636 +ellps=GRS80 +units=m +no_defs

Second Edit: The comments below made me realize I can play with the scale -

+proj=omerc +lat_0=-23.258567543613964 +lonc=117.58903874790323 +alpha=-0.0009318714702833909 +gamma=-18.166493294460672 +k=1.0000628514828176 +x_0=969.710105681703 +y_0=-1213.4835412494535 +ellps=GRS80 +units=m +no_defs

I get an average error of 0.0645m

4
  • You are right, taking gamma instead of alpha leads to better values.
    – AndreJ
    Feb 27 '20 at 15:58
  • I suggest to clean the trial and error from your answer to give a good advice for future visitors; and accept your own answer.
    – AndreJ
    Feb 27 '20 at 19:22
  • I'll test out a few more local projections with the script to make sure it works - I'll provide an update Feb 28 '20 at 1:31
  • 1
    @AndreJ It works on a few other local projections, i'm happy this solution works. For someone who is not familiar with projections I find there is a real lack of good documentation on converting projection systems - or i'm not good at finding them. Feb 29 '20 at 7:48
2

You are almost there, here are my steps:

First, calculate from MGA to local using a plane rotation:

MineX = k ((MGAx-xo) cos phi + (MGAy-yo) sin phi)
MineY = k (-(MGAx-xo) sin phi + (MGAy-yo) cos phi)

with MGAx and MGAy as MGA coordinates. This works perfectly with

k = 1.0004
phi = -18.4
xo = 559714
yo = 7429191

So now we have the center in MGA coordinates, and the angle in degrees.

Put the MGA coordinates in a text file and Convert the MGA to latlon with cs2cs:

cs2cs +init=epsg:28350 +to +init=epsg:4326 -f "%%.5f" <Paraburdoo-center.txt >out.txt

returns

117.58373   -23.24543 0.00000

From that, you can get the PROJ string:

+proj=omerc +lat_0=-23.24543 +lonc=117.58373 +alpha=18.4 +k=1 +x_0=0 +y_0=0 +gamma=0 +ellps=GRS80  +units=m +no_defs

And the sample coordinates in red, displayed at the MGA coordiantes, fit in a grid with the rotated CRS in blue:

enter image description here

Calculating all points, I still get offsets about 50 m.

Keep in mind that the given rotation is plane. The MGA Mercator cylinder is placed at the equator at 117°E, while the rotated Mercator cylinder is placed at 23° South.


In the Hotine definition, alpha is used to rotate the cylinder from true North, and gamma is used to rotated the plane coordinates back to North-up.

So, you can use a different approach: Leave the Mercator cylinder where MGA places it (117°E on the equator), and do the rotation with gamma only.

The local coordinates of 117°E are the false Easting and Northing, and can be calculated with MGAx=500000 and MGAy=1000000 in the formula above:

MineX = -868482
MineY = 2421499

with that, the PROJ string is:

+proj=omerc +lat_0=0 +lonc=117 +alpha=0 +gamma=-18.40009 +k=1.000006 +x_0=-868484 +y_0=2421498 +ellps=GRS80 +to_meter=1 +no_defs

k and gamma (and the false Easting/Northing as a follow-up) are adjusted to reduce distortion to less than 1 meter. You might adjust to_meter as well to get better values.

1
  • I used a script to "Brute force it" see the response above - I think I get less than a meter accuracy. I'll look into your formula. Feb 27 '20 at 14:23

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.