Subversion Repositories Projects

#       Geo::Ellipsoid

#

#       This package implements an Ellipsoid class to perform latitude 

#       and longitude calculations on the surface of an ellipsoid.

#

#       This is a Perl conversion of existing Fortran code (see

#       ACKNOWLEDGEMENTS) and the author of this class makes no 

#       claims of originality. Nor can he even vouch for the 

#       results of the calculations, although they do seem to 

#       work for him and have been tested against other methods.

package Geo::Ellipsoid;

use warnings;

use strict;

use 5.006_00;

use Scalar::Util 'looks_like_number';

use Math::Trig;

use Carp;

=head1 NAME

Geo::Ellipsoid - Longitude and latitude calculations using an ellipsoid model.

=head1 VERSION

Version 1.12, released July 4, 2008.

=cut

our $VERSION = '1.12';

our $DEBUG = 0;

=head1 SYNOPSIS

  use Geo::Ellipsoid;

  $geo = Geo::Ellipsoid->new(ellipsoid=>'NAD27', units=>'degrees');

  @origin = ( 37.619002, -122.374843 );    # SFO

  @dest = ( 33.942536, -118.408074 );      # LAX

  ( $range, $bearing ) = $geo->to( @origin, @dest );

  ($lat,$lon) = $geo->at( @origin, 2000, 45.0 );

  ( $x, $y ) = $geo->displacement( @origin, $lat, $lon );

  @pos = $geo->location( $lat, $lon, $x, $y );

=head1 DESCRIPTION

Geo::Ellipsoid performs geometrical calculations on the surface of

an ellipsoid. An ellipsoid is a three-dimension object formed from

the rotation of an ellipse about one of its axes. The approximate

shape of the earth is an ellipsoid, so Geo::Ellipsoid can accurately

calculate distance and bearing between two widely-separated locations

on the earth's surface.

The shape of an ellipsoid is defined by the lengths of its

semi-major and semi-minor axes. The shape may also be specifed by

the flattening ratio C<f> as:

    f = ( semi-major - semi-minor ) / semi-major

which, since f is a small number, is normally given as the reciprocal

of the flattening C<1/f>.

The shape of the earth has been surveyed and estimated differently

at different times over the years. The two most common sets of values

used to describe the size and shape of the earth in the United States

are 'NAD27', dating from 1927, and 'WGS84', from 1984. United States

Geological Survey topographical maps, for example, use one or the

other of these values, and commonly-available Global Positioning

System (GPS) units can be set to use one or the other.

See L<"DEFINED ELLIPSOIDS"> below for the ellipsoid survey values

that may be selected for use by Geo::Ellipsoid.

=cut

# class data and constants

our $degrees_per_radian = 180/pi;

our $eps = 1.0e-23;

our $max_loop_count = 20;

our $twopi = 2 * pi;

our $halfpi = pi/2;

our %defaults = (

  ellipsoid => 'WGS84',

  units => 'radians',

  distance_units => 'meter',

  longitude => 0,

  latitude => 1,        # allows use of _normalize_output

  bearing => 0,

);

our %distance = (

  'foot'      => 0.3048,

  'kilometer' => 1_000,

  'meter'     => 1.0,

  'mile'      => 1_609.344,

  'nm'        => 1_852,

);

# set of ellipsoids that can be used.

# values are 

#  1) a = semi-major (equatorial) radius of Ellipsoid 

#  2) 1/f = reciprocal of flattening (f), the ratio of the semi-minor 

#     (polar) radius to the semi-major (equatorial) axis, or

#     polar radius = equatorial radius * ( 1 - f )

