Geo::Ellipsoid module

Geo::Ellipsoid module

am 21.10.2005 22:30:06 von Jim Gibson

I have written a module I am tentatively calling Geo::Ellipsoid. The
module performs geographical calculations involving location, range,
and bearing on the surface of an ellipsoid, such as the Earth. I am
preparing the module for submission to CPAN, after I am satisfied with
the testing procedures.

There is already a Geo::Distance module in CPAN, but it uses
approximate calculations using one of several different methods. This
module uses a slower but more exact calculational method.

The module is strictly object-oriented and pure Perl.


Here is the POD. I invite your comments:


NAME
Geo::Ellipsoid - Calculate positions, distances, and bearings on the
surface of an ellipsoid.

VERSION
Version 1.0, released October, 2005.

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, 45.0, 2000 );
( $x, $y ) = $geo->displacement( @origin, $lat, $lon );
@pos = $geo->location( $lat, $lon, $x, $y );

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 "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 "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 "DEFINED ELLIPSOIDS" below for the ellipsoid survey values that
may be selected for use by Geo::Ellipsoid.

CONSTRUCTOR
new
The new() constructor may be called with a hash list to set the
value of the ellipsoid to be used and the value of the units to
be used for angles. The default constructor is equivalent to the
following:

my $geo = Geo::Ellipsoid->new(
ellipsoid => 'WGS84',
units => 'radians'
);

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', ell => 'NAD27' ) is
valid. The constructor returns a blessed object of the class
Geo::Ellipsoid.

METHODS
set_units
Set the 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. An invalid value will be accepted as 'radians'.

$geo->set_units('degrees');

set_ellipsoid
Set the ellipsoid to be used by the Geo::Ellipsoid object. See
"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');

set_custom_ellipsoid
Sets the ellipsoid parameters to the specified ( major semiaxis
and reciprocal flattening.

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

scales
Returns a list consisting of the meters 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);

range
Returns the range in meters between two specified locations
given as latitude, longitude pairs.

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

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

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

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

displacement
Returns the (x,y) displacement in meters 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.

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

DEFINED ELLIPSOIDS
The following ellipsoids are defined in Geo::Ellipsoid, with the
semi-major axis in meters and the reciprocal flattening as
shown. The default ellipsoid is WGS84.

Ellipsoid Semi-Major Axis (m.) 1/Flattening
--------- ------------------- ---------------
WGS84 6378137.0 298.25722210088
NAD27 6378206.4 294.9786982138
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 290.8017
EVEREST-MODIFIED 6377304.063 290.8017
FISHER-1960 6378166.0 298.3
FISHER-1968 6378150.0 298.3
HOUGH-1956 6378270.0 297.0
HAYFORD 6378388.0 297.0
KRASSOVSKY-1938 6378245.0 298.3
NWL-9D 6378145.0 298.25
SOUTHAMERICAN-1969 6378160.0 298.25
SOVIET-1985 6378136.0 298.257
WGS72 6378135.0 298.26

LIMITATIONS
The methods should not be used on points which are too near the
poles (above or below 89 degrees), and should not be used on
points which are antipodal, i.e., exactly on opposite sides of
the geoid. The methods will not return valid results in these
cases.

ACKNOWLEDGEMENTS
The conversion algorithms used here are Perl translations of
Fortran routines written by LCDR L. Pfeifer NGS Rockville MD
that implement T. Vincenty's Modified Rainsford's method with
Helmert's elliptical terms as published in "Direct and Inverse
Solutions of Ellipsoid on the Ellipsoid with Application of
Nested Equations", T. Vincenty, Survey Review, April 1975.

The Fortran source code files inverse.for and forward.for may be
obtained from

ftp://ftp.ngs.noaa.gov/pub/pcsoft/for_inv.3d/source/

BUGS
This version cannot handle points that are too near to the
poles, for some of the methods, or points which are anti-podal,
that is on opposite sides of the earth. In this case, the
iterative algorithms will not converge, a warning message will
be emitted, and undefined values will be returned.

--
Jim Gibson

Re: Geo::Ellipsoid module

am 23.10.2005 05:50:28 von unknown

