I'm downloading this csv from this website. I'm having a lot of trouble with the coordinates in this file, they use easting and northing and so far I haven't found a working solution to covert those values into latitude longitude. I've already found this, but those functions return invalid values when I test them on this website. I've looked up various UTM maps and figured out London should be in zone 30, but the converted latitude, longitude are something like this:
I test the function with the first line of the 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)
}
Could someone please explain to me what I'm doing wrong and how to convert easting, northing to latitude and longitude properly?
I'm using PHP and MySQL for this code, the latitude and longitude shoudl be stores in a database and late rbe used to display markers on Google Maps.
I've tried these functions:
////////////////////////////////////////////////////////////////////////////////////////////
//
// 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
}
And I've used this code:
<?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: pdana@pdana.com, 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
?>
