Восток Север к широте

Я загружаю этот csv с этого веб-сайта. У меня большие проблемы с координатами в этом файле, они используют восток и север, и до сих пор я не нашел рабочего решения для преобразования этих значений в широту и долготу. Я уже нашел это, но эти функции возвращают неверные значения, когда я проверяю их на этом веб-сайте. Я просмотрел различные карты UTM и понял, что Лондон должен быть в зоне 30, но преобразованные широта и долгота выглядят примерно так:

Я тестирую функцию с первой строкой CSV:

TfL Reference (TADS area)   Camera Head Reference   Borough Easting Northing    Speed limit Location                                Fixed Camera Type   Site monitoring start date
CR-B-16012                  16/802013               Barkin  547752  186619      30mph       -  Valence Avenue  -  Becontree Avenue  Red Light           9/1/1994

Easting: 547752 Northing: 186619

array(2) {
  ["lat"]=>
  float(1.688348039115)
  ["lng"]=>
  float(-2.5706837360343)
}

Может ли кто-нибудь объяснить мне, что я делаю неправильно и как правильно преобразовать восточное и северное направление в широту и долготу?

Я использую PHP и MySQL для этого кода, широта и долгота должны храниться в базе данных, а затем использоваться для отображения маркеров на Картах Google.

Я пробовал эти функции:

////////////////////////////////////////////////////////////////////////////////////////////
//
// ToLL - function to compute Latitude and Longitude given UTM Northing and Easting in meters
//
//  Description:
//    This function converts input north and east coordinates (meters)
//    to the corresponding WGS84 Lat/Lon values relative to the defined
//    UTM zone.  
//
//  Parameters:
//    north   - (i) Northing (meters)
//    east    - (i) Easting (meters)
//    utmZone - (i) UTM Zone of the North and East parameters
//    lat     - (o) Latitude in degrees 
//    lon     - (o) Longitude in degrees
//
function ToLL(north,east,utmZone)
{ 
  // This is the lambda knot value in the reference
  var LngOrigin = DegToRad(utmZone * 6 - 183)

  // The following set of class constants define characteristics of the
  // ellipsoid, as defined my the WGS84 datum.  These values need to be
  // changed if a different dataum is used.    

  var FalseNorth = 0.   // South or North?
  //if (lat < 0.) FalseNorth = 10000000.  // South or North?
  //else          FalseNorth = 0.   

  var Ecc = 0.081819190842622       // Eccentricity
  var EccSq = Ecc * Ecc
  var Ecc2Sq = EccSq / (1. - EccSq)
  var Ecc2 = Math.sqrt(Ecc2Sq)      // Secondary eccentricity
  var E1 = ( 1 - Math.sqrt(1-EccSq) ) / ( 1 + Math.sqrt(1-EccSq) )
  var E12 = E1 * E1
  var E13 = E12 * E1
  var E14 = E13 * E1

  var SemiMajor = 6378137.0         // Ellipsoidal semi-major axis (Meters)
  var FalseEast = 500000.0          // UTM East bias (Meters)
  var ScaleFactor = 0.9996          // Scale at natural origin

  // Calculate the Cassini projection parameters

  var M1 = (north - FalseNorth) / ScaleFactor
  var Mu1 = M1 / ( SemiMajor * (1 - EccSq/4.0 - 3.0*EccSq*EccSq/64.0 -
    5.0*EccSq*EccSq*EccSq/256.0) )

  var Phi1 = Mu1 + (3.0*E1/2.0 - 27.0*E13/32.0) * Math.sin(2.0*Mu1)
    + (21.0*E12/16.0 - 55.0*E14/32.0)           * Math.sin(4.0*Mu1)
    + (151.0*E13/96.0)                          * Math.sin(6.0*Mu1)
    + (1097.0*E14/512.0)                        * Math.sin(8.0*Mu1)

  var sin2phi1 = Math.sin(Phi1) * Math.sin(Phi1)
  var Rho1 = (SemiMajor * (1.0-EccSq) ) / Math.pow(1.0-EccSq*sin2phi1,1.5)
  var Nu1 = SemiMajor / Math.sqrt(1.0-EccSq*sin2phi1)

  // Compute parameters as defined in the POSC specification.  T, C and D

  var T1 = Math.tan(Phi1) * Math.tan(Phi1)
  var T12 = T1 * T1
  var C1 = Ecc2Sq * Math.cos(Phi1) * Math.cos(Phi1)
  var C12 = C1 * C1
  var D  = (east - FalseEast) / (ScaleFactor * Nu1)
  var D2 = D * D
  var D3 = D2 * D
  var D4 = D3 * D
  var D5 = D4 * D
  var D6 = D5 * D

  // Compute the Latitude and Longitude and convert to degrees
  var lat = Phi1 - Nu1*Math.tan(Phi1)/Rho1 *
    ( D2/2.0 - (5.0 + 3.0*T1 + 10.0*C1 - 4.0*C12 - 9.0*Ecc2Sq)*D4/24.0
     + (61.0 + 90.0*T1 + 298.0*C1 + 45.0*T12 - 252.0*Ecc2Sq - 3.0*C12)*D6/720.0 )

  lat = RadToDeg(lat)

  var lon = LngOrigin + 
    ( D - (1.0 + 2.0*T1 + C1)*D3/6.0
      + (5.0 - 2.0*C1 + 28.0*T1 - 3.0*C12 + 8.0*Ecc2Sq + 24.0*T12)*D5/120.0) /     Math.cos(Phi1)

  lon = RadToDeg(lon)

  // Create a object to store the calculated Latitude and Longitude values
  var sendLatLon = new PC_LatLon(lat,lon)

  // Returns a PC_LatLon object
  return sendLatLon
}

