The relation between the Mine grid CRS and EPSG:28350 is not affine (presents variable errors). But if you want to define a CRS derived from EPSG:28350, we can do it just for 2-D affine operation methods.
The Mine grid is also not an Oblique Mercator projection, but I know that may be easier to define omerc CRSes. As a last resort, all depends of the tolerances admitted for your work.
In my opinion, the right way is "unproject" EPSG:28350 control points to geocentric coordinates and find the transformation parameters from the local Mine grid to geocentric CRS, but you will need to transform the data to a new CRS instead of trying to define the CRS of the data as it is. Also, elevation values are needed to get a better accuracy.
So, let me show my way to find the similarity (i.e., affine that preserves the shapes) transformation parameters that better adapts EPSG:28350 control points to their correlated pairs in the Mine grid CRS. I will use a Python module that I wrote: https://github.com/gabriel-de-luca/simil
import numpy as np
np.set_printoptions(precision=3,suppress=True)
import simil
# points list in [X,Y]
epsg_28350_points = [[563053.406, 7431461.771],
[557798.872, 7428929.256],
[567416.171, 7434410.368],
[571218.306, 7423848.605],
[568605.287, 7428993.832],
[558433.331, 7433259.449]]
mine_grid_points = [[2453.122, 3210.002],
[-1735.225, -853.24],
[5663.648, 7386.58],
[12607.859, -1438.839],
[8502.84, 2620.24],
[-2500.032, 3457.767]]
# include a zeros Z dimension
epsg_28350_points_z = np.concatenate((epsg_28350_points,
np.zeros(6).reshape(6,1)),
axis=1)
mine_grid_points_z = np.concatenate((mine_grid_points,
np.zeros(6).reshape(6,1)),
axis=1)
# find the similaity transformation parameters
m, r, t = simil.process(epsg_28350_points_z, mine_grid_points_z)
# transpose source points to get coords
epsg_28350_coords_z = np.array(epsg_28350_points_z).T
# get transformed coords
transformed_coords = m * r @ epsg_28350_coords_z + t
# print transformed points without Z dimension
print('Transfomed points = \n' + str(transformed_coords.T[...,:2]))
# scale the rotation matrix
mr = m * r
print()
# print the WKT2:2019 affine coefficients
print('A0 = ' + str(t[0][0]))
print('A1 = ' + str(mr[0][0]))
print('A2 = ' + str(mr[0][1]))
print('B0 = ' + str(t[1][0]))
print('B1 = ' + str(mr[1][0]))
print('B2 = ' + str(mr[1][1]))
Returns:
Transfomed points =
[[ 2453.217 3210.012]
[-1735.147 -853.188]
[ 5663.665 7386.61 ]
[12607.863 -1438.904]
[ 8502.757 2620.277]
[-2500.143 3457.703]]
A0 = 1814469.6967790825
A1 = 0.9492782454965082
A2 = -0.31575360293528776
B0 = -7229101.121198954
B1 = 0.31575360293528776
B2 = 0.9492782454965082
If the difference between translated coordinates and control points coordinates is admisible for your work, you can write now the WKT2:2019 CRS representation of a Derived (from EPSG:28350) Projected CRS:
DERIVEDPROJCRS["mine_derived",
BASEPROJCRS["GDA94 / MGA zone 54",
BASEGEOGCRS["GDA94",
DATUM["Geocentric Datum of Australia 1994",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1
]
]
],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433
]
],
ID["EPSG",4283
]
],
CONVERSION["Map Grid of Australia zone 54",
METHOD["Transverse Mercator",
ID["EPSG",9807
]
],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433
],
ID["EPSG",8801
]
],
PARAMETER["Longitude of natural origin",141,
ANGLEUNIT["degree",0.0174532925199433
],
ID["EPSG",8802
]
],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1
],
ID["EPSG",8805
]
],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1
],
ID["EPSG",8806
]
],
PARAMETER["False northing",10000000,
LENGTHUNIT["metre",1
],
ID["EPSG",8807
]
]
]
],
DERIVINGCONVERSION["Affine",
METHOD["Affine parametric transformation",
ID["EPSG",9624
]
],
PARAMETER["A0",1814469.6967790825,
LENGTHUNIT["metre",1
],
ID["EPSG",8623
]
],
PARAMETER["A1",0.9492782454965082,
SCALEUNIT["coefficient",1
],
ID["EPSG",8624
]
],
PARAMETER["A2",-0.31575360293528776,
SCALEUNIT["coefficient",1
],
ID["EPSG",8625
]
],
PARAMETER["B0",-7229101.121198954,
LENGTHUNIT["metre",1
],
ID["EPSG",8639
]
],
PARAMETER["B1",0.31575360293528776,
SCALEUNIT["coefficient",1
],
ID["EPSG",8640
]
],
PARAMETER["B2",0.9492782454965082,
SCALEUNIT["coefficient",1
],
ID["EPSG",8641
]
]
],
CS[Cartesian,2
],
AXIS["(E)",east,
ORDER[1
],
LENGTHUNIT["metre",1
]
],
AXIS["(N)",north,
ORDER[2
],
LENGTHUNIT["metre",1
]
]
]
This definition is valid and recognized by QGIS, and you can define a custom CRS with it.
What I think that you can't yet do is transform with GDAL to that CRS (but you can transform with the affine pipeline string if you want).
I was able to set this CRS to the points in QGIS and export to other CRS without problem.
"Mine - MGA", 3008, 33,7, 117, 0, 0.9996, 500000, 10000000 , 7, 0.949288, -0.315744, 1814394.28, 0.315737, 0.94926, -7228957.16, -100000, -100000, 100000, 100000
– Garth Cooper Feb 12 '20 at 06:00