 # Datum Conversion

Subject: Datum Conversion Changes necessary for datum conversion 1 Vick Airfull 2012-07-04 12:28:04

 1. Datum Conversion Reply Report abuse Vick Airfull - 2012-07-04 12:28:04
I've implemented a few new functions inside the exsiting code to provide datum conversion using the "Simple Three Parameter XYZ" (according to the equations and instructions found at http://www.colorado.edu/geography/gcraft/notes/datum/datum_f.html)

The basic idea is use the code as following (this example convert from ED50 to WGS84):

<?php
\$pt =& new gPoint("International"); //using "International" ellipsoid (for ED50 datum)
\$pt->setUTM( \$ex, \$ey, \$ez);
//Here you can use \$pt->E(), \$pt->N() and \$pt->Z() functions
\$pt->convertTMtoLL(); // Convert to Lat/Long in original Datum
// Here you can use \$pt->Lat() and \$pt->Long() functions
\$pt->convertLLtoECEF(); // Convert to XYZ in ECEF (Earth Centered Earth Fixed). First stage of conversion
\$pt->ChangeDATUMbyXYZ("ED 50","WGS 84"); //Convert Datum using "simple Three Parameter" method from ED50 (original) to WGS84 (new). Second stage
\$pt->convertECEFtoLL(); // Convert to Lat/long in new Datum. Third stage
// Here you can use \$pt->Lat() and \$pt->Long() functions
\$pt->convertLLtoTM(); //Convert to UTM in new Datum
//Here you can use \$pt->E(), \$pt->N() and \$pt->Z() functions
?>

The code only includes parameters (XYZ) for transformation from European_1950 datums, but it easy to add new transformations adding the parameters needed (check http://www.colorado.edu/geography/gcraft/notes/datum/datum_f.html) to the new defined \$xyz array in the code.

I hope it would help anyone so I would like to post my changes here but as I don�t know if I can upload a new gPoint.php file I will paste here the contents of the new gPoint.php (I've named it gPoint-1.4b.php):

<?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
**
**
** 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
** 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.
** 1.4b 07/04/2012 (changes by perdigano [beta])
** Included new methods and variables for datum conversion based
** on simple three parameter XYZ and geodetic coordinates, according
** to equations and instructions found at
**
*/
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-
** 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)
);

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

var \$xyz = array(//Datum, Ellipsoid, Delta a, Delta f (* 10^4), Delta X, Delta Y, Delta Z
"European_1950" =>array ("International", -251.0, -0.14192702, -87, -96, -120),
"European_1950_mean" =>array ("International", -251.0, -0.14192702, -87, -98, -121),
"ED 50" =>array ("International", -251.0, -0.14192702, -84, -107, -120),
"European_1979" =>array ("International", -251.0, -0.14192702, -86, -98, -119),
"WGS 84" =>array ("WGS 84", 0.0, 0.0, 0, 0, 0)
);

// Properties
var \$e2; // Square of eccentricity
var \$datum; // Selected datum
var \$Xp, \$Yp; // X,Y pixel location
var \$lat, \$long, \$h; // Latitude & Longitude of the point (and height)
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
var \$ecefX, \$ecefY, \$ecefZ; // For Earth Centered Earth Fixed XYZ Cartesian Coordinates

// constructor
function gPoint(\$datum='WGS 84') // Default datum is WGS 84
{
\$this->a = \$this->ellipsoid[\$datum]; // Set datum Equatorial Radius
\$this->e2 = \$this->ellipsoid[\$datum]; // 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; \$this->h = 0;
}

function setLongLatHeight(\$long, \$lat, \$h)
{
\$this->long = \$long; \$this->lat = \$lat; \$this->h = \$h;
}

function Lat() { return \$this->lat; }
function Long() { return \$this->long; }
function Height() { return \$this->h; }

function printLatLong() { printf("Latitude: %1.5f Longitude: %1.5f",\$this->lat, \$this->long); }
function printLatLongHeight() { printf("Latitude: %1.5f Longitude: %1.5f Height: %1.5f",\$this->lat, \$this->long, \$this->h); }

function xyzX() { return \$this->ecefX; }
function xyzY() { return \$this->ecefY; }
function xyzZ() { return \$this->ecefZ; }

function printXYZ() { printf("X: %1.5f Y: %1.5f Z: %1.5f",\$this->ecefX, \$this->ecefY, \$this->ecefZ); }

//
// 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 ECEF XYZ
//
//
function convertLLtoECEF()
{

//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;
\$Height = \$this->h;

\$this->ecefZ = (\$N * (1 - \$this->e2) + \$Height) * sin(\$LatRad);

}

//------------------------------------------------------------------------------
//
// Convert ECEF XYZ to Longitude/Latitude
//
//
function convertECEFtoLL()
{
\$X = \$this->ecefX;
\$Y = \$this->ecefY;
\$Z = \$this->ecefZ;

\$a = \$this->a;
\$e2 = \$this->e2;

\$b = \$a * sqrt(1-\$e2);
\$p = sqrt(\$X * \$X + \$Y * \$Y);
\$ee2 = (\$a * \$a - \$b * \$b) / (\$b * \$b);
\$phi = atan((\$Z * \$a) / (\$p * \$b));

\$this->lat = atan((\$Z + \$ee2 * \$b * sin(\$phi) * sin(\$phi) * sin(\$phi)) / (\$p - \$e2 * \$a * cos(\$phi) * cos(\$phi) * cos(\$phi)) );
\$this->long = atan2(\$Y,\$X);
\$N = \$this->a / sqrt(1-\$this->e2*sin(\$this->lat)*sin(\$this->lat));
\$this->h = (\$p / cos(\$this->lat)) - \$N;

}

