File Coverage

blib/lib/Geo/Ellipsoids.pm

Criterion	Covered	Total	%
statement	113	138	81.8
branch	40	58	68.9
condition	3	5	60.0
subroutine	24	29	82.7
pod	23	23	100.0
total	203	253	80.2

package Geo::Ellipsoids;
use strict;
use warnings;
use Geo::Constants qw{PI};
use Geo::Functions qw{rad_deg};

our $VERSION       = '0.17';
our $DEFAULT_ELIPS = 'WGS84';

=head1 NAME

Geo::Ellipsoids - Package for standard Geo:: ellipsoid a, b, f and 1/f values.

=head1 SYNOPSIS

  use Geo::Ellipsoids;
  my $obj = Geo::Ellipsoids->new();
  $obj->set('WGS84'); #default
  print "a=", $obj->a, "\n";
  print "b=", $obj->b, "\n";
  print "f=", $obj->f, "\n";
  print "i=", $obj->i, "\n";
  print "e=", $obj->e, "\n";
  print "n=", $obj->n(45), "\n";

=head1 DESCRIPTION

Package for standard Geo:: ellipsoid a, b, f and 1/f values.

=head1 CONSTRUCTOR

=head2 new

The new() constructor may be called with any parameter that is appropriate to the set method.

  my $obj = Geo::Ellipsoid->new();

=cut

sub new {
  my $this  = shift;
  my $class = ref($this) ? ref($this) : $this;
  my $self  = {};
  bless $self, $class;
  $self->initialize(@_);
  return $self;
}

=head1 METHODS

=head2 initialize

=cut

sub initialize {
  my $self  = shift;
  my $param = shift;
  $self->set($param);
}

=head2 set

Method sets the current ellipsoid.  This method is called when the object is constructed (default is WGS84).

  $obj->set(); #default WGS84
  $obj->set('Clarke 1866'); #All built in ellipsoids are stored in meters
  $obj->set({a=>1, b=>1});  #Custom Sphere 1 unit radius

=cut

sub set {
  my $self  = shift;
  my $param = shift || $DEFAULT_ELIPS;
  undef($self->{'shortname'});
  undef($self->{'longname'});
  if ('HASH' eq ref($param)) {
    return $self->_setref($param);
  } elsif ('' eq ref($param)) {
    return $self->_setname($param);
  } else {
    die('Error: Parameter must be the name of an ellipsoid or a hash reference');
  }
}

=head2 list

Method returns a list of known ellipsoid names.

  my @list=$obj->list;

  my $list=$obj->list;
  while (@$list) {
    print "$_\n";
  }

=cut

sub list {
  my $self = shift;
  my $data = $self->data;
  my @keys = keys %$data;
  return wantarray ? @keys : \@keys;
}

=head2 a

Method returns the value of the semi-major axis.

  my $a=$obj->a;

=cut

sub a {
  my $self = shift;
  return $self->{'a'} || die(q{Error: $self->{'a'} must be defined here});
}

=head2 b

Method returns the value of the semi-minor axis.

  my $b=$obj->b;  #b=a(1-f)

=cut

sub b {
  my $self = shift;
  if (defined $self->{'b'}) {
    return $self->{'b'};
  } elsif (defined $self->{'f'}) {
    return $self->{'a'}*(1-$self->{'f'});
  } elsif (defined $self->{'i'}) {
    return $self->{'a'}*(1-1/$self->{'i'});
  } else {
    return undef();
  }
}

=head2 f

Method returns the value of flatting

  my $f=$obj->f;  #f=(a-b)/a

=cut

sub f {
  my $self = shift;
  if (defined $self->{'f'}) {
    return $self->{'f'};
  } elsif (defined $self->{'b'}) {
    return ($self->{'a'}-$self->{'b'})/$self->{'a'};
  } elsif (defined $self->{'i'}) {
    return 1/$self->{'i'};
  } else {
    return undef();
  }
}