Jim Gibson wrote:
>
> DEFINED ELLIPSOIDS
> The following ellipsoids are defined in Geo::Ellipsoid, with the
> semi-major axis in meters and the reciprocal flattening as
> shown. The default ellipsoid is WGS84.
>
> Ellipsoid Semi-Major Axis (m.) 1/Flattening
> --------- ------------------- ---------------
> WGS84 6378137.0 298.25722210088
> NAD27 6378206.4 294.9786982138
> 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 290.8017
> EVEREST-MODIFIED 6377304.063 290.8017
> FISHER-1960 6378166.0 298.3
> FISHER-1968 6378150.0 298.3
> HOUGH-1956 6378270.0 297.0
> HAYFORD 6378388.0 297.0
> KRASSOVSKY-1938 6378245.0 298.3
> NWL-9D 6378145.0 298.25
> SOUTHAMERICAN-1969 6378160.0 298.25
> SOVIET-1985 6378136.0 298.257
> WGS72 6378135.0 298.26
>

How about GRS80? Or IAU76 if you're interested in what the astronomers
are doing. Apparently the U.S. Department of commerce has its own
variant of GRS80, with the same semimajor, but more decimal places than
the "standard" standard. This per
http://biology.usgs.gov/fgdc.metadata/version2/spref/horiz/g eodet/faq.htm
Oddly, I can't figure out where to get the International Astronomical
Union's ellipsoid from their web site. Buried under some subcommittee, I
guess. I got it from Jean Meeus' "Astronomical Algorithms", but

Letting the user define his/her own reference ellipsoid is a good idea,
also. Can it be done at the class level as well as the object leve, so
that (e.g.) "sphere" becomes, at least for the life of the Perl script,
a defined ellipsoid?

Could you define/document the behavior of set_ellipsoid() if the given
ellipsoid is invalid?

But these are really all nit-picks. I would go for it.

Thanks for the opportunity to comment,
Tom Wyant

Re: Geo::Ellipsoid module

am 28.10.2005 21:30:23 von Jim Gibson

In article , harryfmudd [AT]
comcast [DOT] net wrote:

> Jim Gibson wrote:
> >
> > DEFINED ELLIPSOIDS
> > The following ellipsoids are defined in Geo::Ellipsoid, with the
> > semi-major axis in meters and the reciprocal flattening as
> > shown. The default ellipsoid is WGS84.
> >
> > Ellipsoid Semi-Major Axis (m.) 1/Flattening
> > --------- ------------------- ---------------
> > WGS84 6378137.0 298.25722210088

[other ellipsoids snipped]

> >
>
> How about GRS80? Or IAU76 if you're interested in what the astronomers
> are doing. Apparently the U.S. Department of commerce has its own
> variant of GRS80, with the same semimajor, but more decimal places than
> the "standard" standard. This per
> http://biology.usgs.gov/fgdc.metadata/version2/spref/horiz/g eodet/faq.htm
> Oddly, I can't figure out where to get the International Astronomical
> Union's ellipsoid from their web site. Buried under some subcommittee, I
> guess. I got it from Jean Meeus' "Astronomical Algorithms", but

Thanks for the suggestions.

As far as I can tell, the GRS80 ellipsoid is the same as the WGS84
ellipsoid, and I have added it as a separate option with the same
values as WGS84.

As for IAU76, the best Google would yield are the values (r=6378140,
1/f=298.257) that I got from:



Are those values the same as Meeus? I have added these as 'IAU76'.

>
> Letting the user define his/her own reference ellipsoid is a good idea,
> also. Can it be done at the class level as well as the object leve, so
> that (e.g.) "sphere" becomes, at least for the life of the Perl script,
> a defined ellipsoid?

Specifying a custom ellipsoid not already defined adds that definition
to the class hash of ellipsoids. I have also added the functionality
that specifying an ellipsoid changes the default, so that in

my $e1 = Geo::Ellipsoid->new( ellipsoid => 'IAU76' );
my $e2 = Geo::Ellipsoid->new();

$e2 will have an ellipsoid 'IAU76'.

>
> Could you define/document the behavior of set_ellipsoid() if the given
> ellipsoid is invalid?

If you specify an ellipsoid that is not in the list, it defaults to
'WGS84'. If you specify a custom ellipsoid, any value for semi-major
radius and reciprocal flattening will be accepted. The only check is
for undefined or zero reciprocal flattening, which results in a sphere.

Posted Via Usenet.com Premium Usenet Newsgroup Services
----------------------------------------------------------
** SPEED ** RETENTION ** COMPLETION ** ANONYMITY **
----------------------------------------------------------
http://www.usenet.com

Re: Geo::Ellipsoid module

am 29.10.2005 15:46:40 von unknown

Jim Gibson wrote:

[snip]
>
> As far as I can tell, the GRS80 ellipsoid is the same as the WGS84
> ellipsoid, and I have added it as a separate option with the same
> values as WGS84.

Yeah. There appear to be a couple cases where different standards use
the same values. Duplication is not strictly necessary, but may help
fend off questions from people like me.