//------------------------------------------------------------------------------
//
// Fix ECEF XYZ Datums
//
//
function ChangeDATUMbyXYZ(\$fromDatum, \$toDatum)
{
\$this->ecefX += \$this->xyz[\$fromDatum];
\$this->ecefY += \$this->xyz[\$fromDatum];
\$this->ecefZ += \$this->xyz[\$fromDatum];

\$this->ecefX -= \$this->xyz[\$toDatum];
\$this->ecefY -= \$this->xyz[\$toDatum];
\$this->ecefZ -= \$this->xyz[\$toDatum];

\$this->a = \$this->ellipsoid[\$this->xyz[\$toDatum]]; // Set datum Equatorial Radius
\$this->e2 = \$this->ellipsoid[\$this->xyz[\$toDatum]]; // Set datum Square of eccentricity
\$this->datum = \$this->xyz[\$toDatum]; // Save the datum

}

//------------------------------------------------------------------------------
//
// 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;

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;
}

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

\$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)

\$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);

+ (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);
}

//------------------------------------------------------------------------------
//
// Convert UTM to Longitude/Latitude
//
// 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
//
// If a value is passed for \$LongOrigin then the function assumes that
// a Local (to the Longitude of Origin passed in) Transverse Mercator
// coordinates is to be converted - not a UTM coordinate. This is the
// complementary function to the previous one. The function cannot
// tell if a set of Northing/Easting coordinates are in the North
// or South hemesphere - they just give distance from the equator not
// direction - so only northern hemesphere lat/long coordinates are returned.
// If you live south of the equator there is a note later in the code
// explaining how to have it just return southern hemesphere lat/longs.
//
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));

+ (21*\$e1*\$e1/16-55*\$e1*\$e1*\$e1*\$e1/32)*sin(4*\$mu)
+(151*\$e1*\$e1*\$e1/96)*sin(6*\$mu);

\$D = \$x/(\$N1*\$k0);

+(61+90*\$T1+298*\$C1+45*\$T1*\$T1-252*\$eccPrimeSquared-3*\$C1*\$C1)*\$D*\$D*\$D*\$D*\$D*\$D/720); // fixed in 1.1

\$tlong = (\$D-(1+2*\$T1+\$C1)*\$D*\$D*\$D/6+(5-2*\$C1+28*\$T1-3*\$C1*\$C1+8*\$eccPrimeSquared+24*\$T1*\$T1)
}

//------------------------------------------------------------------------------
// Configure a Lambert Conic Conformal Projection
//
// falseEasting & falseNorthing are just an offset in meters added to the final
// coordinate calculated.
//
// longOfOrigin & LatOfOrigin are the "center" latitiude and longitude of the
// area being projected. All coordinates will be calculated in meters relative
// to this point on the earth.
//
// firstStdParallel & secondStdParallel are the two lines of longitude (that
// is they run east-west) that define where the "cone" intersects the earth.
// Simply put they should bracket the area being projected.
//
//
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;
}

//------------------------------------------------------------------------------
//
// Convert Longitude/Latitude to Lambert Conic Easting/Northing
//
// This routine will convert a Latitude/Longitude coordinate to an Northing/
// Easting coordinate on a Lambert Conic Projection. The configLambertProjection()
// function should have been called prior to this one to setup the specific
// parameters for the projection. The Northing/Easting parameters calculated are
// in meters (because the datum used is in meters) and are relative to the
// falseNorthing/falseEasting coordinate. Which in turn is relative to the
// Lat/Long of origin The formula were obtained from URL:
// http://www.ihsenergy.com/epsg/guid7_2.html.
// Code was written by Brenor Brophy, brenor dot brophy at gmail dot com
//
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);
}
//------------------------------------------------------------------------------
//
// Convert Easting/Northing on a Lambert Conic projection to Longitude/Latitude
//
// This routine will convert a Lambert Northing/Easting coordinate to an
// Latitude/Longitude coordinate. The configLambertProjection() function should
// have been called prior to this one to setup the specific parameters for the
// projection. The Northing/Easting parameters are in meters (because the datum
// used is in meters) and are relative to the falseNorthing/falseEasting
// coordinate. Which in turn is relative to the Lat/Long of origin The formula
// were obtained from URL http://www.ihsenergy.com/epsg/guid7_2.html. Code
// was written by Brenor Brophy, brenor dot brophy at gmail dot com
//
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 is a useful function that returns the Great Circle distance from the
// gPoint to another Long/Lat coordinate
//
// Result is returned as meters
//
function distanceFrom(\$lon1, \$lat1)
{

\$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
}

//------------------------------------------------------------------------------
// This function also calculates the distance between two points. In this case
// it just uses Pythagoras's theorm using TM coordinates.
//
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;
}

//------------------------------------------------------------------------------
// This function geo-references a geoPoint to a given map. This means that it
// calculates the x,y pixel coordinate that coresponds to the Lat/Long value of
// the geoPoint. The calculation is done using the Transverse Mercator(TM)
// coordinates of the gPoint with respect to the TM coordinates of the center
// point of the map. So this only makes sense if you are using Local TM
// projection.
//
// \$rX & \$rY are the pixel coordinates that corespond to the Northing/Easting
// (\$rE/\$rN) coordinate it is to this coordinate that the point will be
// geo-referenced. The \$LongOrigin is needed to make sure the Easting/Northing
// coordinates of the point are correctly converted.
//
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

?> About us Advertise on this site Site map Newsletter Statistics Site tips Privacy policy Contact
For more information send a message to `info at phpclasses dot org`.