=head2 i

Method returns the value of the inverse flatting

  my $i = $obj->i; #i=1/f=a/(a-b)

=cut

sub i {
  my $self = shift;
  if (defined $self->{'i'}) {
    return $self->{'i'};
  } elsif (defined $self->{'b'}) {
    if ($self->{'a'} == $self->{'b'}) {
      return undef();
    } else {
      return $self->{'a'}/($self->{'a'}-$self->{'b'});
    }
  } elsif (defined $self->{'f'}) {
    return 1/$self->{'f'};
  } else {
    return undef();
  }
}

=head2 invf

Method synonym for the i method

  my $i = $obj->invf; #i=1/f

=cut

sub invf {
  my $self = shift;
  return $self->i(@_);
}

=head2 e

Method returns the value of the first eccentricity, e.  This is the eccentricity of the earth's elliptical cross-section.

  my $e=$obj->e;

=cut

sub e {
  my $self = shift;
  return sqrt($self->e2);
}

=head2 e2

Method returns the value of eccentricity squared (e.g. e^2). This is not the second eccentricity, e' or e-prime see the "ep" method.

  my $e2 = sqrt($obj->e2); #e^2 = f(2-f) = 2f-f^2 = 1-b^2/a^2

=cut

sub e2 {
  my $self = shift;
  my $f    = $self->f();
  return $f*(2 - $f);
}

=head2 ep

Method returns the value of the second eccentricity, e' or e-prime.  The second eccentricity is related to the first eccentricity by the equation: 1=(1-e^2)(1+e'^2).

  my $ep = $obj->ep;

=cut

sub ep {
  my $self = shift;
  return sqrt($self->ep2);
}

=head2 ep2

Method returns the square of value of second eccentricity, e' (e-prime).  This is more useful in almost all equations.

  my $ep=sqrt($obj->ep2);  #ep2=(ea/b)^2=e2/(1-e2)=a^2/b^2-1

=cut

sub ep2 {
  my $self = shift;
  my $a    = $self->a();
  my $b    = $self->b();
  return $a**2/$b**2 - 1;
}

=head2 n

Method returns the value of n given latitude (degrees).  Typically represented by the Greek letter nu, this is the radius of curvature of the ellipsoid perpendicular to the meridian plane.  It is also the distance from the point in question to the polar axis, measured perpendicular to the ellipsoid's surface.

  my $n = $obj->n($lat);

Note: Some define a variable n as (a-b)/(a+b) this is not that variable.

Note: It appears that n can also be calculated as

  n = a^2/sqrt(a^2 * cos($lat)^2 + $b^2 * sin($lat)^2);

=cut

sub n {
  my $self = shift;
  my $lat  = shift; #degrees
  die('Error: Latitude (degrees) required.') unless defined $lat;
  return $self->n_rad(rad_deg($lat));
}

=head2 n_rad

Method returns the value of n given latitude (radians).

  my $n_rad = $obj->n_rad($lat);

Reference: John P. Snyder, "Map Projections: A Working Manual", USGS, page 25, equation (4-20) http://pubs.er.usgs.gov/usgspubs/pp/pp1395

=cut

sub n_rad {
  my $self = shift;
  my $lat  = shift; #radians
  die('Error: Latitude (radians) required.') unless defined $lat;
  my $a    = $self->a;
  my $e2   = $self->e2;
  return $a / sqrt(1 - $e2 * sin($lat)**2);
}

=head2 rho

rho is the radius of curvature of the earth in the meridian plane.

  my $rho=$obj->rho($lat);

=cut

sub rho {
  my $self = shift;
  my $lat  = shift; #degrees
  die('Error: Latitude (degrees) required.') unless defined $lat;
  return $self->rho_rad(rad_deg($lat));
}

=head2 rho_rad

rho is the radius of curvature of the earth in the meridian plane. Sometimes denoted as R'.

  my $rho = $obj->rho_rad($lat);