our %ellipsoids = (

    'AIRY'               => [ 6377563.396, 299.3249646     ],

    'AIRY-MODIFIED'      => [ 6377340.189, 299.3249646     ],

    'AUSTRALIAN'         => [ 6378160.0,   298.25          ],

    'BESSEL-1841'        => [ 6377397.155, 299.1528128     ],

    'CLARKE-1880'        => [ 6378249.145, 293.465         ],

    'EVEREST-1830'       => [ 6377276.345, 300.8017        ],

    'EVEREST-MODIFIED'   => [ 6377304.063, 300.8017        ],

    'FISHER-1960'        => [ 6378166.0,   298.3           ],

    'FISHER-1968'        => [ 6378150.0,   298.3           ],

    'GRS80'              => [ 6378137.0,   298.25722210088 ],

    'HOUGH-1956'         => [ 6378270.0,   297.0           ],

    'HAYFORD'            => [ 6378388.0,   297.0           ],

    'IAU76'              => [ 6378140.0,   298.257         ],

    'KRASSOVSKY-1938'    => [ 6378245.0,   298.3           ],

    'NAD27'              => [ 6378206.4,   294.9786982138  ],

    'NWL-9D'             => [ 6378145.0,   298.25          ],

    'SOUTHAMERICAN-1969' => [ 6378160.0,   298.25          ],

    'SOVIET-1985'        => [ 6378136.0,   298.257         ],

    'WGS72'              => [ 6378135.0,   298.26          ],

    'WGS84'              => [ 6378137.0,   298.257223563   ],

);

=head1 CONSTRUCTOR

=head2 new

The new() constructor may be called with a hash list to set the value of the

ellipsoid to be used, the value of the units to be used for angles and

distances, and whether or not the output range of longitudes and bearing

angles should be symmetric around zero or always greater than zero.  

The initial default constructor is equivalent to the following:

    my $geo = Geo::Ellipsoid->new(

      ellipsoid => 'WGS84',

      units => 'radians' ,

      distance_units => 'meter',

      longitude => 0,

      bearing => 0,

    );

The constructor arguments may be of any case and, with the exception of

the ellipsoid value, abbreviated to their first three characters.

Thus, ( UNI => 'DEG', DIS => 'FEE', Lon => 1, ell => 'NAD27', bEA => 0 )

is valid.

=cut

sub new

{

  my( $class, %args ) = @_;

  my $self = {%defaults};

  print "new: @_\n" if $DEBUG;

  foreach my $key ( keys %args ) {

    my $val = $args{$key};

    if( $key =~ /^ell/i ) {

      $self->{ellipsoid} = uc $args{$key};

    }elsif( $key =~ /^uni/i ) {

      $self->{units} = $args{$key};

    }elsif( $key =~ /^dis/i ) {

      $self->{distance_units} = $args{$key};

    }elsif( $key =~ /^lon/i ) {

      $self->{longitude} = $args{$key};

    }elsif( $key =~ /^bea/i ) {

      $self->{bearing} = $args{$key};

    }else{

      carp("Unknown argument: $key => $args{$key}");

    }

  }

  set_units($self,$self->{units});

  set_ellipsoid($self,$self->{ellipsoid});

  set_distance_unit($self,$self->{distance_units});

  set_longitude_symmetric($self,$self->{longitude});

  set_bearing_symmetric($self,$self->{bearing});

  print

    "Ellipsoid(units=>$self->{units},distance_units=>" .

    "$self->{distance_units},ellipsoid=>$self->{ellipsoid}," .

    "longitude=>$self->{longitude},bearing=>$self->{bearing})\n" if $DEBUG;

  bless $self,$class;

  return $self;

}

=head1 METHODS

=head2 set_units

Set the angle units used by the Geo::Ellipsoid object. The units may

also be set in the constructor of the object. The allowable values are

'degrees' or 'radians'. The default is 'radians'. The units value is

not case sensitive and may be abbreviated to 3 letters. The units of

angle apply to both input and output latitude, longitude, and bearing

values.

    $geo->set_units('degrees');

=cut

sub set_units

{

  my $self = shift;

  my $units = shift;

  if( $units =~ /deg/i ) {

    $units = 'degrees';

  }elsif( $units =~ /rad/i ) {

    $units = 'radians';

  }else{

    croak("Invalid units specifier '$units' - please use either " .

      "degrees or radians (the default)") unless $units =~ /rad/i;

  }

  $self->{units} = $units;

}

=head2 set_distance_unit

Set the distance unit used by the Geo::Ellipsoid object. The unit of

distance may also be set in the constructor of the object. The recognized

values are 'meter', 'kilometer', 'mile', 'nm' (nautical mile), or 'foot'.

The default is 'meter'. The value is not case sensitive and may be

abbreviated to 3 letters.

    $geo->set_distance_unit('kilometer');

For any other unit of distance not recogized by this method, pass a

