Subversion Repositories eFlore/Projets.eflore-projets

Rev

Go to most recent revision | Blame | Compare with Previous | Last modification | View Log | RSS feed 

<?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();