Details | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
365 | rain-er | 1 | # Geo::Ellipsoid |
2 | # |
||
3 | # This package implements an Ellipsoid class to perform latitude |
||
4 | # and longitude calculations on the surface of an ellipsoid. |
||
5 | # |
||
6 | # This is a Perl conversion of existing Fortran code (see |
||
7 | # ACKNOWLEDGEMENTS) and the author of this class makes no |
||
8 | # claims of originality. Nor can he even vouch for the |
||
9 | # results of the calculations, although they do seem to |
||
10 | # work for him and have been tested against other methods. |
||
11 | |||
12 | package Geo::Ellipsoid; |
||
13 | |||
14 | use warnings; |
||
15 | use strict; |
||
16 | use 5.006_00; |
||
17 | |||
18 | use Scalar::Util 'looks_like_number'; |
||
19 | use Math::Trig; |
||
20 | use Carp; |
||
21 | |||
22 | =head1 NAME |
||
23 | |||
24 | Geo::Ellipsoid - Longitude and latitude calculations using an ellipsoid model. |
||
25 | |||
26 | =head1 VERSION |
||
27 | |||
28 | Version 1.12, released July 4, 2008. |
||
29 | |||
30 | =cut |
||
31 | |||
32 | our $VERSION = '1.12'; |
||
33 | our $DEBUG = 0; |
||
34 | |||
35 | =head1 SYNOPSIS |
||
36 | |||
37 | use Geo::Ellipsoid; |
||
38 | $geo = Geo::Ellipsoid->new(ellipsoid=>'NAD27', units=>'degrees'); |
||
39 | @origin = ( 37.619002, -122.374843 ); # SFO |
||
40 | @dest = ( 33.942536, -118.408074 ); # LAX |
||
41 | ( $range, $bearing ) = $geo->to( @origin, @dest ); |
||
42 | ($lat,$lon) = $geo->at( @origin, 2000, 45.0 ); |
||
43 | ( $x, $y ) = $geo->displacement( @origin, $lat, $lon ); |
||
44 | @pos = $geo->location( $lat, $lon, $x, $y ); |
||
45 | |||
46 | =head1 DESCRIPTION |
||
47 | |||
48 | Geo::Ellipsoid performs geometrical calculations on the surface of |
||
49 | an ellipsoid. An ellipsoid is a three-dimension object formed from |
||
50 | the rotation of an ellipse about one of its axes. The approximate |
||
51 | shape of the earth is an ellipsoid, so Geo::Ellipsoid can accurately |
||
52 | calculate distance and bearing between two widely-separated locations |
||
53 | on the earth's surface. |
||
54 | |||
55 | The shape of an ellipsoid is defined by the lengths of its |
||
56 | semi-major and semi-minor axes. The shape may also be specifed by |
||
57 | the flattening ratio C<f> as: |
||
58 | |||
59 | f = ( semi-major - semi-minor ) / semi-major |
||
60 | |||
61 | which, since f is a small number, is normally given as the reciprocal |
||
62 | of the flattening C<1/f>. |
||
63 | |||
64 | The shape of the earth has been surveyed and estimated differently |
||
65 | at different times over the years. The two most common sets of values |
||
66 | used to describe the size and shape of the earth in the United States |
||
67 | are 'NAD27', dating from 1927, and 'WGS84', from 1984. United States |
||
68 | Geological Survey topographical maps, for example, use one or the |
||
69 | other of these values, and commonly-available Global Positioning |
||
70 | System (GPS) units can be set to use one or the other. |
||
71 | See L<"DEFINED ELLIPSOIDS"> below for the ellipsoid survey values |
||
72 | that may be selected for use by Geo::Ellipsoid. |
||
73 | |||
74 | =cut |
||
75 | |||
76 | # class data and constants |
||
77 | our $degrees_per_radian = 180/pi; |
||
78 | our $eps = 1.0e-23; |
||
79 | our $max_loop_count = 20; |
||
80 | our $twopi = 2 * pi; |
||
81 | our $halfpi = pi/2; |
||
82 | our %defaults = ( |
||
83 | ellipsoid => 'WGS84', |
||
84 | units => 'radians', |
||
85 | distance_units => 'meter', |
||
86 | longitude => 0, |
||
87 | latitude => 1, # allows use of _normalize_output |
||
88 | bearing => 0, |
||
89 | ); |
||
90 | our %distance = ( |
||
91 | 'foot' => 0.3048, |
||
92 | 'kilometer' => 1_000, |
||
93 | 'meter' => 1.0, |
||
94 | 'mile' => 1_609.344, |
||
95 | 'nm' => 1_852, |
||
96 | ); |
||
97 | |||
98 | # set of ellipsoids that can be used. |
||
99 | # values are |
||
100 | # 1) a = semi-major (equatorial) radius of Ellipsoid |
||
101 | # 2) 1/f = reciprocal of flattening (f), the ratio of the semi-minor |
||
102 | # (polar) radius to the semi-major (equatorial) axis, or |
||
103 | # polar radius = equatorial radius * ( 1 - f ) |
||
104 | |||
105 | our %ellipsoids = ( |
||
106 | 'AIRY' => [ 6377563.396, 299.3249646 ], |
||
107 | 'AIRY-MODIFIED' => [ 6377340.189, 299.3249646 ], |
||
108 | 'AUSTRALIAN' => [ 6378160.0, 298.25 ], |
||
109 | 'BESSEL-1841' => [ 6377397.155, 299.1528128 ], |
||
110 | 'CLARKE-1880' => [ 6378249.145, 293.465 ], |
||
111 | 'EVEREST-1830' => [ 6377276.345, 300.8017 ], |
||
112 | 'EVEREST-MODIFIED' => [ 6377304.063, 300.8017 ], |
||
113 | 'FISHER-1960' => [ 6378166.0, 298.3 ], |
||
114 | 'FISHER-1968' => [ 6378150.0, 298.3 ], |
||
115 | 'GRS80' => [ 6378137.0, 298.25722210088 ], |
||
116 | 'HOUGH-1956' => [ 6378270.0, 297.0 ], |
||
117 | 'HAYFORD' => [ 6378388.0, 297.0 ], |
||
118 | 'IAU76' => [ 6378140.0, 298.257 ], |
||
119 | 'KRASSOVSKY-1938' => [ 6378245.0, 298.3 ], |
||
120 | 'NAD27' => [ 6378206.4, 294.9786982138 ], |
||
121 | 'NWL-9D' => [ 6378145.0, 298.25 ], |
||
122 | 'SOUTHAMERICAN-1969' => [ 6378160.0, 298.25 ], |
||
123 | 'SOVIET-1985' => [ 6378136.0, 298.257 ], |
||
124 | 'WGS72' => [ 6378135.0, 298.26 ], |
||
125 | 'WGS84' => [ 6378137.0, 298.257223563 ], |
||
126 | ); |
||
127 | |||
128 | =head1 CONSTRUCTOR |
||
129 | |||
130 | =head2 new |
||
131 | |||
132 | The new() constructor may be called with a hash list to set the value of the |
||
133 | ellipsoid to be used, the value of the units to be used for angles and |
||
134 | distances, and whether or not the output range of longitudes and bearing |
||
135 | angles should be symmetric around zero or always greater than zero. |
||
136 | The initial default constructor is equivalent to the following: |
||
137 | |||
138 | my $geo = Geo::Ellipsoid->new( |
||
139 | ellipsoid => 'WGS84', |
||
140 | units => 'radians' , |
||
141 | distance_units => 'meter', |
||
142 | longitude => 0, |
||
143 | bearing => 0, |
||
144 | ); |
||
145 | |||
146 | The constructor arguments may be of any case and, with the exception of |
||
147 | the ellipsoid value, abbreviated to their first three characters. |
||
148 | Thus, ( UNI => 'DEG', DIS => 'FEE', Lon => 1, ell => 'NAD27', bEA => 0 ) |
||
149 | is valid. |
||
150 | |||
151 | =cut |
||
152 | |||
153 | sub new |
||
154 | { |
||
155 | my( $class, %args ) = @_; |
||
156 | my $self = {%defaults}; |
||
157 | print "new: @_\n" if $DEBUG; |
||
158 | foreach my $key ( keys %args ) { |
||
159 | my $val = $args{$key}; |
||
160 | if( $key =~ /^ell/i ) { |
||
161 | $self->{ellipsoid} = uc $args{$key}; |
||
162 | }elsif( $key =~ /^uni/i ) { |
||
163 | $self->{units} = $args{$key}; |
||
164 | }elsif( $key =~ /^dis/i ) { |
||
165 | $self->{distance_units} = $args{$key}; |
||
166 | }elsif( $key =~ /^lon/i ) { |
||
167 | $self->{longitude} = $args{$key}; |
||
168 | }elsif( $key =~ /^bea/i ) { |
||
169 | $self->{bearing} = $args{$key}; |
||
170 | }else{ |
||
171 | carp("Unknown argument: $key => $args{$key}"); |
||
172 | } |
||
173 | } |
||
174 | set_units($self,$self->{units}); |
||
175 | set_ellipsoid($self,$self->{ellipsoid}); |
||
176 | set_distance_unit($self,$self->{distance_units}); |
||
177 | set_longitude_symmetric($self,$self->{longitude}); |
||
178 | set_bearing_symmetric($self,$self->{bearing}); |
||
179 | |||
180 | "Ellipsoid(units=>$self->{units},distance_units=>" . |
||
181 | "$self->{distance_units},ellipsoid=>$self->{ellipsoid}," . |
||
182 | "longitude=>$self->{longitude},bearing=>$self->{bearing})\n" if $DEBUG; |
||
183 | bless $self,$class; |
||
184 | return $self; |
||
185 | } |
||
186 | |||
187 | =head1 METHODS |
||
188 | |||
189 | =head2 set_units |
||
190 | |||
191 | Set the angle units used by the Geo::Ellipsoid object. The units may |
||
192 | also be set in the constructor of the object. The allowable values are |
||
193 | 'degrees' or 'radians'. The default is 'radians'. The units value is |
||
194 | not case sensitive and may be abbreviated to 3 letters. The units of |
||
195 | angle apply to both input and output latitude, longitude, and bearing |
||
196 | values. |
||
197 | |||
198 | $geo->set_units('degrees'); |
||
199 | |||
200 | =cut |
||
201 | |||
202 | sub set_units |
||
203 | { |
||
204 | my $self = shift; |
||
205 | my $units = shift; |
||
206 | if( $units =~ /deg/i ) { |
||
207 | $units = 'degrees'; |
||
208 | }elsif( $units =~ /rad/i ) { |
||
209 | $units = 'radians'; |
||
210 | }else{ |
||
211 | croak("Invalid units specifier '$units' - please use either " . |
||
212 | "degrees or radians (the default)") unless $units =~ /rad/i; |
||
213 | } |
||
214 | $self->{units} = $units; |
||
215 | } |
||
216 | |||
217 | =head2 set_distance_unit |
||
218 | |||
219 | Set the distance unit used by the Geo::Ellipsoid object. The unit of |
||
220 | distance may also be set in the constructor of the object. The recognized |
||
221 | values are 'meter', 'kilometer', 'mile', 'nm' (nautical mile), or 'foot'. |
||
222 | The default is 'meter'. The value is not case sensitive and may be |
||
223 | abbreviated to 3 letters. |
||
224 | |||
225 | $geo->set_distance_unit('kilometer'); |
||
226 | |||
227 | For any other unit of distance not recogized by this method, pass a |
||
228 | numerical argument representing the length of the distance unit in |
||
229 | meters. For example, to use units of furlongs, call |
||
230 | |||
231 | $geo->set_distance_unit(201.168); |
||
232 | |||
233 | The distance conversion factors used by this module are as follows: |
||
234 | |||
235 | Unit Units per meter |
||
236 | -------- --------------- |
||
237 | foot 0.3048 |
||
238 | kilometer 1000.0 |
||
239 | mile 1609.344 |
||
240 | nm 1852.0 |
||
241 | |||
242 | =cut |
||
243 | |||
244 | sub set_distance_unit |
||
245 | { |
||
246 | my $self = shift; |
||
247 | my $unit = shift; |
||
248 | print "distance unit = $unit\n" if $DEBUG; |
||
249 | |||
250 | my $conversion = 0; |
||
251 | |||
252 | if( defined $unit ) { |
||
253 | |||
254 | my( $key, $val ); |
||
255 | while( ($key,$val) = each %distance ) { |
||
256 | my $re = substr($key,0,3); |
||
257 | print "trying ($key,$re,$val)\n" if $DEBUG; |
||
258 | if( $unit =~ /^$re/i ) { |
||
259 | $self->{distance_units} = $unit; |
||
260 | $conversion = $val; |
||
261 | |||
262 | # finish iterating to reset 'each' function call |
||
263 | while( each %distance ) {} |
||
264 | last; |
||
265 | } |
||
266 | } |
||
267 | |||
268 | if( $conversion == 0 ) { |
||
269 | if( looks_like_number($unit) ) { |
||
270 | $conversion = $unit; |
||
271 | }else{ |
||
272 | carp("Unknown argument to set_distance_unit: $unit\nAssuming meters"); |
||
273 | $conversion = 1.0; |
||
274 | } |
||
275 | } |
||
276 | }else{ |
||
277 | carp("Missing or undefined argument to set_distance_unit: ". |
||
278 | "$unit\nAssuming meters"); |
||
279 | $conversion = 1.0; |
||
280 | } |
||
281 | $self->{conversion} = $conversion; |
||
282 | } |
||
283 | |||
284 | =head2 set_ellipsoid |
||
285 | |||
286 | Set the ellipsoid to be used by the Geo::Ellipsoid object. See |
||
287 | L<"DEFINED ELLIPSOIDS"> below for the allowable values. The value |
||
288 | may also be set by the constructor. The default value is 'WGS84'. |
||
289 | |||
290 | $geo->set_ellipsoid('NAD27'); |
||
291 | |||
292 | =cut |
||
293 | |||
294 | sub set_ellipsoid |
||
295 | { |
||
296 | my $self = shift; |
||
297 | my $ellipsoid = uc shift || $defaults{ellipsoid}; |
||
298 | print " set ellipsoid to $ellipsoid\n" if $DEBUG; |
||
299 | unless( exists $ellipsoids{$ellipsoid} ) { |
||
300 | croak("Ellipsoid $ellipsoid does not exist - please use " . |
||
301 | "set_custom_ellipsoid to use an ellipsoid not in valid set"); |
||
302 | } |
||
303 | $self->{ellipsoid} = $ellipsoid; |
||
304 | my( $major, $recip ) = @{$ellipsoids{$ellipsoid}}; |
||
305 | $self->{equatorial} = $major; |
||
306 | if( $recip == 0 ) { |
||
307 | carp("Infinite flattening specified by ellipsoid -- assuming a sphere"); |
||
308 | $self->{polar} = $self->{equatorial}; |
||
309 | $self->{flattening} = 0; |
||
310 | $self->{eccentricity} = 0; |
||
311 | }else{ |
||
312 | $self->{flattening} = ( 1.0 / $ellipsoids{$ellipsoid}[1]); |
||
313 | $self->{polar} = $self->{equatorial} * ( 1.0 - $self->{flattening} ); |
||
314 | $self->{eccentricity} = sqrt( 2.0 * $self->{flattening} - |
||
315 | ( $self->{flattening} * $self->{flattening} ) ); |
||
316 | } |
||
317 | } |
||
318 | |||
319 | =head2 set_custom_ellipsoid |
||
320 | |||
321 | Sets the ellipsoid parameters to the specified ( major semiaxis and |
||
322 | reciprocal flattening. A zero value for the reciprocal flattening |
||
323 | will result in a sphere for the ellipsoid, and a warning message |
||
324 | will be issued. |
||
325 | |||
326 | $geo->set_custom_ellipsoid( 'sphere', 6378137, 0 ); |
||
327 | |||
328 | =cut |
||
329 | |||
330 | sub set_custom_ellipsoid |
||
331 | { |
||
332 | my $self = shift; |
||
333 | my( $name, $major, $recip ) = @_; |
||
334 | $name = uc $name; |
||
335 | $recip = 0 unless defined $recip; |
||
336 | if( $major ) { |
||
337 | $ellipsoids{$name} = [ $major, $recip ]; |
||
338 | }else{ |
||
339 | croak("set_custom_ellipsoid called without semi-major radius parameter"); |
||
340 | } |
||
341 | set_ellipsoid($self,$name); |
||
342 | } |
||
343 | |||
344 | =head2 set_longitude_symmetric |
||
345 | |||
346 | If called with no argument or a true argument, sets the range of output |
||
347 | values for longitude to be in the range [-pi,+pi) radians. If called with |
||
348 | a false or undefined argument, sets the output angle range to be |
||
349 | [0,2*pi) radians. |
||
350 | |||
351 | $geo->set_longitude_symmetric(1); |
||
352 | |||
353 | =cut |
||
354 | |||
355 | sub set_longitude_symmetric |
||
356 | { |
||
357 | my( $self, $sym ) = @_; |
||
358 | # see if argument passed |
||
359 | if( $#_ > 0 ) { |
||
360 | # yes -- use value passed |
||
361 | $self->{longitude} = $sym; |
||
362 | }else{ |
||
363 | # no -- set to true |
||
364 | $self->{longitude} = 1; |
||
365 | } |
||
366 | } |
||
367 | |||
368 | =head2 set_bearing_symmetric |
||
369 | |||
370 | If called with no argument or a true argument, sets the range of output |
||
371 | values for bearing to be in the range [-pi,+pi) radians. If called with |
||
372 | a false or undefined argument, sets the output angle range to be |
||
373 | [0,2*pi) radians. |
||
374 | |||
375 | $geo->set_bearing_symmetric(1); |
||
376 | |||
377 | =cut |
||
378 | |||
379 | sub set_bearing_symmetric |
||
380 | { |
||
381 | my( $self, $sym ) = @_; |
||
382 | # see if argument passed |
||
383 | if( $#_ > 0 ) { |
||
384 | # yes -- use value passed |
||
385 | $self->{bearing} = $sym; |
||
386 | }else{ |
||
387 | # no -- set to true |
||
388 | $self->{bearing} = 1; |
||
389 | } |
||
390 | } |
||
391 | |||
392 | =head2 set_defaults |
||
393 | |||
394 | Sets the defaults for the new method. Call with key, value pairs similar to |
||
395 | new. |
||
396 | |||
397 | $Geo::Ellipsoid->set_defaults( |
||
398 | units => 'degrees', |
||
399 | ellipsoid => 'GRS80', |
||
400 | distance_units => 'kilometer', |
||
401 | longitude => 1, |
||
402 | bearing => 0 |
||
403 | ); |
||
404 | |||
405 | Keys and string values (except for the ellipsoid identifier) may be shortened |
||
406 | to their first three letters and are case-insensitive: |
||
407 | |||
408 | $Geo::Ellipsoid->set_defaults( |
||
409 | uni => 'deg', |
||
410 | ell => 'GRS80', |
||
411 | dis => 'kil', |
||
412 | lon => 1, |
||
413 | bea => 0 |
||
414 | ); |
||
415 | |||
416 | =cut |
||
417 | |||
418 | sub set_defaults |
||
419 | { |
||
420 | my $self = shift; |
||
421 | my %args = @_; |
||
422 | foreach my $key ( keys %args ) { |
||
423 | if( $key =~ /^ell/i ) { |
||
424 | $defaults{ellipsoid} = uc $args{$key}; |
||
425 | }elsif( $key =~ /^uni/i ) { |
||
426 | $defaults{units} = $args{$key}; |
||
427 | }elsif( $key =~ /^dis/i ) { |
||
428 | $defaults{distance_units} = $args{$key}; |
||
429 | }elsif( $key =~ /^lon/i ) { |
||
430 | $defaults{longitude} = $args{$key}; |
||
431 | }elsif( $key =~ /^bea/i ) { |
||
432 | $defaults{bearing} = $args{$key}; |
||
433 | }else{ |
||
434 | croak("Geo::Ellipsoid::set_defaults called with invalid key: $key"); |
||
435 | } |
||
436 | } |
||
437 | print "Defaults set to ($defaults{ellipsoid},$defaults{units}\n" |
||
438 | if $DEBUG; |
||
439 | } |
||
440 | |||
441 | =head2 scales |
||
442 | |||
443 | Returns a list consisting of the distance unit per angle of latitude |
||
444 | and longitude (degrees or radians) at the specified latitude. |
||
445 | These values may be used for fast approximations of distance |
||
446 | calculations in the vicinity of some location. |
||
447 | |||
448 | ( $lat_scale, $lon_scale ) = $geo->scales($lat0); |
||
449 | $x = $lon_scale * ($lon - $lon0); |
||
450 | $y = $lat_scale * ($lat - $lat0); |
||
451 | |||
452 | =cut |
||
453 | |||
454 | sub scales |
||
455 | { |
||
456 | my $self = shift; |
||
457 | my $units = $self->{units}; |
||
458 | my $lat = $_[0]; |
||
459 | if( defined $lat ) { |
||
460 | $lat /= $degrees_per_radian if( $units eq 'degrees' ); |
||
461 | }else{ |
||
462 | carp("scales() method requires latitude argument; assuming 0"); |
||
463 | $lat = 0; |
||
464 | } |
||
465 | |||
466 | my $aa = $self->{equatorial}; |
||
467 | my $bb = $self->{polar}; |
||
468 | my $a2 = $aa*$aa; |
||
469 | my $b2 = $bb*$bb; |
||
470 | my $d1 = $aa * cos($lat); |
||
471 | my $d2 = $bb * sin($lat); |
||
472 | my $d3 = $d1*$d1 + $d2*$d2; |
||
473 | my $d4 = sqrt($d3); |
||
474 | my $n1 = $aa * $bb; |
||
475 | my $latscl = ( $n1 * $n1 ) / ( $d3 * $d4 * $self->{conversion} ); |
||
476 | my $lonscl = ( $aa * $d1 ) / ( $d4 * $self->{conversion} ); |
||
477 | |||
478 | if( $DEBUG ) { |
||
479 | print "lat=$lat, aa=$aa, bb=$bb\nd1=$d1, d2=$d2, d3=$d3, d4=$d4\n"; |
||
480 | print "latscl=$latscl, lonscl=$lonscl\n"; |
||
481 | } |
||
482 | |||
483 | if( $self->{units} eq 'degrees' ) { |
||
484 | $latscl /= $degrees_per_radian; |
||
485 | $lonscl /= $degrees_per_radian; |
||
486 | } |
||
487 | return ( $latscl, $lonscl ); |
||
488 | } |
||
489 | |||
490 | =head2 range |
||
491 | |||
492 | Returns the range in distance units between two specified locations given |
||
493 | as latitude, longitude pairs. |
||
494 | |||
495 | my $dist = $geo->range( $lat1, $lon1, $lat2, $lon2 ); |
||
496 | my $dist = $geo->range( @origin, @destination ); |
||
497 | |||
498 | =cut |
||
499 | |||
500 | sub range |
||
501 | { |
||
502 | my $self = shift; |
||
503 | my @args = _normalize_input($self->{units},@_); |
||
504 | my($range,$bearing) = _inverse($self,@args); |
||
505 | print "inverse(@_[1..4]) returns($range,$bearing)\n" if $DEBUG; |
||
506 | return $range; |
||
507 | } |
||
508 | |||
509 | =head2 bearing |
||
510 | |||
511 | Returns the bearing in degrees or radians from the first location to |
||
512 | the second. Zero bearing is true north. |
||
513 | |||
514 | my $bearing = $geo->bearing( $lat1, $lon1, $lat2, $lon2 ); |
||
515 | |||
516 | =cut |
||
517 | |||
518 | sub bearing |
||
519 | { |
||
520 | my $self = shift; |
||
521 | my $units = $self->{units}; |
||
522 | my @args = _normalize_input($units,@_); |
||
523 | my($range,$bearing) = _inverse($self,@args); |
||
524 | print "inverse(@args) returns($range,$bearing)\n" if $DEBUG; |
||
525 | my $t = $bearing; |
||
526 | $self->_normalize_output('bearing',$bearing); |
||
527 | print "_normalize_output($t) returns($bearing)\n" if $DEBUG; |
||
528 | return $bearing; |
||
529 | } |
||
530 | |||
531 | |||
532 | =head2 at |
||
533 | |||
534 | Returns the list (latitude,longitude) in degrees or radians that is a |
||
535 | specified range and bearing from a given location. |
||
536 | |||
537 | my( $lat2, $lon2 ) = $geo->at( $lat1, $lon1, $range, $bearing ); |
||
538 | |||
539 | =cut |
||
540 | |||
541 | sub at |
||
542 | { |
||
543 | my $self = shift; |
||
544 | my $units = $self->{units}; |
||
545 | my( $lat, $lon, $az ) = _normalize_input($units,@_[0,1,3]); |
||
546 | my $r = $_[2]; |
||
547 | print "at($lat,$lon,$r,$az)\n" if $DEBUG; |
||
548 | my( $lat2, $lon2 ) = _forward($self,$lat,$lon,$r,$az); |
||
549 | print "_forward returns ($lat2,$lon2)\n" if $DEBUG; |
||
550 | $self->_normalize_output('longitude',$lon2); |
||
551 | $self->_normalize_output('latitude',$lat2); |
||
552 | return ( $lat2, $lon2 ); |
||
553 | } |
||
554 | |||
555 | =head2 to |
||
556 | |||
557 | In list context, returns (range, bearing) between two specified locations. |
||
558 | In scalar context, returns just the range. |
||
559 | |||
560 | my( $dist, $theta ) = $geo->to( $lat1, $lon1, $lat2, $lon2 ); |
||
561 | my $dist = $geo->to( $lat1, $lon1, $lat2, $lon2 ); |
||
562 | |||
563 | =cut |
||
564 | |||
565 | sub to |
||
566 | { |
||
567 | my $self = shift; |
||
568 | my $units = $self->{units}; |
||
569 | my @args = _normalize_input($units,@_); |
||
570 | print "to($units,@args)\n" if $DEBUG; |
||
571 | my($range,$bearing) = _inverse($self,@args); |
||
572 | print "to: inverse(@args) returns($range,$bearing)\n" if $DEBUG; |
||
573 | #$bearing *= $degrees_per_radian if $units eq 'degrees'; |
||
574 | $self->_normalize_output('bearing',$bearing); |
||
575 | if( wantarray() ) { |
||
576 | return ( $range, $bearing ); |
||
577 | }else{ |
||
578 | return $range; |
||
579 | } |
||
580 | } |
||
581 | |||
582 | =head2 displacement |
||
583 | |||
584 | Returns the (x,y) displacement in distance units between the two specified |
||
585 | locations. |
||
586 | |||
587 | my( $x, $y ) = $geo->displacement( $lat1, $lon1, $lat2, $lon2 ); |
||
588 | |||
589 | NOTE: The x and y displacements are only approximations and only valid |
||
590 | between two locations that are fairly near to each other. Beyond 10 kilometers |
||
591 | or more, the concept of X and Y on a curved surface loses its meaning. |
||
592 | |||
593 | =cut |
||
594 | |||
595 | sub displacement |
||
596 | { |
||
597 | my $self = shift; |
||
598 | print "displacement(",join(',',@_),"\n" if $DEBUG; |
||
599 | my @args = _normalize_input($self->{units},@_); |
||
600 | print "call _inverse(@args)\n" if $DEBUG; |
||
601 | my( $range, $bearing ) = _inverse($self,@args); |
||
602 | print "disp: _inverse(@args) returns ($range,$bearing)\n" if $DEBUG; |
||
603 | my $x = $range * sin($bearing); |
||
604 | my $y = $range * cos($bearing); |
||
605 | return ($x,$y); |
||
606 | } |
||
607 | |||
608 | =head2 location |
||
609 | |||
610 | Returns the list (latitude,longitude) of a location at a given (x,y) |
||
611 | displacement from a given location. |
||
612 | |||
613 | my @loc = $geo->location( $lat, $lon, $x, $y ); |
||
614 | |||
615 | =cut |
||
616 | |||
617 | sub location |
||
618 | { |
||
619 | my $self = shift; |
||
620 | my $units = $self->{units}; |
||
621 | my($lat,$lon,$x,$y) = @_; |
||
622 | my $range = sqrt( $x*$x+ $y*$y ); |
||
623 | my $bearing = atan2($x,$y); |
||
624 | $bearing *= $degrees_per_radian if $units eq 'degrees'; |
||
625 | print "location($lat,$lon,$x,$y,$range,$bearing)\n" if $DEBUG; |
||
626 | return $self->at($lat,$lon,$range,$bearing); |
||
627 | } |
||
628 | |||
629 | ######################################################################## |
||
630 | # |
||
631 | # internal functions |
||
632 | # |
||
633 | # inverse |
||
634 | # |
||
635 | # Calculate the displacement from origin to destination. |
||
636 | # The input to this subroutine is |
||
637 | # ( latitude-1, longitude-1, latitude-2, longitude-2 ) in radians. |
||
638 | # |
||
639 | # Return the results as the list (range,bearing) with range in the |
||
640 | # current specified distance unit and bearing in radians. |
||
641 | |||
642 | sub _inverse() |
||
643 | { |
||
644 | my $self = shift; |
||
645 | my( $lat1, $lon1, $lat2, $lon2 ) = (@_); |
||
646 | print "_inverse($lat1,$lon1,$lat2,$lon2)\n" if $DEBUG; |
||
647 | |||
648 | my $a = $self->{equatorial}; |
||
649 | my $f = $self->{flattening}; |
||
650 | |||
651 | my $r = 1.0 - $f; |
||
652 | my $tu1 = $r * sin($lat1) / cos($lat1); |
||
653 | my $tu2 = $r * sin($lat2) / cos($lat2); |
||
654 | my $cu1 = 1.0 / ( sqrt(($tu1*$tu1) + 1.0) ); |
||
655 | my $su1 = $cu1 * $tu1; |
||
656 | my $cu2 = 1.0 / ( sqrt( ($tu2*$tu2) + 1.0 )); |
||
657 | my $s = $cu1 * $cu2; |
||
658 | my $baz = $s * $tu2; |
||
659 | my $faz = $baz * $tu1; |
||
660 | my $dlon = $lon2 - $lon1; |
||
661 | |||
662 | if( $DEBUG ) { |
||
663 | printf "lat1=%.8f, lon1=%.8f\n", $lat1, $lon1; |
||
664 | printf "lat2=%.8f, lon2=%.8f\n", $lat2, $lon2; |
||
665 | printf "r=%.8f, tu1=%.8f, tu2=%.8f\n", $r, $tu1, $tu2; |
||
666 | printf "faz=%.8f, dlon=%.8f\n", $faz, $dlon; |
||
667 | } |
||
668 | |||
669 | my $x = $dlon; |
||
670 | my $cnt = 0; |
||
671 | print "enter loop:\n" if $DEBUG; |
||
672 | my( $c2a, $c, $cx, $cy, $cz, $d, $del, $e, $sx, $sy, $y ); |
||
673 | do { |
||
674 | printf " x=%.8f\n", $x if $DEBUG; |
||
675 | $sx = sin($x); |
||
676 | $cx = cos($x); |
||
677 | $tu1 = $cu2*$sx; |
||
678 | $tu2 = $baz - ($su1*$cu2*$cx); |
||
679 | |||
680 | printf " sx=%.8f, cx=%.8f, tu1=%.8f, tu2=%.8f\n", |
||
681 | $sx, $cx, $tu1, $tu2 if $DEBUG; |
||
682 | |||
683 | $sy = sqrt( $tu1*$tu1 + $tu2*$tu2 ); |
||
684 | $cy = $s*$cx + $faz; |
||
685 | $y = atan2($sy,$cy); |
||
686 | my $sa; |
||
687 | if( $sy == 0.0 ) { |
||
688 | $sa = 1.0; |
||
689 | }else{ |
||
690 | $sa = ($s*$sx) / $sy; |
||
691 | } |
||
692 | |||
693 | printf " sy=%.8f, cy=%.8f, y=%.8f, sa=%.8f\n", $sy, $cy, $y, $sa |
||
694 | if $DEBUG; |
||
695 | |||
696 | $c2a = 1.0 - ($sa*$sa); |
||
697 | $cz = $faz + $faz; |
||
698 | if( $c2a > 0.0 ) { |
||
699 | $cz = ((-$cz)/$c2a) + $cy; |
||
700 | } |
||
701 | $e = ( 2.0 * $cz * $cz ) - 1.0; |
||
702 | $c = ( ((( (-3.0 * $c2a) + 4.0)*$f) + 4.0) * $c2a * $f )/16.0; |
||
703 | $d = $x; |
||
704 | $x = ( ($e * $cy * $c + $cz) * $sy * $c + $y) * $sa; |
||
705 | $x = ( 1.0 - $c ) * $x * $f + $dlon; |
||
706 | $del = $d - $x; |
||
707 | |||
708 | if( $DEBUG ) { |
||
709 | printf " c2a=%.8f, cz=%.8f\n", $c2a, $cz; |
||
710 | printf " e=%.8f, d=%.8f\n", $e, $d; |
||
711 | printf " (d-x)=%.8g\n", $del; |
||
712 | } |
||
713 | |||
714 | }while( (abs($del) > $eps) && ( ++$cnt <= $max_loop_count ) ); |
||
715 | |||
716 | $faz = atan2($tu1,$tu2); |
||
717 | $baz = atan2($cu1*$sx,($baz*$cx - $su1*$cu2)) + pi; |
||
718 | $x = sqrt( ((1.0/($r*$r)) -1.0 ) * $c2a+1.0 ) + 1.0; |
||
719 | $x = ($x-2.0)/$x; |
||
720 | $c = 1.0 - $x; |
||
721 | $c = (($x*$x)/4.0 + 1.0)/$c; |
||
722 | $d = ((0.375*$x*$x) - 1.0)*$x; |
||
723 | $x = $e*$cy; |
||
724 | |||
725 | if( $DEBUG ) { |
||
726 | printf "e=%.8f, cy=%.8f, x=%.8f\n", $e, $cy, $x; |
||
727 | printf "sy=%.8f, c=%.8f, d=%.8f\n", $sy, $c, $d; |
||
728 | printf "cz=%.8f, a=%.8f, r=%.8f\n", $cz, $a, $r; |
||
729 | } |
||
730 | |||
731 | $s = 1.0 - $e - $e; |
||
732 | $s = (((((((( $sy * $sy * 4.0 ) - 3.0) * $s * $cz * $d/6.0) - $x) * |
||
733 | $d /4.0) + $cz) * $sy * $d) + $y ) * $c * $a * $r; |
||
734 | |||
735 | printf "s=%.8f\n", $s if $DEBUG; |
||
736 | |||
737 | # adjust azimuth to (0,360) or (-180,180) as specified |
||
738 | if( $self->{symmetric} ) { |
||
739 | $faz += $twopi if $faz < -(pi); |
||
740 | $faz -= $twopi if $faz >= pi; |
||
741 | }else{ |
||
742 | $faz += $twopi if $faz < 0; |
||
743 | $faz -= $twopi if $faz >= $twopi; |
||
744 | } |
||
745 | |||
746 | # return result |
||
747 | my @disp = ( ($s/$self->{conversion}), $faz ); |
||
748 | print "disp = (@disp)\n" if $DEBUG; |
||
749 | return @disp; |
||
750 | } |
||
751 | |||
752 | # _forward |
||
753 | # |
||
754 | # Calculate the location (latitue,longitude) of a point |
||
755 | # given a starting point and a displacement from that |
||
756 | # point as (range,bearing) |
||
757 | # |
||
758 | sub _forward |
||
759 | { |
||
760 | my $self = shift; |
||
761 | my( $lat1, $lon1, $range, $bearing ) = @_; |
||
762 | |||
763 | if( $DEBUG ) { |
||
764 | printf "_forward(lat1=%.8f,lon1=%.8f,range=%.8f,bearing=%.8f)\n", |
||
765 | $lat1, $lon1, $range, $bearing; |
||
766 | } |
||
767 | |||
768 | my $eps = 0.5e-13; |
||
769 | |||
770 | my $a = $self->{equatorial}; |
||
771 | my $f = $self->{flattening}; |
||
772 | my $r = 1.0 - $f; |
||
773 | |||
774 | my $tu = $r * sin($lat1) / cos($lat1); |
||
775 | my $faz = $bearing; |
||
776 | my $s = $self->{conversion} * $range; |
||
777 | my $sf = sin($faz); |
||
778 | my $cf = cos($faz); |
||
779 | |||
780 | my $baz = 0.0; |
||
781 | $baz = 2.0 * atan2($tu,$cf) if( $cf != 0.0 ); |
||
782 | |||
783 | my $cu = 1.0 / sqrt(1.0 + $tu*$tu); |
||
784 | my $su = $tu * $cu; |
||
785 | my $sa = $cu * $sf; |
||
786 | my $c2a = 1.0 - ($sa*$sa); |
||
787 | my $x = 1.0 + sqrt( (((1.0/($r*$r)) - 1.0 )*$c2a) +1.0); |
||
788 | $x = ($x-2.0)/$x; |
||
789 | my $c = 1.0 - $x; |
||
790 | $c = ((($x*$x)/4.0) + 1.0)/$c; |
||
791 | my $d = $x * ((0.375*$x*$x)-1.0); |
||
792 | $tu = (($s/$r)/$a)/$c; |
||
793 | my $y = $tu; |
||
794 | |||
795 | if( $DEBUG ) { |
||
796 | printf "r=%.8f, tu=%.8f, faz=%.8f\n", $r, $tu, $faz; |
||
797 | printf "baz=%.8f, sf=%.8f, cf=%.8f\n", $baz, $sf, $cf; |
||
798 | printf "cu=%.8f, su=%.8f, sa=%.8f\n", $cu, $su, $sa; |
||
799 | printf "x=%.8f, c=%.8f, y=%.8f\n", $x, $c, $y; |
||
800 | } |
||
801 | |||
802 | my( $cy, $cz, $e, $sy ); |
||
803 | do { |
||
804 | $sy = sin($y); |
||
805 | $cy = cos($y); |
||
806 | $cz = cos($baz+$y); |
||
807 | $e = (2.0*$cz*$cz)-1.0; |
||
808 | $c = $y; |
||
809 | $x = $e * $cy; |
||
810 | $y = (2.0 * $e) - 1.0; |
||
811 | $y = ((((((((($sy*$sy*4.0)-3.0)*$y*$cz*$d)/6.0)+$x)*$d)/4.0)-$cz)*$sy*$d) + |
||
812 | $tu; |
||
813 | } while( abs($y-$c) > $eps ); |
||
814 | |||
815 | $baz = ($cu*$cy*$cf) - ($su*$sy); |
||
816 | $c = $r*sqrt(($sa*$sa) + ($baz*$baz)); |
||
817 | $d = $su*$cy + $cu*$sy*$cf; |
||
818 | my $lat2 = atan2($d,$c); |
||
819 | $c = $cu*$cy - $su*$sy*$cf; |
||
820 | $x = atan2($sy*$sf,$c); |
||
821 | $c = (((((-3.0*$c2a)+4.0)*$f)+4.0)*$c2a*$f)/16.0; |
||
822 | $d = (((($e*$cy*$c) + $cz)*$sy*$c)+$y)*$sa; |
||
823 | my $lon2 = $lon1 + $x - (1.0-$c)*$d*$f; |
||
824 | #$baz = atan2($sa,$baz) + pi; |
||
825 | |||
826 | # return result |
||
827 | return ($lat2,$lon2); |
||
828 | |||
829 | } |
||
830 | |||
831 | # _normalize_input |
||
832 | # |
||
833 | # Normalize a set of input angle values by converting to |
||
834 | # radians if given in degrees and by converting to the |
||
835 | # range [0,2pi), i.e. greater than or equal to zero and |
||
836 | # less than two pi. |
||
837 | # |
||
838 | sub _normalize_input |
||
839 | { |
||
840 | my $units = shift; |
||
841 | my @args = @_; |
||
842 | return map { |
||
843 | $_ = deg2rad($_) if $units eq 'degrees'; |
||
844 | while( $_ < 0 ) { $_ += $twopi } |
||
845 | while( $_ >= $twopi ) { $_ -= $twopi } |
||
846 | $_ |
||
847 | } @args; |
||
848 | } |
||
849 | |||
850 | # _normalize_output |
||
851 | # |
||
852 | # Normalize a set of output angle values by converting to |
||
853 | # degrees if needed and by converting to the range [-pi,+pi) or |
||
854 | # [0,2pi) as needed. |
||
855 | # |
||
856 | sub _normalize_output |
||
857 | { |
||
858 | my $self = shift; |
||
859 | my $elem = shift; # 'bearing' or 'longitude' |
||
860 | # adjust remaining input values by reference |
||
861 | for ( @_ ) { |
||
862 | if( $self->{$elem} ) { |
||
863 | # normalize to range [-pi,pi) |
||
864 | while( $_ < -(pi) ) { $_ += $twopi } |
||
865 | while( $_ >= pi ) { $_ -= $twopi } |
||
866 | }else{ |
||
867 | # normalize to range [0,2*pi) |
||
868 | while( $_ < 0 ) { $_ += $twopi } |
||
869 | while( $_ >= $twopi ) { $_ -= $twopi } |
||
870 | } |
||
871 | $_ = rad2deg($_) if $self->{units} eq 'degrees'; |
||
872 | } |
||
873 | } |
||
874 | |||
875 | =head1 DEFINED ELLIPSOIDS |
||
876 | |||
877 | The following ellipsoids are defined in Geo::Ellipsoid, with the |
||
878 | semi-major axis in meters and the reciprocal flattening as shown. |
||
879 | The default ellipsoid is WGS84. |
||
880 | |||
881 | Ellipsoid Semi-Major Axis (m.) 1/Flattening |
||
882 | --------- ------------------- --------------- |
||
883 | AIRY 6377563.396 299.3249646 |
||
884 | AIRY-MODIFIED 6377340.189 299.3249646 |
||
885 | AUSTRALIAN 6378160.0 298.25 |
||
886 | BESSEL-1841 6377397.155 299.1528128 |
||
887 | CLARKE-1880 6378249.145 293.465 |
||
888 | EVEREST-1830 6377276.345 290.8017 |
||
889 | EVEREST-MODIFIED 6377304.063 290.8017 |
||
890 | FISHER-1960 6378166.0 298.3 |
||
891 | FISHER-1968 6378150.0 298.3 |
||
892 | GRS80 6378137.0 298.25722210088 |
||
893 | HOUGH-1956 6378270.0 297.0 |
||
894 | HAYFORD 6378388.0 297.0 |
||
895 | IAU76 6378140.0 298.257 |
||
896 | KRASSOVSKY-1938 6378245.0 298.3 |
||
897 | NAD27 6378206.4 294.9786982138 |
||
898 | NWL-9D 6378145.0 298.25 |
||
899 | SOUTHAMERICAN-1969 6378160.0 298.25 |
||
900 | SOVIET-1985 6378136.0 298.257 |
||
901 | WGS72 6378135.0 298.26 |
||
902 | WGS84 6378137.0 298.257223563 |
||
903 | |||
904 | =head1 LIMITATIONS |
||
905 | |||
906 | The methods should not be used on points which are too near the poles |
||
907 | (above or below 89 degrees), and should not be used on points which |
||
908 | are antipodal, i.e., exactly on opposite sides of the ellipsoid. The |
||
909 | methods will not return valid results in these cases. |
||
910 | |||
911 | =head1 ACKNOWLEDGEMENTS |
||
912 | |||
913 | The conversion algorithms used here are Perl translations of Fortran |
||
914 | routines written by LCDR S<L. Pfeifer> NGS Rockville MD that implement |
||
915 | S<T. Vincenty's> Modified Rainsford's method with Helmert's elliptical |
||
916 | terms as published in "Direct and Inverse Solutions of Ellipsoid on |
||
917 | the Ellipsoid with Application of Nested Equations", S<T. Vincenty,> |
||
918 | Survey Review, April 1975. |
||
919 | |||
920 | The Fortran source code files inverse.for and forward.for |
||
921 | may be obtained from |
||
922 | |||
923 | ftp://ftp.ngs.noaa.gov/pub/pcsoft/for_inv.3d/source/ |
||
924 | |||
925 | =head1 AUTHOR |
||
926 | |||
927 | Jim Gibson, C<< <Jim@Gibson.org> >> |
||
928 | |||
929 | =head1 BUGS |
||
930 | |||
931 | See LIMITATIONS, above. |
||
932 | |||
933 | Please report any bugs or feature requests to |
||
934 | C<bug-geo-ellipsoid@rt.cpan.org>, or through the web interface at |
||
935 | L<http://rt.cpan.org/NoAuth/ReportBug.html?Queue=Geo-Ellipsoid>. |
||
936 | |||
937 | =head1 COPYRIGHT & LICENSE |
||
938 | |||
939 | Copyright 2005-2008 Jim Gibson, all rights reserved. |
||
940 | |||
941 | This program is free software; you can redistribute it and/or modify it |
||
942 | under the same terms as Perl itself. |
||
943 | |||
944 | =head1 SEE ALSO |
||
945 | |||
946 | Geo::Distance, Geo::Ellipsoids |
||
947 | |||
948 | =cut |
||
949 | |||
950 | 1; # End of Geo::Ellipsoid |