numerical argument representing the length of the distance unit in

meters. For example, to use units of furlongs, call

    $geo->set_distance_unit(201.168);

The distance conversion factors used by this module are as follows:

  Unit          Units per meter

  --------      ---------------

  foot             0.3048

  kilometer     1000.0

  mile          1609.344

  nm            1852.0

=cut

sub set_distance_unit

{

  my $self = shift;

  my $unit = shift;

  print "distance unit = $unit\n" if $DEBUG;

  my $conversion = 0;

  if( defined $unit ) {

    my( $key, $val );

    while( ($key,$val) = each %distance ) {

      my $re = substr($key,0,3);

      print "trying ($key,$re,$val)\n" if $DEBUG;

      if( $unit =~ /^$re/i ) {

        $self->{distance_units} = $unit;

        $conversion = $val;    

        # finish iterating to reset 'each' function call

        while( each %distance ) {}

        last;

      }

    }

    if( $conversion == 0 ) {

      if( looks_like_number($unit) ) {

        $conversion = $unit;

      }else{

        carp("Unknown argument to set_distance_unit: $unit\nAssuming meters");

        $conversion = 1.0;

      }

    }

  }else{

    carp("Missing or undefined argument to set_distance_unit: ".

      "$unit\nAssuming meters");

    $conversion = 1.0;

  }

  $self->{conversion} = $conversion;

}

=head2 set_ellipsoid

Set the ellipsoid to be used by the Geo::Ellipsoid object. See

L<"DEFINED ELLIPSOIDS"> below for the allowable values. The value

may also be set by the constructor. The default value is 'WGS84'.

    $geo->set_ellipsoid('NAD27');

=cut

sub set_ellipsoid

{

  my $self = shift;

  my $ellipsoid = uc shift || $defaults{ellipsoid};

  print "  set ellipsoid to $ellipsoid\n" if $DEBUG;

  unless( exists $ellipsoids{$ellipsoid} ) {

    croak("Ellipsoid $ellipsoid does not exist - please use " .

      "set_custom_ellipsoid to use an ellipsoid not in valid set");

  }

  $self->{ellipsoid} = $ellipsoid;

  my( $major, $recip ) = @{$ellipsoids{$ellipsoid}};

  $self->{equatorial} = $major;

  if( $recip == 0 ) {

    carp("Infinite flattening specified by ellipsoid -- assuming a sphere");

    $self->{polar} = $self->{equatorial};

    $self->{flattening} = 0;

    $self->{eccentricity} = 0;

  }else{

    $self->{flattening} = ( 1.0 / $ellipsoids{$ellipsoid}[1]);

    $self->{polar} = $self->{equatorial} * ( 1.0  - $self->{flattening} );

    $self->{eccentricity} = sqrt( 2.0 * $self->{flattening} -

      ( $self->{flattening} * $self->{flattening} ) );

  }

}

=head2 set_custom_ellipsoid

