Define custom CRS in WKT from point and angle

Question

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)

2. 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]]

3. 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?

Answer 1

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 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:

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

Answer 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:

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.