This project aims to georeference historic maps which are of interest to bushwalkers in NSW, and make them available in a format which allows for easy viewing.
General parts involved are:
Here, twelve 1:63,360 scale military topographic maps were georeferenced and overlayed on the NSW LPI basemap. These digitised maps were sourced as low-quality JPGs (large TIFs also available) from the NLA (National Library of Australia) via Trove.
The maps were then georeferenced using the lat,long coordinates of corners of each map. These maps were generally published to cover 0.5 deg of long by 0.25 deg of lat, resulting in a regular grid of maps. Finally, they were clipped to cover only their mapped extent, removing the boards and information contained within.
No attempt was made to correct for the difference between the ANG and current GDA datums, so although the maps generally align with each other well, there can be discrepancies of 100s metres between the georeferenced maps and the LPI basemap.
Further work should be done in exploring variations of this process, particularly if a datum transformation can be applied.
When using old maps which have a CRS (coordinate reference system) significantly different from that which we use today, it can be useful to be able to transform coordinates from one CRS to another.
A brief overview of 'grids & datums' history for Australia can be found in this 2003 publication
Transformations listed before are quite well known, and many GIS packages can be used:
Note that these refer to geographic datums, rather than map grids of a projection
Geographic datum (short) | Map projection |
---|
Some details for pre AGD66 transformations can be found in the AGD technical manual (pp 52-54) which have been included below for reference.
Another way to transform coordinates is through a transformation grid. This is essentially a fancy version of a 'block shift' where the change (dLat and dLong or dE and dN) between the two CRSs was specified over a large region. An example of a rough block shift is adding 100m to eastings, and 200m to northings when transforming between AMG66 and MGA94.
A transformation grid uses the same principle but defines the shift at regular intervals throughout a region (eg, different values every minute). The NTv2 transformation grid is a specific file containing this information developed by Canada but now used throughout the world.
In many cases, you can find an NTv2 grid for the transformation you wish to perform. ICSM has produced multiple grids to transform between GDA94 and GDA202 which they have made freely available online. Grids for transforming between AGD66/84 and GDA94 also exist.
If a grid file does not exist for the transformation you wish to perform, it is possible (but appears difficult) to create your own. Some information about process flows was discussed in 2013 on the maptools mail list.
The Victorian land agency also appear to have a free tool perform transformations called GDAit
Killet Software has an NTvtTools package for developing and processing NTv2 files in both binary and ASCII formats. A trial version is freely available and a full licence currently costs 160 EUR.
A paper on the Generation of a NAD27-NAD83(CSRS) NTv2-type Grid Shift File for New Brunswick is also of interest.
The trial version of SEVENPAR by Killet Software was used to produce local parameters for a 7 parameter Helmann transformation. The trial version is limited to only 25 common points (locations known in both source and target CRS) so ones near Sydney sourced from the NGRS were used.
It's important to note the convention used (either Coordinate-Frame Rotation, or Position Vector) to define the three rotation parameters (signs are flipped between them).
The output in position vector format was:
Parameter | Value | Units |
---|---|---|
dx | 195.056526 | m |
dy | -112.959677 | m |
dz | -172.456394 | m |
rx | 3.59354054 | rad, seconds |
ry | 5.39693644 | rad, seconds |
rz | -1.34857381 | rad, seconds |
s | -3.2756552 | ppm |
Combining these with the known values of the ellipsoid parameters of the ANS, a custom CRS can be defined in PROJ to allow transform between ANG and GDA94 coordinates. This can then be used in QGIS to 'reproject' ANG georeferenced rasters into a GDA94 CRS.
The following PROJ string was used to define the ANG CRS:
+proj=longlat +a=6378339.78 +rf=294.26 +towgs84=195.056526,-112.959677,-172.456394,3.59354054,5.39693644,-1.34857381,-3.2756552 +no_defs
This CRS is then available to be applied to any raster.
PROJ is a useful program here.
Note that PROJ is being updated and dramatically changed - a useful overview here.
Handy tables of parameters can be found here
'Clarke Coordinates' or the 'ANG' (Australian National Grid) have the following parameters:
Parameter | Value | Units |
---|---|---|
ellipsoid | Clarke 1858 | |
a | 20,926,348 | british feet |
a | 6,378,339.78 | metres |
b | 20,855,233 | british feet |
b | 6,356,663.92 | metres |
1/f | 294.26 | |
projection | transverse mercator | |
false easting | 400,000 | yards |
false easting | 365,759.36 | metres |
false northing | 800,000 | yards |
false northing | 731,518.73 | metres |
zones | 1 to 8 | |
zone width | 5 | degrees |
zone 1 central meridian | 116E | degrees |
zone latitude origin | 34S | degrees |
scale factor | 1 |
On these maps, Sydney Observatory was used as a datum for triangulation, with the following values:
Format | Latitude | Longitude |
---|---|---|
DMS | S 33d 51m 41.10s | E 151d 12m 17.85s |
DD | S 33.861417d | E 151.204958d |
The NGDB provides the following values of the observatory in different datums:
Datum | Latitude | Longitude | Easting | Northing |
---|---|---|---|---|
Clarke | S33°51' 41.107“ | E151°12' 17.993” | 420746.19 | 816789.65 |
AGD66 | S33°51' 40.2924“ | E151°12' 12.3136” | 333809.38 | 6251769.51 |
AGD84 | S33°51' 40.3241“ | E151°12' 12.2557” | 333807.91 | 6251768.51 |
GDA94 | S33°51' 34.610232“ | E151°12' 16.470597” | 333913.76 | 6251959.42 |
Note:
proj -lu
, default is metres)+to_meter-parameter
instead and specify the value you multiple your unit by to get 1 metre (eg, +to_meter=100
would output centimetres)As an example, the below PROJ code will convert a lat,long (S 33.75, E 151.5) to E,N for zone 8 (cm=151e) of the ANG:
echo 151d30e 33d45s | proj +proj=tmerc +lat_0=34s +lon_0=151e +k_0=1 +a=6378339.78 +rf=294.26 +to_meter=0.91439841 +x_0=365759.36 +y_0=731518.73
Which returns: E=450,666.78 N=830,202.64
An explaination of parameters is as follows:
Zone Number | Western Longitude | Eastern Longitude | Zone Central Meridian |
---|---|---|---|
1 | 113d30e | 118d30e | 116de |
2 | 118d30e | 123d30e | 121de |
3 | 123d30e | 128d30e | 126de |
4 | 128d30e | 133d30e | 131de |
5 | 133d30e | 138d30e | 136de |
6 | 138d30e | 143d30e | 141de |
7 | 143d30e | 148d30e | 146de |
8 | 148d30e | 153d30e | 151de |