Subversion Repositories eFlore/Applications.cel

Rev

Rev 3473 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
442 jpm 1
<?php
2
/*------------------------------------------------------------------------------
3
** File:		gPoint.php
4
** Description:	PHP class to convert Latitude & Longitude coordinates into
5
**				UTM & Lambert Conic Conformal Northing/Easting coordinates.
6
** Version:		1.3
7
** Author:		Brenor Brophy
8
** Email:		brenor dot brophy at gmail dot com
9
** Homepage:	brenorbrophy.com
10
**------------------------------------------------------------------------------
11
** COPYRIGHT (c) 2005, 2006, 2007, 2008 BRENOR BROPHY
12
**
13
** The source code included in this package is free software; you can
14
** redistribute it and/or modify it under the terms of the GNU General Public
15
** License as published by the Free Software Foundation. This license can be
16
** read at:
17
**
18
** http://www.opensource.org/licenses/gpl-license.php
19
**
20
** This program is distributed in the hope that it will be useful, but WITHOUT
21
** ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
22
** FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
23
**------------------------------------------------------------------------------
24
**
25
** Code for datum and UTM conversion was converted from C++ code written by
26
** Chuck Gantz (chuck dot gantz at globalstar dot com) from
27
** http://www.gpsy.com/gpsinfo/geotoutm/ This URL has many other references to
28
** useful information concerning conversion of coordinates.
29
**
30
** Rev History
31
** -----------------------------------------------------------------------------
32
** 1.0	08/25/2005	Initial Release
33
** 1.1	05/15/2006	Added software license language to header comments
34
**					Fixed an error in the convertTMtoLL() method. The latitude
35
**					calculation had a bunch of variables without $ symbols.
36
**					Fixed an error in convertLLtoTM() method, The $this-> was
37
**					missing in front of a couple of variables. Thanks to Bob
38
**					Robins of Maryland for catching the bugs.
39
** 1.2	05/18/2007	Added default of NULL to $LongOrigin arguement in convertTMtoLL()
40
**					and convertLLtoTM() to eliminate warning messages when the
41
**					methods are called without a value for $LongOrigin.
42
** 1.3	02/21/2008	Fixed a bug in the distanceFrom method, where the input parameters
43
**					were not being converted to radians prior to calculating the
44
**					distance. Thanks to Enrico Benco for finding pointing it out.
45
*/
46
define ("meter2nm", (1/1852));
47
define ("nm2meter", 1852);
48
 
