Custom CRS for historic "site" grid

Question

I have a set of historic drawings which show a local site grid. I'd like to setup a custom CRS for this grid in QGIS so i can convert layers to/from this grid and MGA94 Zone 56 (EPSG:28356).

Presumably the drawing is shown with north to the top of the page, but after georeferencing it's rotated approx 9.4 deg clockwise.

I've read some grid intersections points from the drawing

Point 1: local xy = 7950, 1000 MGA94 Zone 56 = 383276.37, 6359888.50

Point 2: local xy = 7650, 850 MGA94 Zone 56 = 382941.47, 6359791.40

Point 3: Local xy = 7850, 850 MGA94 Zone 56 = 383155.14, 6359756.31

I've spent the day trying to wrap my head around:

but can't make anything work. Can anyone help me out?

It sounds fairly familiar, you have enough information to perform a control point reference with gdaltransform https://gdal.org/programs/gdaltransform.html. It's not unusual for sites to use a local CRS, sometimes it's easy like UTM with a more local false easting and northing to make the numbers smaller, other cases, like this sounds, are set by the initial chalk line somebody laid out on the site (evident by 9 degree rotation) that's roughly E-W or N-S as that was considered close enough at the time. — Michael Stimson, Jul 02 '20 at 05:40

Gabriel De Luca · Accepted Answer · 2020-07-02T16:01:13.120

First, some previous considerations and decisions: Both systems are 2D. EPSG:28356 is a projected (Transverse Mercator method) one. About the local reference system we don't know its procedence, but we can think that it is just a local cartesian 2D system.

We have some options, but let me consider two: One is transform the EPSG:28356 coordinates to geocentric, determine the paremeters of a similitude transformation from the local cartesian system to the geocentric one, in 3D, and then transform the coordinates of the local cartesian system to geocentric and reproject them to EPSG:28356. The other is consider a 2D similitude transformation between the local cartesian system and the projected one. This method will be less precise but has an advantage, we can define a CRS WKT for the 2D affine transformed from the projected reference system, and we can just define that CRS for the local layer.

Since you want to setup a custom CRS we will take the second option. So, lets find the 2D parameters to transform the coordinates from EPSG:28356 to the local CRS.

I will use a Python module that I wrote when I need to do that (https://github.com/gabriel-de-luca/simil). It is made for 3D, so I will create a Z=zero coordinate. There are other ways to get the parameters, but I use this:

import numpy as np
np.set_printoptions(precision=3,suppress=True)
import simil
source_points = [[383276.37, 6359888.50, 0],
                 [382941.47, 6359791.40, 0],
                 [383155.14, 6359756.31, 0]]
target_points = [[7950, 1000, 0],
                 [7650, 850, 0],
                 [7850, 850, 0]]
source_points_array = np.array(source_points)
target_points_array = np.array(target_points)
Get the parameters
m_scalar, r_matrix, t_vector = simil.process(source_points_array,target_points_array)
Print the parameters
print('\n m scalar = \n' + str(m_scalar))
print('\n R matrix = \n' + str(r_matrix))
print('\n T vector = \n' + str(t_vector))
Print them ready for the WKT
print('\n A0 = ' + str(t_vector[0]))
print('\n A1 = ' + str(m_scalar * r_matrix[0][0]))
print('\n A2 = ' + str(m_scalar * r_matrix[0][1]))
print('\n B0 = ' + str(t_vector[1]))
print('\n B1 = ' + str(m_scalar * r_matrix[1][0]))
print('\n B2 = ' + str(m_scalar * r_matrix[1][1]))

Which returns:

 m scalar = 
0.9598327695807208
R matrix = 
[[ 0.985 -0.172  0.   ]
 [ 0.172  0.985  0.   ]
 [ 0.     0.     1.   ]]
T vector = 
[[  694357.794]
 [-6075863.374]
 [       0.   ]]
A0 = [694357.794]
A1 = 0.9455596465082436
A2 = -0.1649117959886513
B0 = [-6075863.374]
B1 = 0.1649117959886513
B2 = 0.9455596465082436

Now, we can create a custom CRS WKT2:2019 definition, which will be a derived (Affine method) from EPSG:28356:

DERIVEDPROJCRS["Historic site grid",
    BASEPROJCRS["GDA94 / MGA zone 56",
        BASEGEOGCRS["GDA94",
            DATUM["Geocentric Datum of Australia 1994",
                ELLIPSOID["GRS 1980",6378137,298.257222101,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]]],
        CONVERSION["Map Grid of Australia zone 56",
            METHOD["Transverse Mercator",
                ID["EPSG",9807]],
            PARAMETER["Latitude of natural origin",0,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8801]],
            PARAMETER["Longitude of natural origin",153,
                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",694357.794,
            LENGTHUNIT["metre",1],
            ID["EPSG",8623]],
        PARAMETER["A1",0.945559646508244,
            SCALEUNIT["coefficient",1],
            ID["EPSG",8624]],
        PARAMETER["A2",-0.164911795988651,
            SCALEUNIT["coefficient",1],
            ID["EPSG",8625]],
        PARAMETER["B0",-6075863.374,
            LENGTHUNIT["metre",1],
            ID["EPSG",8639]],
        PARAMETER["B1",0.164911795988651,
            SCALEUNIT["coefficient",1],
            ID["EPSG",8640]],
        PARAMETER["B2",0.945559646508244,
            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]]]

That's all. Create the new Custom CRS:

and set it for the local system layer:

thank you @Gabriel De Luca. Great answer and it's given me a workable outcome out-of-the-box. If this was eBay I'd give you an A++++++++++ :) — timtim, Jul 02 '20 at 21:40
FWIW, I did have another shot at trying the multiple successive transformations, which I think is what @gabriel-de-luca you described in the other post (https://gis.stackexchange.com/questions/351394/reprojecting-local-crs-using-proj) making use of a PROJ pipeline. I feel I eventually understood each step in the process, but I couldn't make all the PROJ steps work together to give me a result. The WKT approach is accurate enough for 50 year old sketch maps! — timtim, Jul 02 '20 at 21:55
@timtim Hi, that is because we are defining the transformation from EPSG:28356 to the local cartesian one, to get a derived from projected CRS. You can set this CRS to the layer and export it to EPSG:28356, and the transformation applied is the inverse. If you want to transform the vector layer with ogr2ogr, you need to set a +inv parameter before the +proj=affine transformation, and the parameters are the same but named +xoff, +s11, +s12,+yoff, +s21 and +s22 instead of A1, A2, A3, B1, B2, B3, respectively. — Gabriel De Luca, Jul 02 '20 at 22:29
I also just want to thank you for this. I spent weeks trying to solve a similar problem. — Fiet Kleiner, Apr 09 '21 at 12:44

score 2 · Answer 2 · answered May 23 '21 at 18:25

You should use the tool "2D conformal transformation" of the "LF Tools" QGIS plugin.

First, load the points (You shoud use the "Table to point layer"). Observation: All the layers must bem in the EPSG:28356.
Second, create a line layer where each feature must have only 2 vertices (like a vector). Hint: Use the snapping tool to snap each point.
Lastly, run the tool ,as shown bellow: Look the result:

Custom CRS for historic "site" grid

2 Answers2

Get the parameters

Print the parameters

Print them ready for the WKT

Linked