И я использовал этот код:

<?php
/*------------------------------------------------------------------------------
** File:        gPoint.php
** Description: PHP class to convert Latitude & Longitude coordinates into
**              UTM & Lambert Conic Conformal Northing/Easting coordinates. 
** Version:     1.3
** Author:      Brenor Brophy
** Email:       brenor dot brophy at gmail dot com
** Homepage:    brenorbrophy.com 
**------------------------------------------------------------------------------
** COPYRIGHT (c) 2005, 2006, 2007, 2008 BRENOR BROPHY
**
** The source code included in this package is free software; you can
** redistribute it and/or modify it under the terms of the GNU General Public
** License as published by the Free Software Foundation. This license can be
** read at:
**
** http://www.opensource.org/licenses/gpl-license.php
**
** This program is distributed in the hope that it will be useful, but WITHOUT 
** ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 
** FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. 
**------------------------------------------------------------------------------
**
** Code for datum and UTM conversion was converted from C++ code written by
** Chuck Gantz (chuck dot gantz at globalstar dot com) from
** http://www.gpsy.com/gpsinfo/geotoutm/ This URL has many other references to
** useful information concerning conversion of coordinates.
**
** Rev History
** -----------------------------------------------------------------------------
** 1.0  08/25/2005  Initial Release
** 1.1  05/15/2006  Added software license language to header comments
**                  Fixed an error in the convertTMtoLL() method. The latitude
**                  calculation had a bunch of variables without $ symbols.
**                  Fixed an error in convertLLtoTM() method, The $this-> was
**                  missing in front of a couple of variables. Thanks to Bob
**                  Robins of Maryland for catching the bugs.
** 1.2  05/18/2007  Added default of NULL to $LongOrigin arguement in convertTMtoLL()
**                  and convertLLtoTM() to eliminate warning messages when the
**                  methods are called without a value for $LongOrigin.
** 1.3  02/21/2008  Fixed a bug in the distanceFrom method, where the input parameters
**                  were not being converted to radians prior to calculating the
**                  distance. Thanks to Enrico Benco for finding pointing it out.
*/
define ("meter2nm", (1/1852));
define ("nm2meter", 1852);