49
/*------------------------------------------------------------------------------
50
** class gPoint		... for Geographic Point
51
**
52
** This class encapsulates the methods for representing a geographic point on the
53
** earth in three different coordinate systema. Lat/Long, UTM and Lambert Conic
54
** Conformal.
55
*/
56
class gPoint
57
{
58
/* Reference ellipsoids derived from Peter H. Dana's website-
59
** http://www.colorado.edu/geography/gcraft/notes/datum/datum_f.html
60
** email: pdana@pdana.com, web page: www.pdana.com
61
**
62
** Source:
63
** Defense Mapping Agency. 1987b. DMA Technical Report: Supplement to Department
64
** of Defense World Geodetic System 1984 Technical Report. Part I and II.
65
** Washington, DC: Defense Mapping Agency
66
*/
67
	var $ellipsoid = array(//Ellipsoid name, Equatorial Radius, square of eccentricity
68
	"Airy"					=>array (6377563, 0.00667054),
69
	"Australian National"	=>array	(6378160, 0.006694542),
70
	"Bessel 1841"			=>array	(6377397, 0.006674372),
71
	"Bessel 1841 Nambia"	=>array	(6377484, 0.006674372),
72
	"Clarke 1866"			=>array	(6378206, 0.006768658),
73
	"Clarke 1880"			=>array	(6378249, 0.006803511),
74
	"Everest"				=>array	(6377276, 0.006637847),
75
	"Fischer 1960 Mercury"	=>array (6378166, 0.006693422),
76
	"Fischer 1968"			=>array (6378150, 0.006693422),
77
	"GRS 1967"				=>array	(6378160, 0.006694605),
78
	"GRS 1980"				=>array	(6378137, 0.00669438),
79
	"Helmert 1906"			=>array	(6378200, 0.006693422),
80
	"Hough"					=>array	(6378270, 0.00672267),
81
	"International"			=>array	(6378388, 0.00672267),
82
	"Krassovsky"			=>array	(6378245, 0.006693422),
83
	"Modified Airy"			=>array	(6377340, 0.00667054),
84
	"Modified Everest"		=>array	(6377304, 0.006637847),
85
	"Modified Fischer 1960"	=>array	(6378155, 0.006693422),
86
	"South American 1969"	=>array	(6378160, 0.006694542),
87
	"WGS 60"				=>array (6378165, 0.006693422),
88
	"WGS 66"				=>array (6378145, 0.006694542),
89
	"WGS 72"				=>array (6378135, 0.006694318),
90
	"WGS 84"				=>array (6378137, 0.00669438));
91
 
92
	// Properties
93
	var $a;										// Equatorial Radius
94
	var	$e2;									// Square of eccentricity
95
	var	$datum;									// Selected datum
96
	var	$Xp, $Yp;								// X,Y pixel location
97
	var $lat, $long;							// Latitude & Longitude of the point
98
	var $utmNorthing, $utmEasting, $utmZone;	// UTM Coordinates of the point
99
	var	$lccNorthing, $lccEasting;				// Lambert coordinates of the point
100
	var $falseNorthing, $falseEasting;			// Origin coordinates for Lambert Projection
101
	var $latOfOrigin;							// For Lambert Projection
102
	var $longOfOrigin;							// For Lambert Projection
103
	var $firstStdParallel;						// For lambert Projection
104
	var $secondStdParallel;						// For lambert Projection
105
 
106
	// constructor
3473 killian 107
	function __construct($datum='WGS 84')			// Default datum is WGS 84
442 jpm 108
	{
109
		$this->a = $this->ellipsoid[$datum][0];		// Set datum Equatorial Radius
110
		$this->e2 = $this->ellipsoid[$datum][1];	// Set datum Square of eccentricity
111
		$this->datum = $datum;						// Save the datum
112
	}
113
	//
114
	// Set/Get X & Y pixel of the point (used if it is being drawn on an image)
115
	//
116
	function setXY($x, $y)
117
	{
118
		$this->Xp = $x; $this->Yp = $y;
119
	}
120
	function Xp() { return $this->Xp; }
121
	function Yp() { return $this->Yp; }
122
	//
123
	// Set/Get/Output Longitude & Latitude of the point
124
	//
125
	function setLongLat($long, $lat)
126
	{
127
		$this->long = $long; $this->lat = $lat;
128
	}
129
	function Lat() { return $this->lat; }
130
	function Long() { return $this->long; }
131
	function printLatLong() { printf("Latitude: %1.5f Longitude: %1.5f",$this->lat, $this->long); }
132
	//
133
	// Set/Get/Output Universal Transverse Mercator Coordinates
134
	//
135
	function setUTM($easting, $northing, $zone='')	// Zone is optional
136
	{
137
		$this->utmNorthing = $northing;
138
		$this->utmEasting = $easting;
139
		$this->utmZone = $zone;
140
	}
141
	function N() { return $this->utmNorthing; }
142
	function E() { return $this->utmEasting; }
143
	function Z() { return $this->utmZone; }
144
	function printUTM() { print( "Northing: ".(int)$this->utmNorthing.", Easting: ".(int)$this->utmEasting.", Zone: ".$this->utmZone); }
145
	//
146
	// Set/Get/Output Lambert Conic Conformal Coordinates
147
	//
148
	function setLambert($easting, $northing)
149
	{
150
		$this->lccNorthing = $northing;
151
		$this->lccEasting = $easting;
152
	}
153
	function lccN() { return $this->lccNorthing; }
154
	function lccE() { return $this->lccEasting; }
155
	function printLambert() { print( "Northing: ".(int)$this->lccNorthing.", Easting: ".(int)$this->lccEasting); }
156
 
157
//------------------------------------------------------------------------------
158
//
159
// Convert Longitude/Latitude to UTM
160
//
161
// Equations from USGS Bulletin 1532
162
// East Longitudes are positive, West longitudes are negative.
163
// North latitudes are positive, South latitudes are negative
164
// Lat and Long are in decimal degrees
165
// Written by Chuck Gantz- chuck dot gantz at globalstar dot com, converted to PHP by
166
// Brenor Brophy, brenor dot brophy at gmail dot com
167
//
168
// UTM coordinates are useful when dealing with paper maps. Basically the
169
// map will can cover a single UTM zone which is 6 degrees on longitude.
170
// So you really don't care about an object crossing two zones. You just get a
171
// second map of the other zone. However, if you happen to live in a place that
172
// straddles two zones (For example the Santa Babara area in CA straddles zone 10
173
// and zone 11) Then it can become a real pain having to have two maps all the time.
174
// So relatively small parts of the world (like say California) create their own
175
// version of UTM coordinates that are adjusted to conver the whole area of interest
176
// on a single map. These are called state grids. The projection system is the
177
// usually same as UTM (i.e. Transverse Mercator), but the central meridian
178
// aka Longitude of Origin is selected to suit the logitude of the area being
179
// mapped (like being moved to the central meridian of the area) and the grid
180
// may cover more than the 6 degrees of lingitude found on a UTM map. Areas
181
// that are wide rather than long - think Montana as an example. May still
182
// have to have a couple of maps to cover the whole state because TM projection
183
// looses accuracy as you move further away from the Longitude of Origin, 15 degrees
184
// is usually the limit.
185
//
186
// Now, in the case where we want to generate electronic maps that may be
187
// placed pretty much anywhere on the globe we really don't to deal with the
188
// issue of UTM zones in our coordinate system. We would really just like a
189
// grid that is fully contigious over the area of the map we are drawing. Similiar
190
// to the state grid, but local to the area we are interested in. I call this
191
// Local Transverse Mercator and I have modified the function below to also
192
// make this conversion. If you pass a Longitude value to the function as $LongOrigin
193
// then that is the Longitude of Origin that will be used for the projection.
194
// Easting coordinates will be returned (in meters) relative to that line of
195
// longitude - So an Easting coordinate for a point located East of the longitude
196
// of origin will be a positive value in meters, an Easting coordinate for a point
197
// West of the longitude of Origin will have a negative value in meters. Northings
198
// will always be returned in meters from the equator same as the UTM system. The
199
// UTMZone value will be valid for Long/Lat given - thought it is not meaningful
200
// in the context of Local TM. If a NULL value is passed for $LongOrigin
201
// then the standard UTM coordinates are calculated.
202
//
203
	function convertLLtoTM($LongOrigin = NULL)
204
	{
205
		$k0 = 0.9996;
206
		$falseEasting = 0.0;
207
 
208
		//Make sure the longitude is between -180.00 .. 179.9
209
		$LongTemp = ($this->long+180)-(integer)(($this->long+180)/360)*360-180; // -180.00 .. 179.9;
210
		$LatRad = deg2rad($this->lat);
211
		$LongRad = deg2rad($LongTemp);
212
 
213
		if (!$LongOrigin)
214
		{ // Do a standard UTM conversion - so findout what zone the point is in
215
			$ZoneNumber = (integer)(($LongTemp + 180)/6) + 1;
216
			// Special zone for South Norway
217
			if( $this->lat >= 56.0 && $this->lat < 64.0 && $LongTemp >= 3.0 && $LongTemp < 12.0 ) // Fixed 1.1
218
				$ZoneNumber = 32;
219
			// Special zones for Svalbard
220
			if( $this->lat >= 72.0 && $this->lat < 84.0 )
221
			{
222
				if(      $LongTemp >= 0.0  && $LongTemp <  9.0 ) $ZoneNumber = 31;
223
				else if( $LongTemp >= 9.0  && $LongTemp < 21.0 ) $ZoneNumber = 33;
224
				else if( $LongTemp >= 21.0 && $LongTemp < 33.0 ) $ZoneNumber = 35;
225
				else if( $LongTemp >= 33.0 && $LongTemp < 42.0 ) $ZoneNumber = 37;
226
			}
227
			$LongOrigin = ($ZoneNumber - 1)*6 - 180 + 3;  //+3 puts origin in middle of zone
228
			//compute the UTM Zone from the latitude and longitude
229
			$this->utmZone = sprintf("%d%s", $ZoneNumber, $this->UTMLetterDesignator());
230
			// We also need to set the false Easting value adjust the UTM easting coordinate
231
			$falseEasting = 500000.0;
232
		}
233
		$LongOriginRad = deg2rad($LongOrigin);
234
 
235
		$eccPrimeSquared = ($this->e2)/(1-$this->e2);
236
 
237
		$N = $this->a/sqrt(1-$this->e2*sin($LatRad)*sin($LatRad));
238
		$T = tan($LatRad)*tan($LatRad);
239
		$C = $eccPrimeSquared*cos($LatRad)*cos($LatRad);
240
		$A = cos($LatRad)*($LongRad-$LongOriginRad);
241
 
242
		$M = $this->a*((1	- $this->e2/4		- 3*$this->e2*$this->e2/64	- 5*$this->e2*$this->e2*$this->e2/256)*$LatRad
243
							- (3*$this->e2/8	+ 3*$this->e2*$this->e2/32	+ 45*$this->e2*$this->e2*$this->e2/1024)*sin(2*$LatRad)
244
												+ (15*$this->e2*$this->e2/256 + 45*$this->e2*$this->e2*$this->e2/1024)*sin(4*$LatRad)
245
												- (35*$this->e2*$this->e2*$this->e2/3072)*sin(6*$LatRad));
246
 
247
		$this->utmEasting = ($k0*$N*($A+(1-$T+$C)*$A*$A*$A/6
248
						+ (5-18*$T+$T*$T+72*$C-58*$eccPrimeSquared)*$A*$A*$A*$A*$A/120)
249
						+ $falseEasting);
250
 
251
		$this->utmNorthing = ($k0*($M+$N*tan($LatRad)*($A*$A/2+(5-$T+9*$C+4*$C*$C)*$A*$A*$A*$A/24
252
					 + (61-58*$T+$T*$T+600*$C-330*$eccPrimeSquared)*$A*$A*$A*$A*$A*$A/720)));
253
		if($this->lat < 0)
254
			$this->utmNorthing += 10000000.0; //10000000 meter offset for southern hemisphere
255
	}
256
//
257
// This routine determines the correct UTM letter designator for the given latitude
258
// returns 'Z' if latitude is outside the UTM limits of 84N to 80S
259
// Written by Chuck Gantz- chuck dot gantz at globalstar dot com, converted to PHP by
260
// Brenor Brophy, brenor dot brophy at gmail dot com
261
//
262
	function UTMLetterDesignator()
263
	{
264
		if((84 >= $this->lat) && ($this->lat >= 72)) $LetterDesignator = 'X';
265
		else if((72 > $this->lat) && ($this->lat >= 64)) $LetterDesignator = 'W';
266
		else if((64 > $this->lat) && ($this->lat >= 56)) $LetterDesignator = 'V';
267
		else if((56 > $this->lat) && ($this->lat >= 48)) $LetterDesignator = 'U';
268
		else if((48 > $this->lat) && ($this->lat >= 40)) $LetterDesignator = 'T';
269
		else if((40 > $this->lat) && ($this->lat >= 32)) $LetterDesignator = 'S';
270
		else if((32 > $this->lat) && ($this->lat >= 24)) $LetterDesignator = 'R';
271
		else if((24 > $this->lat) && ($this->lat >= 16)) $LetterDesignator = 'Q';
272
		else if((16 > $this->lat) && ($this->lat >= 8)) $LetterDesignator = 'P';
273
		else if(( 8 > $this->lat) && ($this->lat >= 0)) $LetterDesignator = 'N';
274
		else if(( 0 > $this->lat) && ($this->lat >= -8)) $LetterDesignator = 'M';
275
		else if((-8 > $this->lat) && ($this->lat >= -16)) $LetterDesignator = 'L';
276
		else if((-16 > $this->lat) && ($this->lat >= -24)) $LetterDesignator = 'K';
277
		else if((-24 > $this->lat) && ($this->lat >= -32)) $LetterDesignator = 'J';
278
		else if((-32 > $this->lat) && ($this->lat >= -40)) $LetterDesignator = 'H';
279
		else if((-40 > $this->lat) && ($this->lat >= -48)) $LetterDesignator = 'G';
280
		else if((-48 > $this->lat) && ($this->lat >= -56)) $LetterDesignator = 'F';
281
		else if((-56 > $this->lat) && ($this->lat >= -64)) $LetterDesignator = 'E';
282
		else if((-64 > $this->lat) && ($this->lat >= -72)) $LetterDesignator = 'D';
283
		else if((-72 > $this->lat) && ($this->lat >= -80)) $LetterDesignator = 'C';
284
		else $LetterDesignator = 'Z'; //This is here as an error flag to show that the Latitude is outside the UTM limits
285
 
286
		return($LetterDesignator);
287
	}
288
 
289
//------------------------------------------------------------------------------
290
//
291
// Convert UTM to Longitude/Latitude
292
//
293
// Equations from USGS Bulletin 1532
294
// East Longitudes are positive, West longitudes are negative.
295
// North latitudes are positive, South latitudes are negative
296
// Lat and Long are in decimal degrees.
297
// Written by Chuck Gantz- chuck dot gantz at globalstar dot com, converted to PHP by
298
// Brenor Brophy, brenor dot brophy at gmail dot com
299
//
300
// If a value is passed for $LongOrigin then the function assumes that
301
// a Local (to the Longitude of Origin passed in) Transverse Mercator
302
// coordinates is to be converted - not a UTM coordinate. This is the
303
// complementary function to the previous one. The function cannot
304
// tell if a set of Northing/Easting coordinates are in the North
305
// or South hemesphere - they just give distance from the equator not
306
// direction - so only northern hemesphere lat/long coordinates are returned.
307
// If you live south of the equator there is a note later in the code
308
// explaining how to have it just return southern hemesphere lat/longs.
309
//
310
	function convertTMtoLL($LongOrigin = NULL)
311
	{
312
		$k0 = 0.9996;
313
		$e1 = (1-sqrt(1-$this->e2))/(1+sqrt(1-$this->e2));
314
		$falseEasting = 0.0;
315
		$y = $this->utmNorthing;
316
 
317
		if (!$LongOrigin)
318
		{ // It is a UTM coordinate we want to convert
319
			sscanf($this->utmZone,"%d%s",$ZoneNumber,$ZoneLetter);
320
			if($ZoneLetter >= 'N')
321
				$NorthernHemisphere = 1;//point is in northern hemisphere
322
			else
323
			{
324
				$NorthernHemisphere = 0;//point is in southern hemisphere
325
				$y -= 10000000.0;//remove 10,000,000 meter offset used for southern hemisphere
326
			}
327
			$LongOrigin = ($ZoneNumber - 1)*6 - 180 + 3;  //+3 puts origin in middle of zone
328
			$falseEasting = 500000.0;
329
		}
330
 
331
//		$y -= 10000000.0;	// Uncomment line to make LOCAL coordinates return southern hemesphere Lat/Long
332
		$x = $this->utmEasting - $falseEasting; //remove 500,000 meter offset for longitude
333
 
334
		$eccPrimeSquared = ($this->e2)/(1-$this->e2);
335
 
336
		$M = $y / $k0;
337
		$mu = $M/($this->a*(1-$this->e2/4-3*$this->e2*$this->e2/64-5*$this->e2*$this->e2*$this->e2/256));
338
 
339
		$phi1Rad = $mu	+ (3*$e1/2-27*$e1*$e1*$e1/32)*sin(2*$mu)
340
					+ (21*$e1*$e1/16-55*$e1*$e1*$e1*$e1/32)*sin(4*$mu)
341
					+(151*$e1*$e1*$e1/96)*sin(6*$mu);
342
		$phi1 = rad2deg($phi1Rad);
343
 
344
		$N1 = $this->a/sqrt(1-$this->e2*sin($phi1Rad)*sin($phi1Rad));
345
		$T1 = tan($phi1Rad)*tan($phi1Rad);
346
		$C1 = $eccPrimeSquared*cos($phi1Rad)*cos($phi1Rad);
347
		$R1 = $this->a*(1-$this->e2)/pow(1-$this->e2*sin($phi1Rad)*sin($phi1Rad), 1.5);
348
		$D = $x/($N1*$k0);
349
 
350
		$tlat = $phi1Rad - ($N1*tan($phi1Rad)/$R1)*($D*$D/2-(5+3*$T1+10*$C1-4*$C1*$C1-9*$eccPrimeSquared)*$D*$D*$D*$D/24
351
						+(61+90*$T1+298*$C1+45*$T1*$T1-252*$eccPrimeSquared-3*$C1*$C1)*$D*$D*$D*$D*$D*$D/720); // fixed in 1.1
352
		$this->lat = rad2deg($tlat);
353
 
354
		$tlong = ($D-(1+2*$T1+$C1)*$D*$D*$D/6+(5-2*$C1+28*$T1-3*$C1*$C1+8*$eccPrimeSquared+24*$T1*$T1)
355
						*$D*$D*$D*$D*$D/120)/cos($phi1Rad);
356
		$this->long = $LongOrigin + rad2deg($tlong);
357
	}
358
 
359
//------------------------------------------------------------------------------
360
// Configure a Lambert Conic Conformal Projection
361
//
362
// falseEasting & falseNorthing are just an offset in meters added to the final
363
// coordinate calculated.
364
//
365
// longOfOrigin & LatOfOrigin are the "center" latitiude and longitude of the
366
// area being projected. All coordinates will be calculated in meters relative
367
// to this point on the earth.
368
//
369
// firstStdParallel & secondStdParallel are the two lines of longitude (that
370
// is they run east-west) that define where the "cone" intersects the earth.
371
// Simply put they should bracket the area being projected.
372
//
373
// google is your friend to find out more
374
//
375
	function configLambertProjection ($falseEasting, $falseNorthing,
376
									  $longOfOrigin, $latOfOrigin,
377
									  $firstStdParallel, $secondStdParallel)
378
	{
379
		$this->falseEasting = $falseEasting;
380
		$this->falseNorthing = $falseNorthing;
381
		$this->longOfOrigin = $longOfOrigin;
382
		$this->latOfOrigin = $latOfOrigin;
383
		$this->firstStdParallel = $firstStdParallel;
384
		$this->secondStdParallel = $secondStdParallel;
385
	}
386
 
387
//------------------------------------------------------------------------------
388
//
389
// Convert Longitude/Latitude to Lambert Conic Easting/Northing
390
//
391
// This routine will convert a Latitude/Longitude coordinate to an Northing/
392
// Easting coordinate on a Lambert Conic Projection. The configLambertProjection()
393
// function should have been called prior to this one to setup the specific
394
// parameters for the projection. The Northing/Easting parameters calculated are
395
// in meters (because the datum used is in meters) and are relative to the
396
// falseNorthing/falseEasting coordinate. Which in turn is relative to the
397
// Lat/Long of origin The formula were obtained from URL:
398
// http://www.ihsenergy.com/epsg/guid7_2.html.
399
// Code was written by Brenor Brophy, brenor dot brophy at gmail dot com
400
//
401
	function convertLLtoLCC()
402
	{
403
		$e = sqrt($this->e2);
404
 
405
		$phi 	= deg2rad($this->lat);						// Latitude to convert
406
		$phi1	= deg2rad($this->firstStdParallel);			// Latitude of 1st std parallel
407
		$phi2	= deg2rad($this->secondStdParallel);		// Latitude of 2nd std parallel
408
		$lamda	= deg2rad($this->long);						// Lonitude to convert
409
		$phio	= deg2rad($this->latOfOrigin);				// Latitude of  Origin
410
		$lamdao	= deg2rad($this->longOfOrigin);				// Longitude of  Origin
411
 
412
		$m1 = cos($phi1) / sqrt(( 1 - $this->e2*sin($phi1)*sin($phi1)));
413
		$m2 = cos($phi2) / sqrt(( 1 - $this->e2*sin($phi2)*sin($phi2)));
414
		$t1 = tan((pi()/4)-($phi1/2)) / pow(( ( 1 - $e*sin($phi1) ) / ( 1 + $e*sin($phi1) )),$e/2);
415
		$t2 = tan((pi()/4)-($phi2/2)) / pow(( ( 1 - $e*sin($phi2) ) / ( 1 + $e*sin($phi2) )),$e/2);
416
		$to = tan((pi()/4)-($phio/2)) / pow(( ( 1 - $e*sin($phio) ) / ( 1 + $e*sin($phio) )),$e/2);
417
		$t  = tan((pi()/4)-($phi /2)) / pow(( ( 1 - $e*sin($phi ) ) / ( 1 + $e*sin($phi ) )),$e/2);
418
		$n	= (log($m1)-log($m2)) / (log($t1)-log($t2));
419
		$F	= $m1/($n*pow($t1,$n));
420
		$rf	= $this->a*$F*pow($to,$n);
421
		$r	= $this->a*$F*pow($t,$n);
422
		$theta = $n*($lamda - $lamdao);
423
 
424
		$this->lccEasting = $this->falseEasting + $r*sin($theta);
425
		$this->lccNorthing = $this->falseNorthing + $rf - $r*cos($theta);
426
	}
427
//------------------------------------------------------------------------------
428
//
429
// Convert Easting/Northing on a Lambert Conic projection to Longitude/Latitude
430
//
431
// This routine will convert a Lambert Northing/Easting coordinate to an
432
// Latitude/Longitude coordinate.  The configLambertProjection() function should
433
// have been called prior to this one to setup the specific parameters for the
434
// projection. The Northing/Easting parameters are in meters (because the datum
435
// used is in meters) and are relative to the falseNorthing/falseEasting
436
// coordinate. Which in turn is relative to the Lat/Long of origin The formula
437
// were obtained from URL http://www.ihsenergy.com/epsg/guid7_2.html. Code
438
// was written by Brenor Brophy, brenor dot brophy at gmail dot com
439
//
440
	function convertLCCtoLL()
441
	{
442
		$e = sqrt($this->e2);
443
 
444
		$phi1	= deg2rad($this->firstStdParallel);			// Latitude of 1st std parallel
445
		$phi2	= deg2rad($this->secondStdParallel);		// Latitude of 2nd std parallel
446
		$phio	= deg2rad($this->latOfOrigin);				// Latitude of  Origin
447
		$lamdao	= deg2rad($this->longOfOrigin);				// Longitude of  Origin
448
		$E		= $this->lccEasting;
449
		$N		= $this->lccNorthing;
450
		$Ef		= $this->falseEasting;
451
		$Nf		= $this->falseNorthing;
452
 
453
		$m1 = cos($phi1) / sqrt(( 1 - $this->e2*sin($phi1)*sin($phi1)));
454
		$m2 = cos($phi2) / sqrt(( 1 - $this->e2*sin($phi2)*sin($phi2)));
455
		$t1 = tan((pi()/4)-($phi1/2)) / pow(( ( 1 - $e*sin($phi1) ) / ( 1 + $e*sin($phi1) )),$e/2);
456
		$t2 = tan((pi()/4)-($phi2/2)) / pow(( ( 1 - $e*sin($phi2) ) / ( 1 + $e*sin($phi2) )),$e/2);
457
		$to = tan((pi()/4)-($phio/2)) / pow(( ( 1 - $e*sin($phio) ) / ( 1 + $e*sin($phio) )),$e/2);
458
		$n	= (log($m1)-log($m2)) / (log($t1)-log($t2));
459
		$F	= $m1/($n*pow($t1,$n));
460
		$rf	= $this->a*$F*pow($to,$n);
461
		$r_	= sqrt( pow(($E-$Ef),2) + pow(($rf-($N-$Nf)),2) );
462
		$t_	= pow($r_/($this->a*$F),(1/$n));
463
		$theta_ = atan(($E-$Ef)/($rf-($N-$Nf)));
464
 
465
		$lamda	= $theta_/$n + $lamdao;
466
		$phi0	= (pi()/2) - 2*atan($t_);
467
		$phi1	= (pi()/2) - 2*atan($t_*pow(((1-$e*sin($phi0))/(1+$e*sin($phi0))),$e/2));
468
		$phi2	= (pi()/2) - 2*atan($t_*pow(((1-$e*sin($phi1))/(1+$e*sin($phi1))),$e/2));
469
		$phi	= (pi()/2) - 2*atan($t_*pow(((1-$e*sin($phi2))/(1+$e*sin($phi2))),$e/2));
470
 
471
		$this->lat 	= rad2deg($phi);
472
		$this->long = rad2deg($lamda);
473
	}
474
 
475
//------------------------------------------------------------------------------
476
// This is a useful function that returns the Great Circle distance from the
477
// gPoint to another Long/Lat coordinate
478
//
479
// Result is returned as meters
480
//
481
	function distanceFrom($lon1, $lat1)
482
	{
483
		$lon1 = deg2rad($lon1); $lat1 = deg2rad($lat1); // Added in 1.3
484
		$lon2 = deg2rad($this->Long()); $lat2 = deg2rad($this->Lat());
485
 
486
		$theta = $lon2 - $lon1;
487
		$dist = acos(sin($lat1) * sin($lat2) + cos($lat1) * cos($lat2) * cos($theta));
488
 
489
//		Alternative formula supposed to be more accurate for short distances
490
//		$dist = 2*asin(sqrt( pow(sin(($lat1-$lat2)/2),2) + cos($lat1)*cos($lat2)*pow(sin(($lon1-$lon2)/2),2)));
491
		return ( $dist * 6366710 ); // from http://williams.best.vwh.net/avform.htm#GCF
492
	}
493
 
494
//------------------------------------------------------------------------------
495
// This function also calculates the distance between two points. In this case
496
// it just uses Pythagoras's theorm using TM coordinates.
497
//
498
	function distanceFromTM(&$pt)
499
	{
500
		$E1 = $pt->E(); 	$N1 = $pt->N();
501
		$E2 = $this->E(); 	$N2 = $this->N();
502
 
503
		$dist = sqrt(pow(($E1-$E2),2)+pow(($N1-$N2),2));
504
		return $dist;
505
	}
506
 
507
//------------------------------------------------------------------------------
508
// This function geo-references a geoPoint to a given map. This means that it
509
// calculates the x,y pixel coordinate that coresponds to the Lat/Long value of
510
// the geoPoint. The calculation is done using the Transverse Mercator(TM)
511
// coordinates of the gPoint with respect to the TM coordinates of the center
512
// point of the map. So this only makes sense if you are using Local TM
513
// projection.
514
//
515
// $rX & $rY are the pixel coordinates that corespond to the Northing/Easting
516
// ($rE/$rN) coordinate it is to this coordinate that the point will be
517
// geo-referenced. The $LongOrigin is needed to make sure the Easting/Northing
518
// coordinates of the point are correctly converted.
519
//
520
	function gRef($rX, $rY, $rE, $rN, $Scale, $LongOrigin)
521
	{
522
		$this->convertLLtoTM($LongOrigin);
523
		$x = (($this->E() - $rE) / $Scale)		// The easting in meters times the scale to get pixels
524
												// is relative to the center of the image so adjust to
525
			+ ($rX);							// the left coordinate.
526
		$y = $rY -  							// Adjust to bottom coordinate.
527
			(($rN - $this->N()) / $Scale);		// The northing in meters
528
												// relative to the equator. Subtract center point northing
529
												// to get relative to image center and convert meters to pixels
530
		$this->setXY((int)$x,(int)$y);			// Save the geo-referenced result.
531
	}
532
} // end of class gPoint
533
 
534
?>