QGIS – Defining 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