Sets the ellipsoid parameters to the specified ( major semiaxis and

reciprocal flattening. A zero value for the reciprocal flattening

will result in a sphere for the ellipsoid, and a warning message

will be issued.

    $geo->set_custom_ellipsoid( 'sphere', 6378137, 0 );

=cut

sub set_custom_ellipsoid

{

  my $self = shift;

  my( $name, $major, $recip ) = @_;

  $name = uc $name;

  $recip = 0 unless defined $recip;

  if( $major ) {

    $ellipsoids{$name} = [ $major, $recip ];

  }else{  

    croak("set_custom_ellipsoid called without semi-major radius parameter");

  }

  set_ellipsoid($self,$name);

}

=head2 set_longitude_symmetric

If called with no argument or a true argument, sets the range of output

values for longitude to be in the range [-pi,+pi) radians.  If called with

a false or undefined argument, sets the output angle range to be

[0,2*pi) radians.

    $geo->set_longitude_symmetric(1);

=cut

sub set_longitude_symmetric

{

  my( $self, $sym ) = @_;

  # see if argument passed

  if( $#_ > 0 ) {

    # yes -- use value passed

    $self->{longitude} = $sym;

  }else{

    # no -- set to true

    $self->{longitude} = 1;

  }

}

=head2 set_bearing_symmetric

If called with no argument or a true argument, sets the range of output

values for bearing to be in the range [-pi,+pi) radians.  If called with

a false or undefined argument, sets the output angle range to be

[0,2*pi) radians.

    $geo->set_bearing_symmetric(1);

=cut

sub set_bearing_symmetric

{

  my( $self, $sym ) = @_;

  # see if argument passed

  if( $#_ > 0 ) {

    # yes -- use value passed

    $self->{bearing} = $sym;

  }else{

    # no -- set to true

    $self->{bearing} = 1;

  }

}

=head2 set_defaults

Sets the defaults for the new method. Call with key, value pairs similar to

new.

    $Geo::Ellipsoid->set_defaults(

      units => 'degrees',

      ellipsoid => 'GRS80',

      distance_units => 'kilometer',

      longitude => 1,

      bearing => 0

    );

Keys and string values (except for the ellipsoid identifier) may be shortened

to their first three letters and are case-insensitive:

    $Geo::Ellipsoid->set_defaults(

      uni => 'deg',

      ell => 'GRS80',

      dis => 'kil',

      lon => 1,

      bea => 0

    );

=cut

sub set_defaults

{

  my $self = shift;

  my %args = @_;

  foreach my $key ( keys %args ) {

    if( $key =~ /^ell/i ) {

      $defaults{ellipsoid} = uc $args{$key};

    }elsif( $key =~ /^uni/i ) {

      $defaults{units} = $args{$key};

    }elsif( $key =~ /^dis/i ) {

      $defaults{distance_units} = $args{$key};

    }elsif( $key =~ /^lon/i ) {

      $defaults{longitude} = $args{$key};

    }elsif( $key =~ /^bea/i ) {

      $defaults{bearing} = $args{$key};

    }else{  

      croak("Geo::Ellipsoid::set_defaults called with invalid key: $key");

    }

  }

  print "Defaults set to ($defaults{ellipsoid},$defaults{units}\n"

    if $DEBUG;

}

=head2 scales

Returns a list consisting of the distance unit per angle of latitude

and longitude (degrees or radians) at the specified latitude.

These values may be used for fast approximations of distance

calculations in the vicinity of some location.

    ( $lat_scale, $lon_scale ) = $geo->scales($lat0);

    $x = $lon_scale * ($lon - $lon0);

    $y = $lat_scale * ($lat - $lat0);

=cut

sub scales

{

  my $self = shift;

  my $units = $self->{units};

  my $lat = $_[0];

  if( defined $lat ) {

    $lat /= $degrees_per_radian if( $units eq 'degrees' );

  }else{

    carp("scales() method requires latitude argument; assuming 0");

    $lat = 0;

  }

  my $aa = $self->{equatorial};

  my $bb = $self->{polar};

  my $a2 = $aa*$aa;

  my $b2 = $bb*$bb;

  my $d1 = $aa * cos($lat);

  my $d2 = $bb * sin($lat);

  my $d3 = $d1*$d1 + $d2*$d2;

  my $d4 = sqrt($d3);

  my $n1 = $aa * $bb;

  my $latscl = ( $n1 * $n1 ) / ( $d3 * $d4 * $self->{conversion} );

  my $lonscl = ( $aa * $d1 ) / ( $d4 * $self->{conversion} );

  if( $DEBUG ) {

    print "lat=$lat, aa=$aa, bb=$bb\nd1=$d1, d2=$d2, d3=$d3, d4=$d4\n";

    print "latscl=$latscl, lonscl=$lonscl\n";

  }

  if( $self->{units} eq 'degrees' ) {

    $latscl /= $degrees_per_radian;

    $lonscl /= $degrees_per_radian;

  }

  return ( $latscl, $lonscl );

}

=head2 range

Returns the range in distance units between two specified locations given

as latitude, longitude pairs.

    my $dist = $geo->range( $lat1, $lon1, $lat2, $lon2 );

    my $dist = $geo->range( @origin, @destination );

=cut

sub range

{

  my $self = shift;

  my @args = _normalize_input($self->{units},@_);

  my($range,$bearing) = _inverse($self,@args);

  print "inverse(@_[1..4]) returns($range,$bearing)\n" if $DEBUG;

  return $range;

}

=head2 bearing

Returns the bearing in degrees or radians from the first location to

the second. Zero bearing is true north.

    my $bearing = $geo->bearing( $lat1, $lon1, $lat2, $lon2 );

=cut

sub bearing

{

  my $self = shift;

  my $units = $self->{units};

  my @args = _normalize_input($units,@_);

  my($range,$bearing) = _inverse($self,@args);

  print "inverse(@args) returns($range,$bearing)\n" if $DEBUG;

  my $t = $bearing;

  $self->_normalize_output('bearing',$bearing);

  print "_normalize_output($t) returns($bearing)\n" if $DEBUG;

  return $bearing;

}

=head2 at

Returns the list (latitude,longitude) in degrees or radians that is a

specified range and bearing from a given location.

    my( $lat2, $lon2 ) = $geo->at( $lat1, $lon1, $range, $bearing );

=cut

sub at

{

  my $self = shift;

  my $units = $self->{units};

  my( $lat, $lon, $az ) = _normalize_input($units,@_[0,1,3]);

  my $r = $_[2];

  print "at($lat,$lon,$r,$az)\n" if $DEBUG;

  my( $lat2, $lon2 ) = _forward($self,$lat,$lon,$r,$az);

  print "_forward returns ($lat2,$lon2)\n" if $DEBUG;

  $self->_normalize_output('longitude',$lon2);

  $self->_normalize_output('latitude',$lat2);

  return ( $lat2, $lon2 );  

}

=head2 to

In list context, returns (range, bearing) between two specified locations.

In scalar context, returns just the range.

    my( $dist, $theta ) = $geo->to( $lat1, $lon1, $lat2, $lon2 );

    my $dist = $geo->to( $lat1, $lon1, $lat2, $lon2 );

=cut

sub to

{

  my $self = shift;

  my $units = $self->{units};

  my @args = _normalize_input($units,@_);

  print "to($units,@args)\n" if $DEBUG;

  my($range,$bearing) = _inverse($self,@args);

  print "to: inverse(@args) returns($range,$bearing)\n" if $DEBUG;

  #$bearing *= $degrees_per_radian if $units eq 'degrees';

  $self->_normalize_output('bearing',$bearing);

  if( wantarray() ) {

    return ( $range, $bearing );

  }else{

    return $range;

  }

}

=head2 displacement

Returns the (x,y) displacement in distance units between the two specified

locations.

    my( $x, $y ) = $geo->displacement( $lat1, $lon1, $lat2, $lon2 );

NOTE: The x and y displacements are only approximations and only valid

between two locations that are fairly near to each other. Beyond 10 kilometers

or more, the concept of X and Y on a curved surface loses its meaning.

=cut

sub displacement

{

  my $self = shift;

  print "displacement(",join(',',@_),"\n" if $DEBUG;

  my @args = _normalize_input($self->{units},@_);

  print "call _inverse(@args)\n" if $DEBUG;

  my( $range, $bearing ) = _inverse($self,@args);

  print "disp: _inverse(@args) returns ($range,$bearing)\n" if $DEBUG;

  my $x = $range * sin($bearing);

  my $y = $range * cos($bearing);

  return ($x,$y);

}

=head2 location

Returns the list (latitude,longitude) of a location at a given (x,y)

displacement from a given location.

        my @loc = $geo->location( $lat, $lon, $x, $y );

=cut

sub location

{

  my $self = shift;

  my $units = $self->{units};

  my($lat,$lon,$x,$y) = @_;

  my $range = sqrt( $x*$x+ $y*$y );

  my $bearing = atan2($x,$y);

  $bearing *= $degrees_per_radian if $units eq 'degrees';

  print "location($lat,$lon,$x,$y,$range,$bearing)\n" if $DEBUG;

  return $self->at($lat,$lon,$range,$bearing);

}

########################################################################

#

#      internal functions

#

#       inverse

#

#       Calculate the displacement from origin to destination.

#       The input to this subroutine is 

#         ( latitude-1, longitude-1, latitude-2, longitude-2 ) in radians.

#

#       Return the results as the list (range,bearing) with range in the

#       current specified distance unit and bearing in radians.

sub _inverse()

{

  my $self = shift;

  my( $lat1, $lon1, $lat2, $lon2 ) = (@_);

  print "_inverse($lat1,$lon1,$lat2,$lon2)\n" if $DEBUG;

  my $a = $self->{equatorial};

  my $f = $self->{flattening};

  my $r = 1.0 - $f;

  my $tu1 = $r * sin($lat1) / cos($lat1);

  my $tu2 = $r * sin($lat2) / cos($lat2);

  my $cu1 = 1.0 / ( sqrt(($tu1*$tu1) + 1.0) );

  my $su1 = $cu1 * $tu1;

  my $cu2 = 1.0 / ( sqrt( ($tu2*$tu2) + 1.0 ));

  my $s = $cu1 * $cu2;

  my $baz = $s * $tu2;

  my $faz = $baz * $tu1;

  my $dlon = $lon2 - $lon1;

  if( $DEBUG ) {

    printf "lat1=%.8f, lon1=%.8f\n", $lat1, $lon1;

    printf "lat2=%.8f, lon2=%.8f\n", $lat2, $lon2;

    printf "r=%.8f, tu1=%.8f, tu2=%.8f\n", $r, $tu1, $tu2;

    printf "faz=%.8f, dlon=%.8f\n", $faz, $dlon;

  }

  my $x = $dlon;

  my $cnt = 0;

  print "enter loop:\n" if $DEBUG;

  my( $c2a, $c, $cx, $cy, $cz, $d, $del, $e, $sx, $sy, $y );

  do {

    printf "  x=%.8f\n", $x if $DEBUG;

    $sx = sin($x);

    $cx = cos($x);

    $tu1 = $cu2*$sx;

    $tu2 = $baz - ($su1*$cu2*$cx);

    printf "    sx=%.8f, cx=%.8f, tu1=%.8f, tu2=%.8f\n",

      $sx, $cx, $tu1, $tu2 if $DEBUG;

    $sy = sqrt( $tu1*$tu1 + $tu2*$tu2 );

    $cy = $s*$cx + $faz;

    $y = atan2($sy,$cy);

    my $sa;

    if( $sy == 0.0 ) {

      $sa = 1.0;

    }else{

      $sa = ($s*$sx) / $sy;

    }

    printf "    sy=%.8f, cy=%.8f, y=%.8f, sa=%.8f\n", $sy, $cy, $y, $sa

      if $DEBUG;

    $c2a = 1.0 - ($sa*$sa);

    $cz = $faz + $faz;

    if( $c2a > 0.0 ) {

      $cz = ((-$cz)/$c2a) + $cy;

    }

    $e = ( 2.0 * $cz * $cz ) - 1.0;

    $c = ( ((( (-3.0 * $c2a) + 4.0)*$f) + 4.0) * $c2a * $f )/16.0;

    $d = $x;

    $x = ( ($e * $cy * $c + $cz) * $sy * $c + $y) * $sa;

    $x = ( 1.0 - $c ) * $x * $f + $dlon;

    $del = $d - $x;

    if( $DEBUG ) {

      printf "    c2a=%.8f, cz=%.8f\n", $c2a, $cz;

      printf "    e=%.8f, d=%.8f\n", $e, $d;

      printf "    (d-x)=%.8g\n", $del;

    }

  }while( (abs($del) > $eps) && ( ++$cnt <= $max_loop_count ) );

  $faz = atan2($tu1,$tu2);

  $baz = atan2($cu1*$sx,($baz*$cx - $su1*$cu2)) + pi;

  $x = sqrt( ((1.0/($r*$r)) -1.0 ) * $c2a+1.0 ) + 1.0;

  $x = ($x-2.0)/$x;

  $c = 1.0 - $x;

  $c = (($x*$x)/4.0 + 1.0)/$c;

  $d = ((0.375*$x*$x) - 1.0)*$x;

  $x = $e*$cy;

  if( $DEBUG ) {

    printf "e=%.8f, cy=%.8f, x=%.8f\n", $e, $cy, $x;

    printf "sy=%.8f, c=%.8f, d=%.8f\n", $sy, $c, $d;

    printf "cz=%.8f, a=%.8f, r=%.8f\n", $cz, $a, $r;

  }

  $s = 1.0 - $e - $e;

  $s = (((((((( $sy * $sy * 4.0 ) - 3.0) * $s * $cz * $d/6.0) - $x) *