/*------------------------------------------------------------------------------
** class gPoint     ... for Geographic Point
**
** This class encapsulates the methods for representing a geographic point on the
** earth in three different coordinate systema. Lat/Long, UTM and Lambert Conic
** Conformal.
*/
class gPoint
{
/* Reference ellipsoids derived from Peter H. Dana's website- 
** http://www.colorado.edu/geography/gcraft/notes/datum/datum_f.html
** email: [email protected], web page: www.pdana.com
**
** Source:
** Defense Mapping Agency. 1987b. DMA Technical Report: Supplement to Department
** of Defense World Geodetic System 1984 Technical Report. Part I and II.
** Washington, DC: Defense Mapping Agency
*/
    var $ellipsoid = array(//Ellipsoid name, Equatorial Radius, square of eccentricity  
    "Airy"                  =>array (6377563, 0.00667054),
    "Australian National"   =>array (6378160, 0.006694542),
    "Bessel 1841"           =>array (6377397, 0.006674372),
    "Bessel 1841 Nambia"    =>array (6377484, 0.006674372),
    "Clarke 1866"           =>array (6378206, 0.006768658),
    "Clarke 1880"           =>array (6378249, 0.006803511),
    "Everest"               =>array (6377276, 0.006637847),
    "Fischer 1960 Mercury"  =>array (6378166, 0.006693422),
    "Fischer 1968"          =>array (6378150, 0.006693422),
    "GRS 1967"              =>array (6378160, 0.006694605),
    "GRS 1980"              =>array (6378137, 0.00669438),
    "Helmert 1906"          =>array (6378200, 0.006693422),
    "Hough"                 =>array (6378270, 0.00672267),
    "International"         =>array (6378388, 0.00672267),
    "Krassovsky"            =>array (6378245, 0.006693422),
    "Modified Airy"         =>array (6377340, 0.00667054),
    "Modified Everest"      =>array (6377304, 0.006637847),
    "Modified Fischer 1960" =>array (6378155, 0.006693422),
    "South American 1969"   =>array (6378160, 0.006694542),
    "WGS 60"                =>array (6378165, 0.006693422),
    "WGS 66"                =>array (6378145, 0.006694542),
    "WGS 72"                =>array (6378135, 0.006694318),
    "WGS 84"                =>array (6378137, 0.00669438));

    // Properties
    var $a;                                     // Equatorial Radius
    var $e2;                                    // Square of eccentricity
    var $datum;                                 // Selected datum
    var $Xp, $Yp;                               // X,Y pixel location
    var $lat, $long;                            // Latitude & Longitude of the point
    var $utmNorthing, $utmEasting, $utmZone;    // UTM Coordinates of the point
    var $lccNorthing, $lccEasting;              // Lambert coordinates of the point
    var $falseNorthing, $falseEasting;          // Origin coordinates for Lambert Projection
    var $latOfOrigin;                           // For Lambert Projection
    var $longOfOrigin;                          // For Lambert Projection
    var $firstStdParallel;                      // For lambert Projection
    var $secondStdParallel;                     // For lambert Projection

