view release on metacpan or  search on metacpan

lib/Geo/Coordinates/OSGB.pm  view on Meta::CPAN

    my $af = $a * CONVERGENCE_FACTOR;

    my $dn = $northing - ORIGIN_NORTHING;
    my $de = $easting - ORIGIN_EASTING;

    my $phi = ORIGIN_LATITUDE + $dn/$af;
    my $lam = ORIGIN_LONGITUDE; 

    my ($M, $p_plus, $p_minus);
    while (1) {
        $p_plus  = $phi + ORIGIN_LATITUDE;
        $p_minus = $phi - ORIGIN_LATITUDE;
        $M = $b * CONVERGENCE_FACTOR * (
               (1 + $n * (1 + 5/4*$n*(1 + $n)))*$p_minus
             - 3*$n*(1+$n*(1+7/8*$n))  * sin(  $p_minus) * cos(  $p_plus)
             + (15/8*$n * ($n*(1+$n))) * sin(2*$p_minus) * cos(2*$p_plus)
             - 35/24*$n**3             * sin(3*$p_minus) * cos(3*$p_plus)
               );
        last if abs($dn-$M) < HUNDREDTH_MM;
        $phi = $phi + ($dn-$M)/$af;
    }

    my $cp = cos $phi; 
    my $sp = sin $phi; 
    my $tp  = $sp / $cp; # cos phi cannot be zero in GB

    my $splat = 1 - $e2 * $sp * $sp;
    my $sqrtsplat = sqrt $splat;
    my $nu  = $af / $sqrtsplat;
    my $rho = $af * (1 - $e2) / ( $splat * $sqrtsplat );
    my $eta2 = $nu / $rho - 1;

    my $VII  = $tp / (2 * $rho * $nu);
    my $VIII = $tp / (24 * $rho * $nu**3) * (5 + $eta2 + ( 3 - 9 * $eta2 ) * $tp * $tp );
    my $IX   = $tp / (720 * $rho * $nu**5) * (61 + ( 90 + 45 * $tp * $tp ) * $tp * $tp );

    my $secp = 1/$cp;

    my $X    = $secp / $nu;
    my $XI   = $secp / (    6 * $nu**3 ) * ( $nu / $rho + 2 * $tp * $tp );
    my $XII  = $secp / (  120 * $nu**5 ) * ( 5 + ( 28 + 24 * $tp * $tp ) * $tp * $tp );
    my $XIIA = $secp / ( 5040 * $nu**7 ) * ( 61 + ( 662 + ( 1320 + 720 * $tp * $tp ) * $tp * $tp ) * $tp * $tp );

    $phi = $phi +        ( -$VII + ( $VIII - $IX * $de * $de ) * $de * $de) * $de * $de;
    $lam = $lam + ( $X + ( -$XI + ( $XII - $XIIA * $de * $de ) * $de * $de) * $de * $de) * $de;

    # now put into degrees & return
    return ($phi * 57.29577951308232087679815481410517,
            $lam * 57.29577951308232087679815481410517);
}

1;

=pod

=head1 NAME

Geo::Coordinates::OSGB - Convert coordinates between Lat/Lon and the British National Grid

An implementation of co-ordinate conversion for England, Wales, and Scotland
based on formulae and data published by the Ordnance Survey of Great Britain.

=head1 VERSION

2.20

=for HTML <a href="https://travis-ci.org/thruston/perl-geo-coordinates-osgb">
<img src="https://travis-ci.org/thruston/perl-geo-coordinates-osgb.svg?branch=master"></a>

=head1 SYNOPSIS

  use Geo::Coordinates::OSGB qw(ll_to_grid grid_to_ll);

  ($easting,$northing) = ll_to_grid($lat,$lon);
  ($lat,$lon) = grid_to_ll($easting,$northing);

=head1 DESCRIPTION

These modules convert accurately between OSGB national grid references and
coordinates given in latitude and longitude.

The default "ellipsoid model" used for the conversions is the I<de facto>
international standard WGS84.  This means that you can take latitude and
longitude readings from your GPS receiver, or read them from Wikipedia, or
Google Earth, or your car's sat-nav, and use this module to convert them to
accurate British National grid references for use with one of the Ordnance
Survey's paper maps.  And I<vice versa>, of course.

The module is implemented purely in Perl, and should run on any platform with
Perl version 5.8 or better.

In this description, the abbreviations `OS' and `OSGB' mean `the Ordnance
Survey of Great Britain': the British government agency that produces the
standard maps of England, Wales, and Scotland.  Any mention of `sheets' or
`maps' refers to one or more of the map sheets defined in the accompanying maps
module.

This code is written for the British national grid system.  It is of no use
outside Britain.  In fact it's only really useful in the areas covered by the
OS's main series of maps, which exclude the Channel Islands and Northern
Ireland.

=head1 SUBROUTINES/METHODS

The following functions can be exported from the
C<Geo::Coordinates::OSGB> module:

    grid_to_ll 
    ll_to_grid