Reference: John P. Snyder, "Map Projections: A Working Manual", USGS, page 24, equation (4-18) http://pubs.er.usgs.gov/usgspubs/pp/pp1395

=cut

sub rho_rad {
  my $self = shift;
  my $lat  = shift; #radians
  die('Error: Latitude (radians) required.') unless defined $lat;
  my $a    = $self->a;
  my $e2   = $self->e2;
  return $a * (1-$e2) / ( 1 - $e2 * sin($lat)**2 )**(3/2)
  #return $a * (1-$e2) / sqrt(1 - $e2 * sin($lat)**(3/2)); #Bad formula from somewhere
}

=head2 polar_circumference

Method returns the value of the semi-minor axis times 2*PI.

  my $polar_circumference=$obj->polar_circumference;

=cut

sub polar_circumference {
  my $self = shift;
  return 2 * PI() * $self->b();
}

=head2 equatorial_circumference

Method returns the value of the semi-major axis times 2*PI.

  my $equatorial_circumference=$obj->equatorial_circumference;

=cut

sub equatorial_circumference {
  my $self = shift;
  return 2 * PI() * $self->a();
}

sub _setref {
  my $self  = shift;
  my $param = shift;
  if ('HASH' eq ref($param)) {
    if (defined($param->{'a'})) {
      $self->{'a'}          = $param->{'a'};
      $self->{'shortname'}  = 'Custom' unless defined($self->shortname);
      if (defined $param->{'i'}) {
        $self->{'i'}        = $param->{'i'};
        undef($self->{'b'});
        undef($self->{'f'});
        $self->{'longname'} = 'Custom Ellipsoid {a=>'.$self->a.',i=>'.$self->i.'}'  unless defined($self->longname);
      } elsif (defined $param->{'b'}){
        $self->{'b'}        = $param->{'b'};
        undef($self->{'i'});
        undef($self->{'f'});
        $self->{'longname'} = 'Custom Ellipsoid {a=>'.$self->a.',b=>'.$self->b.'}'  unless defined($self->longname);
      } elsif (defined $param->{'f'}){
        $self->{'f'}        = $param->{'f'};
        undef($self->{'b'});
        undef($self->{'i'});
        $self->{'longname'} = 'Custom Ellipsoid {a=>'.$self->a.',f=>'.$self->f.'}'  unless defined($self->longname);
      } else {
        $self->{'b'}        = $param->{'a'};
        undef($self->{'f'});
        undef($self->{'i'});
        $self->{'longname'} = 'Custom Sphere {a=>'.$self->a.'}' unless defined($self->longname);
      }
    } else {
      die('Error: a must be defined');
    }
  } else {
    die('Error: a hash reference e.g. {a=>###, i=>###} must be define');
  }
  return 1;
}

sub _setname {
  my $self  = shift;
  my $param = shift;
  my $ref   = $self->name2ref($param);
  if ('HASH' eq ref($ref)) {
    $self->{'shortname'}=$param;
    my $data            = $self->data;
    my %data            = map {$_, $data->{$_}->{'name'}} (keys %$data);
    $self->{'longname'} = $data{$param};
    return $self->_setref($ref);
  } else {
    die("Error: Ellipsoid $param was not found");
  }
}

=head2 shortname

Method returns the shortname, which is the hash key, of the current ellipsoid

  my $shortname = $obj->shortname;

=cut

sub shortname {
  my $self = shift;
  return $self->{'shortname'};
}

=head2 longname

Method returns the long name of the current ellipsoid

  my $longname = $obj->longname;

=cut

sub longname {
  my $self = shift;
  return $self->{'longname'};
}

=head2 data

Method returns a hash reference for the ellipsoid definition data structure.

  my $datastructure = $obj->data;

=cut

