# Define custom projection in QGIS from some control points

I have six control points in a local mine grid which I want to use to define a custom projection in QGIS. The reference points are in GDA94 / MGA zone 50 (EPSG:28350) (https://spatialreference.org/ref/epsg/gda94-mga-zone-50/)

Is there a way to do this in QGIS?

For clarity I want to define a CRS in the options below. Bonus points if someone can explain how to define Proj4 or WKT from custom points this would be even more useful.

Points below, Reference points is MGA 94 Zone 50:

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

Mapinfo reference sytle is here: "Mine to 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

• How can we know? This is just the distance matrix: i.stack.imgur.com/a7qU7.png, between each point and the others, in both systems. The differences are not constant. We don't have a numerical method to know the custom projection and parameters. You need to know something about it. How were the points measured? How do you get their coordinates? Feb 11, 2020 at 4:41
• From my understanding you can use the Proj4 or WKT and create a false Easting and Nothing from a known projection system. The points were surveyed. Is there any other parameters/information you need to know? Feb 11, 2020 at 4:57
• If you know the projection system please include it in your post. Feb 11, 2020 at 21:21
• From MGA to Mine... I am trying with the following WKT pastebin.com/raw/iit4NN4Z, QGIS (must be 3.10.2) recognizes it as valid but then can't reproject the layer. Help from the developers is needed I think. This features are form PROJ.6. This is the link to the Affine Transformation: proj.org/operations/transformations/affine.html. This is for the pipelines: proj.org/usage/transformation.html#transformation-pipelines Feb 12, 2020 at 7:52
• For a solution see gis.stackexchange.com/questions/352168/… Feb 27, 2020 at 19:24

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))
print('A1 = ' + str(mr))
print('A2 = ' + str(mr))
print('B0 = ' + str(t))
print('B1 = ' + str(mr))
print('B2 = ' + str(mr))
``````

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.

• Thanks for your response - as far as I can tell it's the most comprehensive guide I can find. I'll have a look at the code when I get some spare time and test a few points to see what average error I get but I'm sure it'll be spot on. Mar 16, 2020 at 7:13