Subversion Repositories eFlore/Projets.eflore-projets

<?php

/**

 * Author : Julien Moquet

 *

 * Inspired by Proj4js from Mike Adair madairATdmsolutions.ca

 *                      and Richard Greenwood rich@greenwoodma$p->com

 * License: LGPL as per: http://www.gnu.org/copyleft/lesser.html

 */

/** datum object

 */

class proj4phpDatum {

    public $datum_type;

    public $datum_params;

    /**

     *

     * @param type $proj

     */

    public function __construct( $proj ) {

        $this->datum_type = Proj4php::$common->PJD_WGS84;   //default setting

        if( isset($proj->datumCode) && $proj->datumCode == 'none' ) {

            $this->datum_type = Proj4php::$common->PJD_NODATUM;

        }

        if( isset( $proj->datum_params ) ) {

            for( $i = 0; $i < sizeof( $proj->datum_params ); $i++ ) {

                $proj->datum_params[$i] = floatval( $proj->datum_params[$i] );

            }

            if( $proj->datum_params[0] != 0 || $proj->datum_params[1] != 0 || $proj->datum_params[2] != 0 ) {

                $this->datum_type = Proj4php::$common->PJD_3PARAM;

            }

            if( sizeof( $proj->datum_params ) > 3 ) {

                if( $proj->datum_params[3] != 0 || $proj->datum_params[4] != 0 ||

                    $proj->datum_params[5] != 0 || $proj->datum_params[6] != 0 ) {

                    $this->datum_type = Proj4php::$common->PJD_7PARAM;

                    $proj->datum_params[3] *= Proj4php::$common->SEC_TO_RAD;

                    $proj->datum_params[4] *= Proj4php::$common->SEC_TO_RAD;

                    $proj->datum_params[5] *= Proj4php::$common->SEC_TO_RAD;

                    $proj->datum_params[6] = ($proj->datum_params[6] / 1000000.0) + 1.0;

                }

            }

            $this->datum_params = $proj->datum_params;

        }

        if( isset( $proj ) ) {

            $this->a = $proj->a;    //datum object also uses these values

            $this->b = $proj->b;

            $this->es = $proj->es;

            $this->ep2 = $proj->ep2;

            #$this->datum_params = $proj->datum_params;

        }

    }

    /**

     *

     * @param type $dest

     * @return boolean Returns TRUE if the two datums match, otherwise FALSE.

     * @throws type

     */

    public function compare_datums( $dest ) {

        if( $this->datum_type != $dest->datum_type ) {

            return false; // false, datums are not equal

        } else if( $this->a != $dest->a || abs( $this->es - $dest->es ) > 0.000000000050 ) {

            // the tolerence for es is to ensure that GRS80 and WGS84

            // are considered identical

            return false;

        } else if( $this->datum_type == Proj4php::$common->PJD_3PARAM ) {

            return ($this->datum_params[0] == $dest->datum_params[0]

                      && $this->datum_params[1] == $dest->datum_params[1]

                      && $this->datum_params[2] == $dest->datum_params[2]);

        } else if( $this->datum_type == Proj4php::$common->PJD_7PARAM ) {

            return ($this->datum_params[0] == $dest->datum_params[0]

                      && $this->datum_params[1] == $dest->datum_params[1]

                      && $this->datum_params[2] == $dest->datum_params[2]

                      && $this->datum_params[3] == $dest->datum_params[3]

                      && $this->datum_params[4] == $dest->datum_params[4]

                      && $this->datum_params[5] == $dest->datum_params[5]

                      && $this->datum_params[6] == $dest->datum_params[6]);

        } else if( $this->datum_type == Proj4php::$common->PJD_GRIDSHIFT ||

            $dest->datum_type == Proj4php::$common->PJD_GRIDSHIFT ) {

            throw(new Exception( "ERROR: Grid shift transformations are not implemented." ));

            return false;

        }

        return true; // datums are equal

    }

    /*

     * The function Convert_Geodetic_To_Geocentric converts geodetic coordinates

     * (latitude, longitude, and height) to geocentric coordinates (X, Y, Z),

     * according to the current ellipsoid parameters.

     *

     *    Latitude  : Geodetic latitude in radians                     (input)

     *    Longitude : Geodetic longitude in radians                    (input)

     *    Height    : Geodetic height, in meters                       (input)

     *    X         : Calculated Geocentric X coordinate, in meters    (output)

     *    Y         : Calculated Geocentric Y coordinate, in meters    (output)

     *    Z         : Calculated Geocentric Z coordinate, in meters    (output)

     *

     */

