WKT for local mine grid

Question

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

Answer 1

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.

---

Late Edit

A better solution using gamma instead of alpha can be found here:

Define custom CRS in WKT from point and angle