2

I have been trying to create a WKT entry to use in gdalwarp etc. I need an affine transform from coordinates in MGA94 Zone 54 to our local mine grid. All efforts have so far failed, including using the bits and pieces I have found here.

I have a MapInfow.prj file entry that works beautifully but I need to convert our imagery from MGA to mine grid to supply to mining consultants. This entry is below with the affine parameters highlighted:

"SRG/AISTR", 3008, 116,7, 141, 0, 0.9996, 500000, 10000000, 7,0.999770876952, -0.031210618094, -193025.578125, 0.031224025413, 0.999768821523, -7589075.9140625, -10000, -10000, 500000, 500000

The standard scale, rotate parameters as supplied by our surveyors are here:

Rotation: 358.2111111 degrees
Scale: 0.999741893
Local East: 16000
Local North: 20000
MGAz54 E: 446231.629
MGAz54 N: 7596899.020

My latest attempt using the affine parameters from MapInfo are here:

FITTEDCS["SRG", 
   PARAM_MT["Affine",
       PARAMETER["num_row, 3], 
       PARAMETER["num_row, 3], 
       PARAMETER["elt_0_0", -193025.578125],
       PARAMETER["elt_0_1", -7589075.9140625], 
       PARAMETER["elt_1_0", 0.999770876952], 
       PARAMETER["elt_1_1", 0.031224025413], 
       PARAMETER["elt_2_0", -0.031210618094], 
       PARAMETER["elt_2_1", 0.999768821523], 
       PARAMETER["Central Meridian", 140.4800753], 
       PARAMETER["Latitude_Of_Origin", -21.73061814]],  
   PROJCS["GDA_1994_MGA_Zone_54",
       GEOGCS["GCS_GDA_1994",DATUM["D_GDA_1994",
            SPHEROID["GRS_1980",6378137.0,298.257222101]],
       PRIMEM["Greenwich",0.0],
       UNIT["Degree",0.0174532925199433]],
   PROJECTION["Transverse_Mercator"],
   PARAMETER["False_Easting",500000.0],
   PARAMETER["False_Northing",10000000.0],
   PARAMETER["Central_Meridian",141.0],
   PARAMETER["Scale_Factor",0.9996],
   PARAMETER["Latitude_Of_Origin",0.0],
   UNIT["Meter",1.0],AUTHORITY["EPSG",28354]]

As requested by Andre here are some control points:

Local East    Local North MGAz54 East MGAz54 North
15896 27965   446376.2852 7604861.329
12458 28452   442956.0464 7605455.262
35698 27412   466146.2673 7603690.743
15698 29654   446231.1445 7606555.249
16547 23456   446886.0792 7600335.373
  • can you supply a few sample coords in both systems? – AndreJ Feb 21 '14 at 5:11
  • Sure thing...16444.297 27316.822 1384.4 446903.945 7604196.522 15896 27965 446376.285 7604861.329 12458 28452 442956.046 7605455.262 35698 27412 466146.267 7603690.743 15698 29654 446231.145 7606555.249 16547 23456 446886.079 7600335.373 – Tonyb Feb 21 '14 at 5:12
  • Apologies for the formatting Andre. The enter key didn't work as expected. – Tonyb Feb 21 '14 at 5:14
  • 16444.297, 27316.822 - 446903.945, 7604196.522 15896, 27965 - 446376.285, 7604861.329 12458,28452 - 442956.046, 7605455.262 35698, 27412 - 466146.267, 7603690.743 15698, 29654 - 446231.145, 7606555.249 16547, 23456 - 446886.079, 7600335.373 – Tonyb Feb 21 '14 at 5:18
2

A first try of mine:

+proj=omerc +lat_0=-21.73061814 +lonc=140.4800753 +alpha=-358.2111111 +k=0.9997418930000001 +x_0=16000 +y_0=20000 +gamma=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

using the parameters you gave, but not yet perfect.

For the Affine implementation, see also:

PROJ.4 custom projection that is Transverse Mercator with Affine post-process

GRASS v.transform should do the affine transformation, but I have not yet used it:

http://grass.osgeo.org/grass64/manuals/v.transform.html


EDIT

an improved version, after doing some statistics:

+proj=omerc +lat_0=-21.906057 +lonc=140.31864 +alpha=2.04 +k=1 +x_0=0 +y_0=0 +gamma=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

or as WKT:

PROJCS["Hotine_Oblique_Mercator_Azimuth_Center",
GEOGCS["GCS_GRS 1980(IUGG, 1980)", 
  DATUM["D_unknown",
  SPHEROID["GRS80",6378137,298.257222101]], 
  PRIMEM["Greenwich",0],
  UNIT["Degree",0.017453292519943295]],
PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],
  PARAMETER["latitude_of_center",-21.906057],
  PARAMETER["longitude_of_center",140.31864],
  PARAMETER["azimuth",2.04],
  PARAMETER["scale_factor",1], 
  PARAMETER["false_easting",0],
  PARAMETER["false_northing",0], 
  UNIT["Meter",1]]

which gives me less than 3 metres offset to the East for your most eastern point.

  • Thanks Andre, Just tried in gdalwarp and the resulting image was only about 4-5 metres east of true according when other data was overlaid. I'll continue to refine based on your version. Thanks a lot – Tonyb Feb 23 '14 at 23:04
  • You could apply a false Easting of +/- 5 meters, but that would spoil the GCP's you provided. Alternatively, try georeferemcing with some of your GCP to a vrt as an output. Then edit the vrt with a text editor, adding as much GCP as you want. – AndreJ Feb 24 '14 at 6:05
  • I had a good look at the original image and based on the GCP's it seems to be about 3.5 metres east of where it should be. Therefore the reprojected image based on your parameters is pretty good. I have since applied a false easting to both and the align well. Thanks – Tonyb Feb 25 '14 at 22:01

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.