    public function geodetic_to_geocentric( $p ) {

        $Longitude = $p->x;

        $Latitude = $p->y;

        $Height = isset( $p->z ) ? $p->z : 0;   //Z value not always supplied

        $Error_Code = 0;  //  GEOCENT_NO_ERROR;

        /*

         * * Don't blow up if Latitude is just a little out of the value

         * * range as it may just be a rounding issue.  Also removed longitude

         * * test, it should be wrapped by cos() and sin().  NFW for PROJ.4, Sep/2001.

         */

        if( $Latitude < -Proj4php::$common->HALF_PI && $Latitude > -1.001 * Proj4php::$common->HALF_PI ) {

            $Latitude = -Proj4php::$common->HALF_PI;

        } else if( $Latitude > Proj4php::$common->HALF_PI && $Latitude < 1.001 * Proj4php::$common->HALF_PI ) {

            $Latitude = Proj4php::$common->HALF_PI;

        } else if( ($Latitude < -Proj4php::$common->HALF_PI) || ($Latitude > Proj4php::$common->HALF_PI) ) {

            /* Latitude out of range */

            Proj4php::reportError( 'geocent:lat out of range:' . $Latitude );

            return null;

        }

        if( $Longitude > Proj4php::$common->PI )

            $Longitude -= (2 * Proj4php::$common->PI);

        $Sin_Lat = sin( $Latitude ); /*  sin(Latitude)  */

        $Cos_Lat = cos( $Latitude ); /*  cos(Latitude)  */

        $Sin2_Lat = $Sin_Lat * $Sin_Lat; /*  Square of sin(Latitude)  */

        $Rn = $this->a / (sqrt( 1.0e0 - $this->es * $Sin2_Lat )); /*  Earth radius at location  */

        $p->x = ($Rn + $Height) * $Cos_Lat * cos( $Longitude );

        $p->y = ($Rn + $Height) * $Cos_Lat * sin( $Longitude );

        $p->z = (($Rn * (1 - $this->es)) + $Height) * $Sin_Lat;

        return $Error_Code;

    }

    /**

     *

     * @param object $p

     * @return type

     */