>
> As for IAU76, the best Google would yield are the values (r=6378140,
> 1/f=298.257) that I got from:
>
>
>
> Are those values the same as Meeus? I have added these as 'IAU76'.
>

Yes. Jean Meeus' "Astronomical Algorithms," 2nd edition, gives a =
6378.14 km, f = 1/298.257, but gives no reference.

[snip]

>
> Specifying a custom ellipsoid not already defined adds that definition
> to the class hash of ellipsoids. I have also added the functionality
> that specifying an ellipsoid changes the default, so that in
>
> my $e1 = Geo::Ellipsoid->new( ellipsoid => 'IAU76' );
> my $e2 = Geo::Ellipsoid->new();
>
> $e2 will have an ellipsoid 'IAU76'.
>

Interesting. I would probably have had a separate method to change the
default, because it would not have occurred to me to do it this way. But
this approach certainly seems to embody Larry Wall's first
characteristic of a good programmer (the full list is "laziness,
impatience, and hubris").

>
>>Could you define/document the behavior of set_ellipsoid() if the given
>>ellipsoid is invalid?
>
>
> If you specify an ellipsoid that is not in the list, it defaults to
> 'WGS84'. If you specify a custom ellipsoid, any value for semi-major
> radius and reciprocal flattening will be accepted. The only check is
> for undefined or zero reciprocal flattening, which results in a sphere.
>

It might be a Good Thing if the user could find out somehow that he or
she didn't get the ellipsoid specified, just to guard against typos. I
tend to croak for things like this, but it looks like you would prefer
not to. Other possibilities I can think of include:

carp

fail to instantiate, and provide a way to figure out what went wrong.

provide an accessor, so the user can find out which ellipsoid is
actually in use.

You may well come up with something better.


Thanks,
Tom Wyant

Re: Geo::Ellipsoid module

am 31.10.2005 19:07:43 von Jim Gibson

In article <4IadnXvEpq7c4P7enZ2dnUVZ_sydnZ2d@comcast.com>, harryfmudd
[AT] comcast [DOT] net wrote:

> Jim Gibson wrote:


>
> It might be a Good Thing if the user could find out somehow that he or
> she didn't get the ellipsoid specified, just to guard against typos. I
> tend to croak for things like this, but it looks like you would prefer
> not to. Other possibilities I can think of include:
>
> carp
>
> fail to instantiate, and provide a way to figure out what went wrong.
>
> provide an accessor, so the user can find out which ellipsoid is
> actually in use.

Indeed. I did not mention it, but the module uses carp to issue the
following message:

"Ellipsoid $ellipsoid does not exist; defaulting to WGS84"

and similarly when the user specifies a zero flattening reciprocal:

"Infinite flattening specified by ellipsoid -- assuming a sphere"

and for specifying angle units other than radians or degrees:

"Invalid units specifier '$units' - using radians"

Posted Via Usenet.com Premium Usenet Newsgroup Services
----------------------------------------------------------
** SPEED ** RETENTION ** COMPLETION ** ANONYMITY **
----------------------------------------------------------
http://www.usenet.com

Re: Geo::Ellipsoid module

am 02.11.2005 10:41:31 von Ilya Zakharevich

[A complimentary Cc of this posting was sent to
harryfmudd [AT] comcast [DOT] net
<"harryfmudd [AT] comcast [DOT] net">], who wrote in article <4IadnXvEpq7c4P7enZ2dnUVZ_sydnZ2d@comcast.com>:

> > my $e1 = Geo::Ellipsoid->new( ellipsoid => 'IAU76' );
> > my $e2 = Geo::Ellipsoid->new();
> >
> > $e2 will have an ellipsoid 'IAU76'.

> Interesting. I would probably have had a separate method to change the
> default, because it would not have occurred to me to do it this way. But
> this approach certainly seems to embody Larry Wall's first
> characteristic of a good programmer (the full list is "laziness,
> impatience, and hubris").

Well, he (obviously) did not mean designing a module API at that
moment. Having side effects of a ->new() call is IMO absolutely disgusting.

> > If you specify an ellipsoid that is not in the list, it defaults to
> > 'WGS84'. If you specify a custom ellipsoid, any value for semi-major
> > radius and reciprocal flattening will be accepted. The only check is
> > for undefined or zero reciprocal flattening, which results in a sphere.

> It might be a Good Thing if the user could find out somehow that he or
> she didn't get the ellipsoid specified, just to guard against typos.

Data-related modules should never try to correct wrong input.
Guesstimates should be done on much higher level, like in UI. Unknown
input should better result in a croak().

Hope this helps,
Ilya