# Define custom CRS in WKT from point and angle

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?

• 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. Commented Feb 27, 2020 at 3:55
• Do you have any other suggestions for projections to use? Commented Feb 27, 2020 at 5:46
• This might help gis.stackexchange.com/questions/63107/… seeing as you have local x,y from MGA94/Zone50 x,y transformation matrix. Commented Feb 27, 2020 at 5:56
• In the linked article, false Easting and Northing are not zero for the "first try". That might give the error when using your parameters. Commented Feb 27, 2020 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. Commented Feb 27, 2020 at 8:26

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)

x_0 = 950
y_0 = -1200
gamma = -18.39841101
alpha=-0
lat_0 = -23.2583926082939
lonc = 117.589084840039

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

• You are right, taking gamma instead of alpha leads to better values. Commented Feb 27, 2020 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. Commented Feb 27, 2020 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 Commented Feb 28, 2020 at 1:31
• @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. Commented Feb 29, 2020 at 7:48

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.

• 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. Commented Feb 27, 2020 at 14:23