Define custom CRS in WKT from point and angle
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
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.