    // constructor
    function gPoint($datum='WGS 84')            // Default datum is WGS 84
    {
        $this->a = $this->ellipsoid[$datum][0];     // Set datum Equatorial Radius
        $this->e2 = $this->ellipsoid[$datum][1];    // Set datum Square of eccentricity
        $this->datum = $datum;                      // Save the datum
    }
    //
    // Set/Get X & Y pixel of the point (used if it is being drawn on an image)
    //
    function setXY($x, $y)
    {
        $this->Xp = $x; $this->Yp = $y;
    }
    function Xp() { return $this->Xp; }
    function Yp() { return $this->Yp; }
    //
    // Set/Get/Output Longitude & Latitude of the point
    //
    function setLongLat($long, $lat)
    {
        $this->long = $long; $this->lat = $lat;
    }
    function Lat() { return $this->lat; }
    function Long() { return $this->long; }
    function printLatLong() { printf("Latitude: %1.5f Longitude: %1.5f",$this->lat, $this->long); }
    //
    // Set/Get/Output Universal Transverse Mercator Coordinates
    //
    function setUTM($easting, $northing, $zone='')  // Zone is optional
    {
        $this->utmNorthing = $northing;
        $this->utmEasting = $easting;
        $this->utmZone = $zone;
    }
    function N() { return $this->utmNorthing; }
    function E() { return $this->utmEasting; }
    function Z() { return $this->utmZone; }
    function printUTM() { print( "Northing: ".(int)$this->utmNorthing.", Easting: ".(int)$this->utmEasting.", Zone: ".$this->utmZone); }
    //
    // Set/Get/Output Lambert Conic Conformal Coordinates
    //
    function setLambert($easting, $northing)
    {
        $this->lccNorthing = $northing;
        $this->lccEasting = $easting;
    }
    function lccN() { return $this->lccNorthing; }
    function lccE() { return $this->lccEasting; }
    function printLambert() { print( "Northing: ".(int)$this->lccNorthing.", Easting: ".(int)$this->lccEasting); }

//------------------------------------------------------------------------------
//
// Convert Longitude/Latitude to UTM
//
// Equations from USGS Bulletin 1532 
// East Longitudes are positive, West longitudes are negative. 
// North latitudes are positive, South latitudes are negative
// Lat and Long are in decimal degrees
// Written by Chuck Gantz- chuck dot gantz at globalstar dot com, converted to PHP by
// Brenor Brophy, brenor dot brophy at gmail dot com
//
// UTM coordinates are useful when dealing with paper maps. Basically the
// map will can cover a single UTM zone which is 6 degrees on longitude.
// So you really don't care about an object crossing two zones. You just get a
// second map of the other zone. However, if you happen to live in a place that
// straddles two zones (For example the Santa Babara area in CA straddles zone 10
// and zone 11) Then it can become a real pain having to have two maps all the time.
// So relatively small parts of the world (like say California) create their own
// version of UTM coordinates that are adjusted to conver the whole area of interest
// on a single map. These are called state grids. The projection system is the
// usually same as UTM (i.e. Transverse Mercator), but the central meridian
// aka Longitude of Origin is selected to suit the logitude of the area being
// mapped (like being moved to the central meridian of the area) and the grid
// may cover more than the 6 degrees of lingitude found on a UTM map. Areas
// that are wide rather than long - think Montana as an example. May still
// have to have a couple of maps to cover the whole state because TM projection
// looses accuracy as you move further away from the Longitude of Origin, 15 degrees
// is usually the limit.
//
// Now, in the case where we want to generate electronic maps that may be
// placed pretty much anywhere on the globe we really don't to deal with the
// issue of UTM zones in our coordinate system. We would really just like a
// grid that is fully contigious over the area of the map we are drawing. Similiar
// to the state grid, but local to the area we are interested in. I call this
// Local Transverse Mercator and I have modified the function below to also
// make this conversion. If you pass a Longitude value to the function as $LongOrigin
// then that is the Longitude of Origin that will be used for the projection.
// Easting coordinates will be returned (in meters) relative to that line of
// longitude - So an Easting coordinate for a point located East of the longitude
// of origin will be a positive value in meters, an Easting coordinate for a point
// West of the longitude of Origin will have a negative value in meters. Northings
// will always be returned in meters from the equator same as the UTM system. The
// UTMZone value will be valid for Long/Lat given - thought it is not meaningful
// in the context of Local TM. If a NULL value is passed for $LongOrigin
// then the standard UTM coordinates are calculated.
//
    function convertLLtoTM($LongOrigin = NULL)
    {
        $k0 = 0.9996;
        $falseEasting = 0.0;

        //Make sure the longitude is between -180.00 .. 179.9
        $LongTemp = ($this->long+180)-(integer)(($this->long+180)/360)*360-180; // -180.00 .. 179.9;
        $LatRad = deg2rad($this->lat);
        $LongRad = deg2rad($LongTemp);

        if (!$LongOrigin)
        { // Do a standard UTM conversion - so findout what zone the point is in
            $ZoneNumber = (integer)(($LongTemp + 180)/6) + 1;
            // Special zone for South Norway
            if( $this->lat >= 56.0 && $this->lat < 64.0 && $LongTemp >= 3.0 && $LongTemp < 12.0 ) // Fixed 1.1
                $ZoneNumber = 32;
            // Special zones for Svalbard
            if( $this->lat >= 72.0 && $this->lat < 84.0 ) 
            {
                if(      $LongTemp >= 0.0  && $LongTemp <  9.0 ) $ZoneNumber = 31;
                else if( $LongTemp >= 9.0  && $LongTemp < 21.0 ) $ZoneNumber = 33;
                else if( $LongTemp >= 21.0 && $LongTemp < 33.0 ) $ZoneNumber = 35;
                else if( $LongTemp >= 33.0 && $LongTemp < 42.0 ) $ZoneNumber = 37;
            }
            $LongOrigin = ($ZoneNumber - 1)*6 - 180 + 3;  //+3 puts origin in middle of zone
            //compute the UTM Zone from the latitude and longitude
            $this->utmZone = sprintf("%d%s", $ZoneNumber, $this->UTMLetterDesignator());
            // We also need to set the false Easting value adjust the UTM easting coordinate
            $falseEasting = 500000.0;
        }
        $LongOriginRad = deg2rad($LongOrigin);

        $eccPrimeSquared = ($this->e2)/(1-$this->e2);

        $N = $this->a/sqrt(1-$this->e2*sin($LatRad)*sin($LatRad));
        $T = tan($LatRad)*tan($LatRad);
        $C = $eccPrimeSquared*cos($LatRad)*cos($LatRad);
        $A = cos($LatRad)*($LongRad-$LongOriginRad);

        $M = $this->a*((1   - $this->e2/4       - 3*$this->e2*$this->e2/64  - 5*$this->e2*$this->e2*$this->e2/256)*$LatRad 
                            - (3*$this->e2/8    + 3*$this->e2*$this->e2/32  + 45*$this->e2*$this->e2*$this->e2/1024)*sin(2*$LatRad)
                                                + (15*$this->e2*$this->e2/256 + 45*$this->e2*$this->e2*$this->e2/1024)*sin(4*$LatRad) 
                                                - (35*$this->e2*$this->e2*$this->e2/3072)*sin(6*$LatRad));

        $this->utmEasting = ($k0*$N*($A+(1-$T+$C)*$A*$A*$A/6
                        + (5-18*$T+$T*$T+72*$C-58*$eccPrimeSquared)*$A*$A*$A*$A*$A/120)
                        + $falseEasting);

        $this->utmNorthing = ($k0*($M+$N*tan($LatRad)*($A*$A/2+(5-$T+9*$C+4*$C*$C)*$A*$A*$A*$A/24
                     + (61-58*$T+$T*$T+600*$C-330*$eccPrimeSquared)*$A*$A*$A*$A*$A*$A/720)));
        if($this->lat < 0)
            $this->utmNorthing += 10000000.0; //10000000 meter offset for southern hemisphere
    }
//
// This routine determines the correct UTM letter designator for the given latitude
// returns 'Z' if latitude is outside the UTM limits of 84N to 80S
// Written by Chuck Gantz- chuck dot gantz at globalstar dot com, converted to PHP by
// Brenor Brophy, brenor dot brophy at gmail dot com
//
    function UTMLetterDesignator()
    {   
        if((84 >= $this->lat) && ($this->lat >= 72)) $LetterDesignator = 'X';
        else if((72 > $this->lat) && ($this->lat >= 64)) $LetterDesignator = 'W';
        else if((64 > $this->lat) && ($this->lat >= 56)) $LetterDesignator = 'V';
        else if((56 > $this->lat) && ($this->lat >= 48)) $LetterDesignator = 'U';
        else if((48 > $this->lat) && ($this->lat >= 40)) $LetterDesignator = 'T';
        else if((40 > $this->lat) && ($this->lat >= 32)) $LetterDesignator = 'S';
        else if((32 > $this->lat) && ($this->lat >= 24)) $LetterDesignator = 'R';
        else if((24 > $this->lat) && ($this->lat >= 16)) $LetterDesignator = 'Q';
        else if((16 > $this->lat) && ($this->lat >= 8)) $LetterDesignator = 'P';
        else if(( 8 > $this->lat) && ($this->lat >= 0)) $LetterDesignator = 'N';
        else if(( 0 > $this->lat) && ($this->lat >= -8)) $LetterDesignator = 'M';
        else if((-8 > $this->lat) && ($this->lat >= -16)) $LetterDesignator = 'L';
        else if((-16 > $this->lat) && ($this->lat >= -24)) $LetterDesignator = 'K';
        else if((-24 > $this->lat) && ($this->lat >= -32)) $LetterDesignator = 'J';
        else if((-32 > $this->lat) && ($this->lat >= -40)) $LetterDesignator = 'H';
        else if((-40 > $this->lat) && ($this->lat >= -48)) $LetterDesignator = 'G';
        else if((-48 > $this->lat) && ($this->lat >= -56)) $LetterDesignator = 'F';
        else if((-56 > $this->lat) && ($this->lat >= -64)) $LetterDesignator = 'E';
        else if((-64 > $this->lat) && ($this->lat >= -72)) $LetterDesignator = 'D';
        else if((-72 > $this->lat) && ($this->lat >= -80)) $LetterDesignator = 'C';
        else $LetterDesignator = 'Z'; //This is here as an error flag to show that the Latitude is outside the UTM limits

        return($LetterDesignator);
    }
    function convertTMtoLL($LongOrigin = NULL)
    {
        $k0 = 0.9996;
        $e1 = (1-sqrt(1-$this->e2))/(1+sqrt(1-$this->e2));
        $falseEasting = 0.0;
        $y = $this->utmNorthing;

        if (!$LongOrigin)
        { // It is a UTM coordinate we want to convert
            sscanf($this->utmZone,"%d%s",$ZoneNumber,$ZoneLetter);
            if($ZoneLetter >= 'N')
                $NorthernHemisphere = 1;//point is in northern hemisphere
            else
            {
                $NorthernHemisphere = 0;//point is in southern hemisphere
                $y -= 10000000.0;//remove 10,000,000 meter offset used for southern hemisphere
            }
            $LongOrigin = ($ZoneNumber - 1)*6 - 180 + 3;  //+3 puts origin in middle of zone
            $falseEasting = 500000.0;
        }

//      $y -= 10000000.0;   // Uncomment line to make LOCAL coordinates return southern hemesphere Lat/Long
        $x = $this->utmEasting - $falseEasting; //remove 500,000 meter offset for longitude

        $eccPrimeSquared = ($this->e2)/(1-$this->e2);

        $M = $y / $k0;
        $mu = $M/($this->a*(1-$this->e2/4-3*$this->e2*$this->e2/64-5*$this->e2*$this->e2*$this->e2/256));

        $phi1Rad = $mu  + (3*$e1/2-27*$e1*$e1*$e1/32)*sin(2*$mu) 
                    + (21*$e1*$e1/16-55*$e1*$e1*$e1*$e1/32)*sin(4*$mu)
                    +(151*$e1*$e1*$e1/96)*sin(6*$mu);
        $phi1 = rad2deg($phi1Rad);

        $N1 = $this->a/sqrt(1-$this->e2*sin($phi1Rad)*sin($phi1Rad));
        $T1 = tan($phi1Rad)*tan($phi1Rad);
        $C1 = $eccPrimeSquared*cos($phi1Rad)*cos($phi1Rad);
        $R1 = $this->a*(1-$this->e2)/pow(1-$this->e2*sin($phi1Rad)*sin($phi1Rad), 1.5);
        $D = $x/($N1*$k0);

        $tlat = $phi1Rad - ($N1*tan($phi1Rad)/$R1)*($D*$D/2-(5+3*$T1+10*$C1-4*$C1*$C1-9*$eccPrimeSquared)*$D*$D*$D*$D/24
                        +(61+90*$T1+298*$C1+45*$T1*$T1-252*$eccPrimeSquared-3*$C1*$C1)*$D*$D*$D*$D*$D*$D/720); // fixed in 1.1
        $this->lat = rad2deg($tlat);

        $tlong = ($D-(1+2*$T1+$C1)*$D*$D*$D/6+(5-2*$C1+28*$T1-3*$C1*$C1+8*$eccPrimeSquared+24*$T1*$T1)
                        *$D*$D*$D*$D*$D/120)/cos($phi1Rad);
        $this->long = $LongOrigin + rad2deg($tlong);
    }
    function configLambertProjection ($falseEasting, $falseNorthing,
                                      $longOfOrigin, $latOfOrigin,
                                      $firstStdParallel, $secondStdParallel)
    {
        $this->falseEasting = $falseEasting;
        $this->falseNorthing = $falseNorthing;
        $this->longOfOrigin = $longOfOrigin;
        $this->latOfOrigin = $latOfOrigin;
        $this->firstStdParallel = $firstStdParallel;
        $this->secondStdParallel = $secondStdParallel;
    }
    function convertLLtoLCC()
    {
        $e = sqrt($this->e2);

        $phi    = deg2rad($this->lat);                      // Latitude to convert
        $phi1   = deg2rad($this->firstStdParallel);         // Latitude of 1st std parallel
        $phi2   = deg2rad($this->secondStdParallel);        // Latitude of 2nd std parallel
        $lamda  = deg2rad($this->long);                     // Lonitude to convert
        $phio   = deg2rad($this->latOfOrigin);              // Latitude of  Origin
        $lamdao = deg2rad($this->longOfOrigin);             // Longitude of  Origin

        $m1 = cos($phi1) / sqrt(( 1 - $this->e2*sin($phi1)*sin($phi1)));
        $m2 = cos($phi2) / sqrt(( 1 - $this->e2*sin($phi2)*sin($phi2)));
        $t1 = tan((pi()/4)-($phi1/2)) / pow(( ( 1 - $e*sin($phi1) ) / ( 1 + $e*sin($phi1) )),$e/2);
        $t2 = tan((pi()/4)-($phi2/2)) / pow(( ( 1 - $e*sin($phi2) ) / ( 1 + $e*sin($phi2) )),$e/2);
        $to = tan((pi()/4)-($phio/2)) / pow(( ( 1 - $e*sin($phio) ) / ( 1 + $e*sin($phio) )),$e/2);
        $t  = tan((pi()/4)-($phi /2)) / pow(( ( 1 - $e*sin($phi ) ) / ( 1 + $e*sin($phi ) )),$e/2);
        $n  = (log($m1)-log($m2)) / (log($t1)-log($t2));
        $F  = $m1/($n*pow($t1,$n));
        $rf = $this->a*$F*pow($to,$n);
        $r  = $this->a*$F*pow($t,$n);
        $theta = $n*($lamda - $lamdao);

        $this->lccEasting = $this->falseEasting + $r*sin($theta);
        $this->lccNorthing = $this->falseNorthing + $rf - $r*cos($theta);
    }
    function convertLCCtoLL()
    {
        $e = sqrt($this->e2);

        $phi1   = deg2rad($this->firstStdParallel);         // Latitude of 1st std parallel
        $phi2   = deg2rad($this->secondStdParallel);        // Latitude of 2nd std parallel
        $phio   = deg2rad($this->latOfOrigin);              // Latitude of  Origin
        $lamdao = deg2rad($this->longOfOrigin);             // Longitude of  Origin
        $E      = $this->lccEasting;
        $N      = $this->lccNorthing;
        $Ef     = $this->falseEasting;
        $Nf     = $this->falseNorthing;

        $m1 = cos($phi1) / sqrt(( 1 - $this->e2*sin($phi1)*sin($phi1)));
        $m2 = cos($phi2) / sqrt(( 1 - $this->e2*sin($phi2)*sin($phi2)));
        $t1 = tan((pi()/4)-($phi1/2)) / pow(( ( 1 - $e*sin($phi1) ) / ( 1 + $e*sin($phi1) )),$e/2);
        $t2 = tan((pi()/4)-($phi2/2)) / pow(( ( 1 - $e*sin($phi2) ) / ( 1 + $e*sin($phi2) )),$e/2);
        $to = tan((pi()/4)-($phio/2)) / pow(( ( 1 - $e*sin($phio) ) / ( 1 + $e*sin($phio) )),$e/2);
        $n  = (log($m1)-log($m2)) / (log($t1)-log($t2));
        $F  = $m1/($n*pow($t1,$n));
        $rf = $this->a*$F*pow($to,$n);
        $r_ = sqrt( pow(($E-$Ef),2) + pow(($rf-($N-$Nf)),2) );
        $t_ = pow($r_/($this->a*$F),(1/$n));
        $theta_ = atan(($E-$Ef)/($rf-($N-$Nf)));

        $lamda  = $theta_/$n + $lamdao;
        $phi0   = (pi()/2) - 2*atan($t_);
        $phi1   = (pi()/2) - 2*atan($t_*pow(((1-$e*sin($phi0))/(1+$e*sin($phi0))),$e/2));
        $phi2   = (pi()/2) - 2*atan($t_*pow(((1-$e*sin($phi1))/(1+$e*sin($phi1))),$e/2));
        $phi    = (pi()/2) - 2*atan($t_*pow(((1-$e*sin($phi2))/(1+$e*sin($phi2))),$e/2));

        $this->lat  = rad2deg($phi);
        $this->long = rad2deg($lamda);
    }
    function distanceFrom($lon1, $lat1)
    { 
        $lon1 = deg2rad($lon1); $lat1 = deg2rad($lat1); // Added in 1.3
        $lon2 = deg2rad($this->Long()); $lat2 = deg2rad($this->Lat());

        $theta = $lon2 - $lon1;
        $dist = acos(sin($lat1) * sin($lat2) + cos($lat1) * cos($lat2) * cos($theta));

//      Alternative formula supposed to be more accurate for short distances
//      $dist = 2*asin(sqrt( pow(sin(($lat1-$lat2)/2),2) + cos($lat1)*cos($lat2)*pow(sin(($lon1-$lon2)/2),2)));
        return ( $dist * 6366710 ); // from http://williams.best.vwh.net/avform.htm#GCF
    }
    function distanceFromTM(&$pt)
    { 
        $E1 = $pt->E();     $N1 = $pt->N();
        $E2 = $this->E();   $N2 = $this->N();

        $dist = sqrt(pow(($E1-$E2),2)+pow(($N1-$N2),2));
        return $dist; 
    }
    function gRef($rX, $rY, $rE, $rN, $Scale, $LongOrigin)
    {
        $this->convertLLtoTM($LongOrigin);
        $x = (($this->E() - $rE) / $Scale)      // The easting in meters times the scale to get pixels
                                                // is relative to the center of the image so adjust to
            + ($rX);                            // the left coordinate.
        $y = $rY -                              // Adjust to bottom coordinate.
            (($rN - $this->N()) / $Scale);      // The northing in meters
                                                // relative to the equator. Subtract center point northing
                                                // to get relative to image center and convert meters to pixels
        $this->setXY((int)$x,(int)$y);          // Save the geo-referenced result.
    }
} // end of class gPoint

?>