    public function geocentric_to_geodetic( $p ) {

        /* local defintions and variables */

        /* end-criterium of loop, accuracy of sin(Latitude) */

        $genau = 1.E-12;

        $genau2 = ($genau * $genau);

        $maxiter = 30;

        $X = $p->x;

        $Y = $p->y;

        $Z = $p->z ? $p->z : 0.0;   //Z value not always supplied

        /*

        $P;        // distance between semi-minor axis and location

        $RR;       // distance between center and location

        $CT;       // sin of geocentric latitude

        $ST;       // cos of geocentric latitude

        $RX;

        $RK;

        $RN;       // Earth radius at location

        $CPHI0;    // cos of start or old geodetic latitude in iterations

        $SPHI0;    // sin of start or old geodetic latitude in iterations

        $CPHI;     // cos of searched geodetic latitude

        $SPHI;     // sin of searched geodetic latitude

        $SDPHI;    // end-criterium: addition-theorem of sin(Latitude(iter)-Latitude(iter-1))

        $At_Pole;     // indicates location is in polar region

        $iter;        // of continous iteration, max. 30 is always enough (s.a.)

        $Longitude;

        $Latitude;

        $Height;

        */

        $At_Pole = false;

        $P = sqrt( $X * $X + $Y * $Y );

        $RR = sqrt( $X * $X + $Y * $Y + $Z * $Z );

        /*      special cases for latitude and longitude */

        if( $P / $this->a < $genau ) {

            /*  special case, if P=0. (X=0., Y=0.) */

            $At_Pole = true;

            $Longitude = 0.0;

            /*  if (X,Y,Z)=(0.,0.,0.) then Height becomes semi-minor axis

             *  of ellipsoid (=center of mass), Latitude becomes PI/2 */

            if( $RR / $this->a < $genau ) {

                $Latitude = Proj4php::$common->HALF_PI;

                $Height = -$this->b;

                return;

            }

        } else {

            /*  ellipsoidal (geodetic) longitude

             *  interval: -PI < Longitude <= +PI */

            $Longitude = atan2( $Y, $X );

        }

        /* --------------------------------------------------------------

         * Following iterative algorithm was developped by

         * "Institut für Erdmessung", University of Hannover, July 1988.

         * Internet: www.ife.uni-hannover.de

         * Iterative computation of CPHI,SPHI and Height.

         * Iteration of CPHI and SPHI to 10**-12 radian res$p->

         * 2*10**-7 arcsec.

         * --------------------------------------------------------------

         */

        $CT = $Z / $RR;

        $ST = $P / $RR;

        $RX = 1.0 / sqrt( 1.0 - $this->es * (2.0 - $this->es) * $ST * $ST );

        $CPHI0 = $ST * (1.0 - $this->es) * $RX;

        $SPHI0 = $CT * $RX;

        $iter = 0;

        /* loop to find sin(Latitude) res$p-> Latitude

         * until |sin(Latitude(iter)-Latitude(iter-1))| < genau */

        do {

            ++$iter;

            $RN = $this->a / sqrt( 1.0 - $this->es * $SPHI0 * $SPHI0 );

            /*  ellipsoidal (geodetic) height */

            $Height = $P * $CPHI0 + $Z * $SPHI0 - $RN * (1.0 - $this->es * $SPHI0 * $SPHI0);

            $RK = $this->es * $RN / ($RN + $Height);

            $RX = 1.0 / sqrt( 1.0 - $RK * (2.0 - $RK) * $ST * $ST );

            $CPHI = $ST * (1.0 - $RK) * $RX;

            $SPHI = $CT * $RX;

            $SDPHI = $SPHI * $CPHI0 - $CPHI * $SPHI0;

            $CPHI0 = $CPHI;

            $SPHI0 = $SPHI;

        } while( $SDPHI * $SDPHI > $genau2 && $iter < $maxiter );

        /*      ellipsoidal (geodetic) latitude */

        $Latitude = atan( $SPHI / abs( $CPHI ) );

        $p->x = $Longitude;

        $p->y = $Latitude;

        $p->z = $Height;

        return $p;

    }

    /**

     * Convert_Geocentric_To_Geodetic

     * The method used here is derived from 'An Improved Algorithm for

     * Geocentric to Geodetic Coordinate Conversion', by Ralph Toms, Feb 1996

     *

     * @param object Point $p

     * @return object Point $p

     */

    public function geocentric_to_geodetic_noniter( $p ) {

        /*

        $Longitude;

        $Latitude;

        $Height;

        $W;        // distance from Z axis

        $W2;       // square of distance from Z axis

        $T0;       // initial estimate of vertical component

        $T1;       // corrected estimate of vertical component

        $S0;       // initial estimate of horizontal component

        $S1;       // corrected estimate of horizontal component

        $Sin_B0;   // sin(B0), B0 is estimate of Bowring aux variable

        $Sin3_B0;  // cube of sin(B0)

        $Cos_B0;   // cos(B0)

        $Sin_p1;   // sin(phi1), phi1 is estimated latitude

        $Cos_p1;   // cos(phi1)

        $Rn;       // Earth radius at location

        $Sum;      // numerator of cos(phi1)

        $At_Pole;  // indicates location is in polar region

        */

        $X = floatval( $p->x );  // cast from string to float

        $Y = floatval( $p->y );

        $Z = floatval( $p->z ? $p->z : 0 );

        $At_Pole = false;

        if( $X <> 0.0 ) {

            $Longitude = atan2( $Y, $X );

        } else {

            if( $Y > 0 ) {

                $Longitude = Proj4php::$common->HALF_PI;

            } else if( Y < 0 ) {

                $Longitude = -Proj4php::$common->HALF_PI;

            } else {

                $At_Pole = true;

                $Longitude = 0.0;

                if( $Z > 0.0 ) { /* north pole */

                    $Latitude = Proj4php::$common->HALF_PI;

                } else if( Z < 0.0 ) { /* south pole */

                    $Latitude = -Proj4php::$common->HALF_PI;

                } else { /* center of earth */

                    $Latitude = Proj4php::$common->HALF_PI;

                    $Height = -$this->b;

                    return;

                }

            }

        }

        $W2 = $X * $X + $Y * $Y;

        $W = sqrt( $W2 );

        $T0 = $Z * Proj4php::$common->AD_C;

        $S0 = sqrt( $T0 * $T0 + $W2 );

        $Sin_B0 = $T0 / $S0;

        $Cos_B0 = $W / $S0;

        $Sin3_B0 = $Sin_B0 * $Sin_B0 * $Sin_B0;

        $T1 = $Z + $this->b * $this->ep2 * $Sin3_B0;

        $Sum = $W - $this->a * $this->es * $Cos_B0 * $Cos_B0 * $Cos_B0;

        $S1 = sqrt( $T1 * $T1 + $Sum * $Sum );

        $Sin_p1 = $T1 / $S1;

        $Cos_p1 = $Sum / $S1;

        $Rn = $this->a / sqrt( 1.0 - $this->es * $Sin_p1 * $Sin_p1 );

        if( $Cos_p1 >= Proj4php::$common->COS_67P5 ) {

            $Height = $W / $Cos_p1 - $Rn;

        } else if( $Cos_p1 <= -Proj4php::$common->COS_67P5 ) {

            $Height = $W / -$Cos_p1 - $Rn;

        } else {

            $Height = $Z / $Sin_p1 + $Rn * ($this->es - 1.0);

        }

        if( $At_Pole == false ) {

            $Latitude = atan( $Sin_p1 / $Cos_p1 );

        }

        $p->x = $Longitude;

        $p->y = $Latitude;

        $p->z = $Height;

        return $p;

    }

