2

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. Custom Projection System Definition QGIS

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

17
1

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.

1
  • 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 '20 at 7:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.