Neither of these is exported by default.

=head2 Main subroutines

=head3 C<ll_to_grid(lat, lon)>

C<ll_to_grid> translates a latitude and longitude pair into a grid
easting and northing pair.

When called in a list context, C<ll_to_grid> returns the easting and
northing as a list of two.  When called in a scalar context, it returns

lib/Geo/Coordinates/OSGB.pm  view on Meta::CPAN

conversion algorithm implemented in most consumer-level GPS devices, which 
generally do not have enough memory space for the whole of OSTN.  It
is based on parameters supplied by the OS; they suggest that in most of
the UK this conversion is accurate to within about 5m.

    my ($e, $n) = ll_to_grid_helmert(51.477811, -0.001475);  # RO Greenwich

The input must be decimal degrees in the WGS84 model, with latitude
first and longitude second.  The results are rounded to the nearest
whole metre.  They can be used with C<format_grid> in the same way as
the results from C<ll_to_grid>.

This function is called automatically by C<ll_to_grid> if your
coordinates are WGS84 and lie outside the OSTN polygon.

=head3 C<grid_to_ll_helmert(e,n)>

You can use this function to do a slightly quicker conversion from OS grid
references to WGS84 latitude and longitude coordinates without using the
whole OSTN data set.  The algorithm used is known as a Helmert
transformation.  This is the usual coordinate conversion algorithm
implemented in most consumer-level GPS devices.  It is based on
parameters supplied by the OS; they suggest that in most of the UK this
conversion is accurate to within about 5m.

    my ($lat, $lon) = grid_to_ll_helmert(538885, 177322);

The input must be in metres from false point of origin (as produced by
C<parse_grid>) and the results are in decimal degrees using the WGS84
model.

The results are returned with the full Perl precision in the same way as
C<grid_to_ll> so that you can choose an appropriate rounding for your
needs.  Four or five decimal places is probably appropriate in most
cases.  This represents somewhere between 1 and 10 m on the ground.

This function is called automatically by C<grid_to_ll> if the grid reference
you supply lies outside the OSTN polygon.  (All such spots are far out to sea).
The results are only useful close to mainland Britain.

=head3 Importing all the functions

You can import all the functions defined in C<OSGB.pm> with an C<:all> tag.

    use Geo::Coordinates::OSGB ':all';

=head1 EXAMPLES

  use Geo::Coordinates::OSGB qw/ll_to_grid grid_to_ll/;

  # Latitude and longitude according to the WGS84 model
  ($lat, $lon) = grid_to_ll($e, $n);

  # and to go the other way
  ($e, $n) = ll_to_grid($lat,$lon);

See the test files for more examples of usage.

=head1 BUGS AND LIMITATIONS

The formulae supplied by the OS and used for the conversion routines are
specifically designed to be close floating-point approximations rather
than exact mathematical equivalences.  So after round-trips like these:

  ($lat1,$lon1) = grid_to_ll(ll_to_grid($lat0,$lon0));
  ($e1,$n1)     = ll_to_grid(grid_to_ll($e0,$n0));

neither C<$lat1 == $lat0> nor C<$lon1 == $lon0> nor C<$e1 == $e0> nor
C<$n1 == $n0> exactly.  However the differences should be very small.

The OS formulae were designed to give an accuracy of about 1 mm of
error.  This means that you can rely on the third decimal place for grid
references and about the seventh or eighth for latitude and longitude
(although the OS themselves only provide six decimal places in their
results).

For all of England, Wales, Scotland, and the Isle of Man the error will be
tiny.  All other areas, like Northern Ireland, the Channel Islands or Rockall,
and any areas of sea more than a few miles off shore, are outside the coverage
of OSTN, so the simpler, less accurate transformation is used.  The OS state
that this is accurate to about 5m but that the parameters used are only valid
in the reasonably close vicinity of the British Isles.

Not enough testing has been done.  I am always grateful for the feedback
I get from users, but especially for problem reports that help me to
make this a better module.

=head1 DIAGNOSTICS

The only error message you will get from this module is about the
ellipsoid shape used for the transformation.  If you try to set C<<
{shape => 'blah'} >> the module will croak with a message saying
C<Unknown shape: blah>.  The shape should be one of the shapes defined:
WGS84 or OSGB36.

Should this software not do what you expect, then please first read this
documentation, secondly verify that you have installed it correctly and
that it passes all the installation tests on your set up, thirdly study
the source code to see what it's supposed to be doing, fourthly get in
touch to ask me about it.

=head1 CONFIGURATION AND ENVIRONMENT

There is no configuration required either of these modules or your
environment.  It should work on any recent version of Perl, on any
platform.

=head1 DEPENDENCIES

Perl 5.10 or better.

=head1 INCOMPATIBILITIES

None known.

=head1 LICENSE AND COPYRIGHT

Copyright (C) 2002-2017 Toby Thurston

OSTN transformation data included in this module is freely available
from the Ordnance Survey but remains Crown Copyright (C) 2002

This program is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.

This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
General Public License for more details.