sub data {
#Information from
#  http://earth-info.nga.mil/GandG/coordsys/datums/datumorigins.html
#  http://www.ngs.noaa.gov/PC_PROD/Inv_Fwd/

  return {

    WGS84=>{name=>'World Geodetic System of 1984',
            data=>{a=>6378137,i=>298.257223563},
            alias=>[qw{WGS-84 NAD83 NAD-83}]},

    GRS80=>{name=>'Geodetic Reference System of 1980',
            data=>{a=>6378137,i=>298.25722210088},
            alias=>['GRS-80','GDA','Geocentric Datum of Australia']},

    'Clarke 1866'=>{name=>'Clarke Ellipsoid of 1866',
                    data=>{a=>6378206.4,i=>294.9786982138},
                    alias=>[qw{NAD27 NAD-27}]},

    'Airy 1858'=>{name=>'Airy 1858 Ellipsoid',
                  data=>{a=>6377563.396,i=>299.3249646}},


    'Airy Modified'=>{name=>'Modified Airy Spheroid',
                      data=>{a=>6377340.189,b=>6356034.448}},

    'Australian National'=>{name=>'Australian National Spheroid of 1965',
                            data=>{a=>6378160,i=>298.25},
                            alias=>['Australian 1965']},

    'Bessel 1841'=>{name=>'Bessel 1841 Ellipsoid',
                    data=>{a=>6377397.155,i=>299.1528128}},

    'Clarke 1880'=>{name=>'Clarke Ellipsoid of 1880',
                    data=>{a=>6378249.145,b=>6356514.966}},

    'Clarke 1866'=>{name=>'Clarke Ellipsoid of 1866',
                    data=>{a=>6378206.4,b=>6356583.8}},

    'Danish 1876'=>{name=>'Danish Spheroid of 1876',
                     data=>{a=>3271883.25*1.94903631,i=>300.00}},

    'Everest 1830'=>{name=>'Everest Spheroid of 1830',
                     data=>{a=>6377276.345,i=>300.8017}},

    'Everest Modified'=>{name=>'Modified Everest Spheroid',
                         data=>{a=>6377304.063,i=>300.8017}},

    'Fisher 1960'=>{name=>'Fisher 1960',
                    data=>{a=>6378166,i=>298.3}},

    'Fisher 1968'=>{name=>'Fisher 1968',
                    data=>{a=>6378150,i=>298.3}},

    'Hough 1956'=>{name=>'Hough 1956',
                   data=>{a=>6378270,i=>297}},

    'International (Hayford)'=>{name=>'International - 1924 (Hayford - 1909)',
                                data=>{a=>6378388,i=>297}},

    'Krassovsky 1938'=>{name=>'Krassovsky 1938',
                        data=>{a=>6378245,i=>298.3},
                        alias=>['Krasovsky 1940']},

    'NWL-9D'=>{name=>'NWL-9D Ellipsoid',
               data=>{a=>6378145,i=>298.25},
               alias=>['WGS-66'=>'World Geodetic System 1966']},

    'SA69'=>{name=>'South American 1969',
             data=>{a=>6378160,i=>298.25},
             alias=>['SA-69']},

    'SGS85'=>{name=>'Soviet Geodetic System 1985',
              data=>{a=>6378136,i=>298.257},
              alias=>['SGS-85']},

    'WGS72'=>{name=>'World Geodetic System 1972',
              data=>{a=>6378135,i=>298.26},
              alias=>['WGS-72']},

    'WOS'=>{name=>'War Office Spheroid',
            data=>{a=>6378300.58,i=>296}},

    'UTM'=>{name=>'Department of the Army Universal Transverse Mercator',
            data=>{a=>6378249.2,b=>6356515.0}},
  };
}

=head2 name2ref

Method returns a hash reference (e.g. {a=>6378137,i=>298.257223563}) when passed a valid ellipsoid name (e.g. 'WGS84').

  my $ref=$obj->name2ref('WGS84')

=cut

sub name2ref {
  my $self = shift;
  my $key  = shift;
  my $data = $self->data;
  return $data->{$key}->{'data'};
}

1;

__END__