New file |
0,0 → 1,338 |
<?php |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
/******************************************************************************* |
NAME NEW ZEALAND MAP GRID |
|
PURPOSE: Transforms input longitude and latitude to Easting and |
Northing for the New Zealand Map Grid projection. The |
longitude and latitude must be in radians. The Easting |
and Northing values will be returned in meters. |
|
|
ALGORITHM REFERENCES |
|
1. Department of Land and Survey Technical Circular 1973/32 |
http://www.linz.govt.nz/docs/miscellaneous/nz-map-definition.pdf |
|
2. OSG Technical Report 4.1 |
http://www.linz.govt.nz/docs/miscellaneous/nzmg.pdf |
|
|
IMPLEMENTATION NOTES |
|
The two references use different symbols for the calculated values. This |
implementation uses the variable names similar to the symbols in reference [1]. |
|
The alogrithm uses different units for delta latitude and delta longitude. |
The delta latitude is assumed to be in units of seconds of arc x 10^-5. |
The delta longitude is the usual radians. Look out for these conversions. |
|
The algorithm is described using complex arithmetic. There were three |
options: |
* find and use a Javascript library for complex arithmetic |
* write my own complex library |
* expand the complex arithmetic by hand to simple arithmetic |
|
This implementation has expanded the complex multiplication operations |
into parallel simple arithmetic operations for the real and imaginary parts. |
The imaginary part is way over to the right of the display; this probably |
violates every coding standard in the world, but, to me, it makes it much |
more obvious what is going on. |
|
The following complex operations are used: |
- addition |
- multiplication |
- division |
- complex number raised to integer power |
- summation |
|
A summary of complex arithmetic operations: |
(from http://en.wikipedia.org/wiki/Complex_arithmetic) |
addition: (a + bi) + (c + di) = (a + c) + (b + d)i |
subtraction: (a + bi) - (c + di) = (a - c) + (b - d)i |
multiplication: (a + bi) x (c + di) = (ac - bd) + (bc + ad)i |
division: (a + bi) / (c + di) = [(ac + bd)/(cc + dd)] + [(bc - ad)/(cc + dd)]i |
|
The algorithm needs to calculate summations of simple and complex numbers. This is |
implemented using a for-loop, pre-loading the summed value to zero. |
|
The algorithm needs to calculate theta^2, theta^3, etc while doing a summation. |
There are three possible implementations: |
- use pow in the summation loop - except for complex numbers |
- precalculate the values before running the loop |
- calculate theta^n = theta^(n-1) * theta during the loop |
This implementation uses the third option for both real and complex arithmetic. |
|
For example |
psi_n = 1; |
sum = 0; |
for (n = 1; n <=6; n++) { |
psi_n1 = psi_n * psi; // calculate psi^(n+1) |
psi_n = psi_n1; |
sum = sum + A[n] * psi_n; |
} |
|
|
TEST VECTORS |
|
NZMG E, N: 2487100.638 6751049.719 metres |
NZGD49 long, lat: 172.739194 -34.444066 degrees |
|
NZMG E, N: 2486533.395 6077263.661 metres |
NZGD49 long, lat: 172.723106 -40.512409 degrees |
|
NZMG E, N: 2216746.425 5388508.765 metres |
NZGD49 long, lat: 169.172062 -46.651295 degrees |
|
Note that these test vectors convert from NZMG metres to lat/long referenced |
to NZGD49, not the more usual WGS84. The difference is about 70m N/S and about |
10m E/W. |
|
These test vectors are provided in reference [1]. Many more test |
vectors are available in |
http://www.linz.govt.nz/docs/topography/topographicdata/placenamesdatabase/nznamesmar08.zip |
which is a catalog of names on the 260-series maps. |
|
|
EPSG CODES |
|
NZMG EPSG:27200 |
NZGD49 EPSG:4272 |
|
http://spatialreference.org/ defines these as |
Proj4php.defs["EPSG:4272"] = "+proj=longlat +ellps=intl +datum=nzgd49 +no_defs "; |
Proj4php.defs["EPSG:27200"] = "+proj=nzmg +lat_0=-41 +lon_0=173 +x_0=2510000 +y_0=6023150 +ellps=intl +datum=nzgd49 +units=m +no_defs "; |
|
|
LICENSE |
Copyright: Stephen Irons 2008 |
Released under terms of the LGPL as per: http://www.gnu.org/copyleft/lesser.html |
|
* ***************************************************************************** */ |
|
/** |
Initialize New Zealand Map Grip projection |
*/ |
class Proj4phpProjNzmg { |
|
/** |
* iterations: Number of iterations to refine inverse transform. |
* 0 -> km accuracy |
* 1 -> m accuracy -- suitable for most mapping applications |
* 2 -> mm accuracy |
*/ |
protected $iterations = 1; |
|
/** |
* |
*/ |
public function init() { |
$this->A = array( ); |
$this->A[1] = +0.6399175073; |
$this->A[2] = -0.1358797613; |
$this->A[3] = +0.063294409; |
$this->A[4] = -0.02526853; |
$this->A[5] = +0.0117879; |
$this->A[6] = -0.0055161; |
$this->A[7] = +0.0026906; |
$this->A[8] = -0.001333; |
$this->A[9] = +0.00067; |
$this->A[10] = -0.00034; |
|
$this->B_re = array( ); |
$this->B_im = array( ); |
$this->B_re[1] = +0.7557853228; |
$this->B_im[1] = 0.0; |
$this->B_re[2] = +0.249204646; |
$this->B_im[2] = +0.003371507; |
$this->B_re[3] = -0.001541739; |
$this->B_im[3] = +0.041058560; |
$this->B_re[4] = -0.10162907; |
$this->B_im[4] = +0.01727609; |
$this->B_re[5] = -0.26623489; |
$this->B_im[5] = -0.36249218; |
$this->B_re[6] = -0.6870983; |
$this->B_im[6] = -1.1651967; |
|
$this->C_re = array( ); |
$this->C_im = array( ); |
$this->C_re[1] = +1.3231270439; |
$this->C_im[1] = 0.0; |
$this->C_re[2] = -0.577245789; |
$this->C_im[2] = -0.007809598; |
$this->C_re[3] = +0.508307513; |
$this->C_im[3] = -0.112208952; |
$this->C_re[4] = -0.15094762; |
$this->C_im[4] = +0.18200602; |
$this->C_re[5] = +1.01418179; |
$this->C_im[5] = +1.64497696; |
$this->C_re[6] = +1.9660549; |
$this->C_im[6] = +2.5127645; |
|
$this->D = array( ); |
$this->D[1] = +1.5627014243; |
$this->D[2] = +0.5185406398; |
$this->D[3] = -0.03333098; |
$this->D[4] = -0.1052906; |
$this->D[5] = -0.0368594; |
$this->D[6] = +0.007317; |
$this->D[7] = +0.01220; |
$this->D[8] = +0.00394; |
$this->D[9] = -0.0013; |
} |
|
/** |
New Zealand Map Grid Forward - long/lat to x/y |
long/lat in radians |
*/ |
public function forward( $p ) { |
|
$lon = $p->x; |
$lat = $p->y; |
|
$delta_lat = $lat - $this->lat0; |
$delta_lon = $lon - $this->long0; |
|
// 1. Calculate d_phi and d_psi ... // and d_lambda |
// For this algorithm, delta_latitude is in seconds of arc x 10-5, so we need to scale to those units. Longitude is radians. |
$d_phi = $delta_lat / Proj4php::$common->SEC_TO_RAD * 1E-5; |
$d_lambda = $delta_lon; |
$d_phi_n = 1; // d_phi^0 |
|
$d_psi = 0; |
for( $n = 1; $n <= 10; $n++ ) { |
$d_phi_n = $d_phi_n * $d_phi; |
$d_psi = $d_psi + $this->A[$n] * $d_phi_n; |
} |
|
// 2. Calculate theta |
$th_re = $d_psi; |
$th_im = $d_lambda; |
|
// 3. Calculate z |
$th_n_re = 1; |
$th_n_im = 0; // theta^0 |
#$th_n_re1; |
#$th_n_im1; |
|
$z_re = 0; |
$z_im = 0; |
for( $n = 1; $n <= 6; $n++ ) { |
$th_n_re1 = $th_n_re * $th_re - $th_n_im * $th_im; |
$th_n_im1 = $th_n_im * $th_re + $th_n_re * $th_im; |
$th_n_re = $th_n_re1; |
$th_n_im = $th_n_im1; |
$z_re = $z_re + $this->B_re[$n] * $th_n_re - $this->B_im[$n] * $th_n_im; |
$z_im = $z_im + $this->B_im[$n] * $th_n_re + $this->B_re[$n] * $th_n_im; |
} |
|
// 4. Calculate easting and northing |
$p->x = ($z_im * $this->a) + $this->x0; |
$p->y = ($z_re * $this->a) + $this->y0; |
|
return $p; |
} |
|
/** |
New Zealand Map Grid Inverse - x/y to long/lat |
*/ |
public function inverse( $p ) { |
|
$x = $p->x; |
$y = $p->y; |
|
$delta_x = $x - $this->x0; |
$delta_y = $y - $this->y0; |
|
// 1. Calculate z |
$z_re = $delta_y / $this->a; |
$z_im = $delta_x / $this->a; |
|
// 2a. Calculate theta - first approximation gives km accuracy |
$z_n_re = 1; |
$z_n_im = 0; // z^0 |
$z_n_re1; |
$z_n_im1; |
|
$th_re = 0; |
$th_im = 0; |
for( $n = 1; $n <= 6; $n++ ) { |
$z_n_re1 = $z_n_re * $z_re - $z_n_im * $z_im; |
$z_n_im1 = $z_n_im * $z_re + $z_n_re * $z_im; |
$z_n_re = $z_n_re1; |
$z_n_im = $z_n_im1; |
$th_re = $th_re + $this->C_re[$n] * $z_n_re - $this->C_im[$n] * $z_n_im; |
$th_im = $th_im + $this->C_im[$n] * $z_n_re + $this->C_re[$n] * $z_n_im; |
} |
|
// 2b. Iterate to refine the accuracy of the calculation |
// 0 iterations gives km accuracy |
// 1 iteration gives m accuracy -- good enough for most mapping applications |
// 2 iterations bives mm accuracy |
for( $i = 0; $i < $this->iterations; $i++ ) { |
$th_n_re = $th_re; |
$th_n_im = $th_im; |
$th_n_re1; |
$th_n_im1; |
|
$num_re = $z_re; |
$num_im = $z_im; |
for( $n = 2; $n <= 6; $n++ ) { |
$th_n_re1 = $th_n_re * th_re - $th_n_im * $th_im; |
$th_n_im1 = $th_n_im * $th_re + $th_n_re * $th_im; |
$th_n_re = $th_n_re1; |
$th_n_im = $th_n_im1; |
$num_re = $num_re + ($n - 1) * ($this->B_re[$n] * $th_n_re - $this->B_im[$n] * $th_n_im); |
$num_im = $num_im + (n - 1) * ($this->B_im[$n] * $th_n_re + $this->B_re[$n] * $th_n_im); |
} |
|
$th_n_re = 1; |
$th_n_im = 0; |
$den_re = $this->B_re[1]; |
$den_im = $this->B_im[1]; |
for( $n = 2; $n <= 6; $n++ ) { |
$th_n_re1 = $th_n_re * $th_re - $th_n_im * $th_im; |
$th_n_im1 = $th_n_im * $th_re + $th_n_re * $th_im; |
$th_n_re = $th_n_re1; |
$th_n_im = $th_n_im1; |
$den_re = $den_re + $n * ($this->B_re[$n] * $th_n_re - $this->B_im[$n] * $th_n_im); |
$den_im = $den_im + $n * ($this->B_im[n] * $th_n_re + $this->B_re[$n] * $th_n_im); |
} |
|
// Complex division |
$den2 = $den_re * $den_re + $den_im * $den_im; |
$th_re = ($num_re * $den_re + $num_im * $den_im) / $den2; |
$th_im = ($num_im * $den_re - $num_re * $den_im) / $den2; |
} |
|
// 3. Calculate d_phi ... // and d_lambda |
$d_psi = $th_re; |
$d_lambda = $th_im; |
$d_psi_n = 1; // d_psi^0 |
|
$d_phi = 0; |
for( $n = 1; $n <= 9; $n++ ) { |
$d_psi_n = $d_psi_n * $d_psi; |
$d_phi = $d_phi + $this->D[$n] * $d_psi_n; |
} |
|
// 4. Calculate latitude and longitude |
// d_phi is calcuated in second of arc * 10^-5, so we need to scale back to radians. d_lambda is in radians. |
$lat = $this->lat0 + ($d_phi * Proj4php::$common->SEC_TO_RAD * 1E5); |
$lon = $this->long0 + $d_lambda; |
|
$p->x = $lon; |
$p->y = $lat; |
|
return $p; |
} |
|
} |
|
Proj4php::$proj['nzmg'] = new Proj4phpProjNzmg(); |