| New file |
|---|
| 0,0 → 1,401 |
|---|
| <?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() |
| } |
| |
| } |
| |