    /************************************************************** */

    // pj_geocentic_to_wgs84( p )

    //  p = point to transform in geocentric coordinates (x,y,z)

    public function geocentric_to_wgs84( $p ) {

        if( $this->datum_type == Proj4php::$common->PJD_3PARAM ) {

            // if( x[io] == HUGE_VAL )

            //    continue;

            $p->x += $this->datum_params[0];

            $p->y += $this->datum_params[1];

            $p->z += $this->datum_params[2];

        } else if( $this->datum_type == Proj4php::$common->PJD_7PARAM ) {

            $Dx_BF = $this->datum_params[0];

            $Dy_BF = $this->datum_params[1];

            $Dz_BF = $this->datum_params[2];

            $Rx_BF = $this->datum_params[3];

            $Ry_BF = $this->datum_params[4];

            $Rz_BF = $this->datum_params[5];

            $M_BF = $this->datum_params[6];

            // if( x[io] == HUGE_VAL )

            //    continue;

            $p->x = $M_BF * ( $p->x - $Rz_BF * $p->y + $Ry_BF * $p->z) + $Dx_BF;

            $p->y = $M_BF * ( $Rz_BF * $p->x + $p->y - $Rx_BF * $p->z) + $Dy_BF;

            $p->z = $M_BF * (-$Ry_BF * $p->x + $Rx_BF * $p->y + $p->z) + $Dz_BF;

        }

    }

    /*************************************************************** */

    // pj_geocentic_from_wgs84()

    //  coordinate system definition,

    //  point to transform in geocentric coordinates (x,y,z)

    public function geocentric_from_wgs84( $p ) {

        if( $this->datum_type == Proj4php::$common->PJD_3PARAM ) {

            //if( x[io] == HUGE_VAL )

            //    continue;

            $p->x -= $this->datum_params[0];

            $p->y -= $this->datum_params[1];

            $p->z -= $this->datum_params[2];

        } else if( $this->datum_type == Proj4php::$common->PJD_7PARAM ) {

            $Dx_BF = $this->datum_params[0];

            $Dy_BF = $this->datum_params[1];

            $Dz_BF = $this->datum_params[2];

            $Rx_BF = $this->datum_params[3];

            $Ry_BF = $this->datum_params[4];

            $Rz_BF = $this->datum_params[5];

            $M_BF = $this->datum_params[6];

            $x_tmp = ($p->x - $Dx_BF) / $M_BF;

            $y_tmp = ($p->y - $Dy_BF) / $M_BF;

            $z_tmp = ($p->z - $Dz_BF) / $M_BF;

            //if( x[io] == HUGE_VAL )

            //    continue;

            $p->x = $x_tmp + $Rz_BF * $y_tmp - $Ry_BF * $z_tmp;

            $p->y = -$Rz_BF * $x_tmp + $y_tmp + $Rx_BF * $z_tmp;

            $p->z = $Ry_BF * $x_tmp - $Rx_BF * $y_tmp + $z_tmp;

        } //cs_geocentric_from_wgs84()

    }

}