Nsryjdtyk

I've been using this wgs84togrid program for a few years. It converts in both directions between National Grid coordinates for GB or Ireland (beginning with a letter or letter-pair identifying the 100km square) and latitude/longitude positions (in decimal degrees, decimal minutes or decimal seconds) on a WGS84 ellipsoid.

It acts as a filter, expecting one point per line, copying unchanged any unrecognised parts of the line.

Program options (all may be shortened, provided that's unambiguous):

• 
  -grid: choose a grid: GB (default) or IE


• 
  -reverse: reverse direction - convert National Grid positions to WGS84


• 
  -lonlat: geodesic positions are longitude first


• 
  -datum: use alternative datum instead of WGS84 (National Grid coordinates are always on the appropriate fixed datum)


• 
  -precision: how many digits to include in northings and eastings (default: 5, which gives 1-metre resolution)


• 
  -verbose: extra output (to confirm that lat/lon are parsed as you expect).

Example use (in Bash):

$ wgs84togrid -p 3 <<<"55°56′55″N 3°12′03″W"

NT251734

$ wgs84togrid -r <<<NT251734

55.9482278708547  -3.20011121889597

The heavy work of coordinate transformation is performed by the PROJ.4 library; all I do is manage the grid letters and I/O formats.

I assume the presence of scotland.gsb and england-wales.gsb grid corrections files, but that option may be removed if you don't have them, at the cost of a few metres of accuracy (< 10m, I'm sure).

Specifically out of scope:

• I don't check that the point is within the valid area of the chosen grid (and certainly don't think of auto-selecting the correct grid).


• No plans to support any other grids elsewhere in the world.

#!/usr/bin/perl -w



use strict;



use Getopt::Long;



use Geo::Proj4;







my %squares = (A=>'04', B=>'14', C=>'24', D=>'34', E=>'44',

               F=>'03', G=>'13', H=>'23', J=>'33', K=>'43',

               L=>'02', M=>'12', N=>'22', O=>'32', P=>'42',

               Q=>'01', R=>'11', S=>'21', T=>'31', U=>'41',

               V=>'00', W=>'10', X=>'20', Y=>'30', Z=>'40');



my %tosquare = map { ($squares{$_}, $_) } keys %squares;



my $grid = 'GB';

my $lonlat;

my $datum = 'WGS84';

my $precision = 5;

my $reverse;

my $verbose;





GetOptions('grid=s'      => $grid,

       'reverse!'    => $reverse,

           'lonlat!'     => $lonlat,

       'datum=s'     => $datum,

           'precision=i' => $precision,

       'verbose!'    => $verbose) or die "Option parsing failuren";



sub any2xy($$$) {

    my ($x, $y, $numbers) = @_;

    my $len = length $numbers;

    die "Odd gridref length - '$_' ($len)n" if $len % 2;

    $len /= 2;

    $x = 100000 * ("$x.".substr($numbers, 0, $len).'5');

    $y = 100000 * ("$y.".substr($numbers, $len).'5');

    return [$x, $y];

}



sub osgb2xy($) {

    local $_ = shift;

    my ($letters, $numbers) = m/^(D{2})(d+)$/ or die "Malformed OSGB ref '$_'n";

    my $x = 0;

    my $y = 0;

    foreach (split '', $letters) {

        my @sq = split '', $squares{$_} or die "Non-grid square '$_'n";

        $x = 5 * $x + $sq[0];

        $y = 5 * $y + $sq[1];

    }

    $x -= 10;

    $y -=  5;

    return any2xy($x, $y, $numbers);

}



sub osi2xy($) {

    $_ = shift;

    my ($letters, $numbers) = m/^(D)(d+)$/ or die "Malformed OSI ref '$_'n";

    my ($x, $y) = split '', $squares{$letters} or die "Non-grid square '$_'n";

    return any2xy($x, $y, $numbers);

}



sub togrid($$$$) {

    my ($sq, $x, $y, $prec) = @_;

    return sprintf('%s%s%s', $sq, map { substr(100000 + $_%100000, 1, $prec) } ($x, $y));

}



sub xy2osi($$$) {

    my ($x, $y, $prec) = @_;

    my $sq = $tosquare{int($x/100000) . int($y/100000)} or die "No square for $x,$yn";

    return togrid($sq, $x, $y, $prec);

}



sub xy2osgb($$$) {

    my ($x, $y, $prec) = @_;

    $x += 1000000;

    $y +=  500000;

    my $sq = $tosquare{int($x/500000) . int($y/500000)} . $tosquare{int($x/100000)%5 . int($y/100000)%5} or die "No square for $x,$yn";

    return togrid($sq, $x, $y, $prec);

}



my $inputs;

sub getnext();

sub getnext() {

    if ($inputs) {

        $_ = <$inputs>;

        return $_ if $_;

        $inputs = undef;

    }

    if (@ARGV) {

        $_ = shift @ARGV;

        if ($_ eq '-') {

            $inputs = *STDIN;

            return getnext();

        }

        return $_;

    }

    return undef;

}





my $wgs84 = Geo::Proj4->new(proj => 'latlon', datum => $datum) or die Geo::Proj4->error;



my ($proj, $xy2grid, $grid2xy);

if (uc $grid eq 'GB') {

    $proj = Geo::Proj4->new(init => 'epsg:27700 +nadgrids=scotland.gsb,england-wales.gsb') or die Geo::Proj4->error;

    $xy2grid = &xy2osgb;

    $grid2xy = &osgb2xy;

} elsif (uc $grid eq 'IE') {

    $proj = Geo::Proj4->new(init => 'epsg:29901') or die Geo::Proj4->error;

    $xy2grid = &xy2osi;

    $grid2xy = &osi2xy;

} else {

    die "Unknown grid '$grid'n";

}



my $numpat = '[+-]?d+(?:.d+)?s*';



@ARGV=('-') unless @ARGV;

while ($_ = getnext()) {

    chomp;

    if ($reverse) {

        my $point = $grid2xy->($_);

        my ($lon, $lat) = @{$proj->transform($wgs84, $point)};

        print $lonlat ? "$lon  $latn" : "$lat  $lonn";

    } else {

        tr/,'"/ ms/;

        s/°/d/g;

        s/′/m/g;

        s/″/s/g;

        s/($numpat)ms*($numpat)s?/($1 + $2/60.0) . "m"/oeg;

        s/($numpat)d(?:eg)?s*($numpat)(?:ms*)?/($1 + $2/60.0)/oeg;

        tr/d//d;

        s/s*b([nsew])bs*/$1/i;

        tr!/,!  !;

        s/($numpat[ew ])s*($numpat[ns]?)/$2 $1/oi;

        s/($numpat)s+($numpat[ns]|[ns]$numpat)/$2 $1/oi;



        my ($lat, $ns, $lon, $ew) = m/^s*($numpat)([ns ]?)s*($numpat)([ew]?)s*$/i

            or die "Malformed input: $_n";

        $lat *= -1 if lc $ns eq 's';

        $lon *= -1 if lc $ew eq 'w';

        print STDERR "$lat, $lonn" if $verbose;

        my $point = ($ns || $ew || $lonlat) ? [$lon, $lat] : [$lat, $lon];

        my ($x, $y) = @{$wgs84->transform($proj, $point)};

        print $xy2grid->($x, $y, $precision), "n";

    }

}