Convert British and Irish National Grid references to or from WGS84 geodetic coordinates

Question

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 failure\n";

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,$y\n";
    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,$y\n";
    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  $lat\n" : "$lat  $lon\n";
    } else {
        tr/,'"/ ms/;          # ' # for prettify
        s/°/d/g;      # UTF-8 multibyte chars don't work with 'tr'
        s/′/m/g;
        s/″/s/g;
        s/($numpat)m\s*($numpat)s?/($1 + $2\/60.0) . "m"/oeg;
        s/($numpat)d(?:eg)?\s*($numpat)(?:m\s*)?/($1 + $2\/60.0)/oeg;
        tr/d//d;
        s/\s*\b([nsew])\b\s*/$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, $lon\n" if $verbose;
        my $point = ($ns || $ew || $lonlat) ? [$lon, $lat] : [$lat, $lon];
        my ($x, $y) = @{$wgs84->transform($proj, $point)};
        print $xy2grid->($x, $y, $precision), "\n";
    }
}

Answer 1

NB: This review assumes Perl5, specifically the Unicode features in 5.12 (released 2010) and later.

---

1. the parsing could be simpler and more featureful

Much code is devoted to handling delimiters that we're only going to throw away.

Explicit N/S/E/W should override -lonlat but don't.

The sole error message ("malformed input") is vague and happens at the very end, after a series of transformations on the input. The mangled string—which may not resemble original input much anymore—is included in the error message and only adds to the confusion.

In general: modifying an input string to impart meaning is usually a mistake. Modify to remove noise, extract the meaningful parts as structured data, and deal with them there.

2. there is a fair amount of duplicated or nearly-duplicated code

A dispatch table is the standard way to choose code based on data. Your "a2b" functions have a lot of common code, and can be merged once the unique parts are moved into a data structure.

3. the data representations could be more suitable

squares and tosquare use 2-digit values, but you never need values in that format. You always need a pair of single digits, and this complicates the conversion functions. Restructure to suit the need, such that $squares{A} == [ 0, 4 ] (hash of arrays) and $tosquare[0][4] == 'A' (array of arrays).

100000 is better written as 100_000 or 1e5.

$numpat can be simplified to qr/[+-]? \d+ \.?\d* \s* /x. Write regular expression pieces with the qr/REGEXP/ quoting construct, so that they are only compiled once; you then won't need /o modifiers when you reference them. The /x modifier allows the use of whitespace in regular expressions, and makes long expressions more readable. Space within [ ] is still recognized; other whitespace is ignored.

4. Unicode handling is haphazard

This is an artifact of writing in Perl4, which had no Unicode facilities. In Perl5, UTF-8 source code (s/°/d/g; etc.) should inform Perl of the source encoding via use utf8;.

To accept UTF-8 input, STDIN should be placed in :utf8 mode, via binmode STDIN, ":utf8". As you're including user input in die messages, STDERR should get the same treatment.

5. tricks and minor stuff

getnext() is about three times longer and more confusing than it ought to be; see below for a revised version.

Every output ends in a newline; use the -l switch instead.

%tosquare = reverse %squares is the idiomatic version of %tosquare = map { ($squares{$_}, $_) } keys %squares.

local $_ = shift; is usually what you want when assigning to $_ in a sub, else it will be clobbered in the calling scope. (The rewrite contravenes this advice and clobbers $_ on purpose.)

nadgrids= can be adjusted at setup time to ignore missing files. Calls to ->transform() should print error on failure (due to, say, a missing nadgrids file :)

A long series of synonym-to-canonical-value substitutions, as you're doing with s/°/d/g, etc., can be replaced by a hash table where the keys combine into a regex, as in:

    my %decoratives=(qw( ' m   " s   ° d   ′ m   ″ s ),  ("," => " ") );
    s/([@{[ join '', keys %decoratives ]}])/$decoratives{$1}/g; 

revision

Here's my response to my own criticisms. It's not much shorter—about 75% of the original's size—but does improve the error messages and is (perhaps) more clear in its intent.

#!/usr/bin/perl -wl
use strict;
use Getopt::Long;
use Geo::Proj4;
use utf8;
binmode STDIN, ":utf8";
binmode STDERR, ":utf8";
sub grid2xy(_);
sub xy2grid($$$);
sub getnext();

my %squares = qw( 
    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;
$tosquare[ int $squares{$_}/10 ][ $squares{$_}%10 ] = $_ for keys %squares;
$_ = [ split '' ] for values %squares; 

my %howto=(
    GB => {
        setup => 'epsg:27700 +nadgrids=' . join(',' => grep { -f } qw( scotland.gsb england-wales.gsb )),
        parse => qr/^(\D\D)(\d+)$/,
        xy2os => sub {  [ map { int($_[$_]/5e5) + 2 - $_ } 0..1 ], [ map { ($_[$_]/1e5) % 5 } 0..1 ] },
        os2xy => sub { map { 5*$_[0][$_] + $_[1][$_] - 10 + 5*$_ } 0..1 }
    },
    IE => {
        setup => 'epsg:29901',
        parse => qr/^(\D)(\d+)$/,
        xy2os => sub { [ map int($_/1e5) => @_ ] },
        os2xy => sub { @{ $_[0] } }
    }
);

my ($grid, $datum, $precision,$lonlat,$reverse,$verbose) = ('GB', 'WGS84', 5);
GetOptions(
    'grid=s'      => \$grid,
    'reverse!'    => \$reverse,
    'lonlat!'     => \$lonlat,
    'datum=s'     => \$datum,
    'precision=i' => \$precision,
    'verbose!'    => \$verbose
) or die "Option parsing failure\n";

our %do=%{ $howto{$grid} or die "Unknown grid $grid\n" };

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

@ARGV=('-') unless @ARGV;
while (getnext) { 
    if ($reverse) {
        my @lola = @{ $proj->transform($wgs84, grid2xy) or die $proj->error };
        local $,="  ";
        print $lonlat ? @lola : reverse @lola;
    } else {
        my @tokens= map {uc} /( [+-]? \d+ \.?\d* | [NSEW] )/gix;
        print "tokens: @tokens" if $verbose;
        my @lalo=(0,0);
        my @dms=( 1, 60, 3600 );
        my ($unit,$ll,$seenNS, $seenEW)=(0,0,0,0);
        my %seen=( N => \$seenNS, S => \$seenNS, E => \$seenEW, W => \$seenEW );
        my %sign=( N => 1, S => -1, E => 1, W => -1 );
        while (@tokens) { 
            my $tok=shift @tokens;
            if ($sign{$tok}) { 
                die "Repeated or conflicting direction '$tok'\n" if ${ $seen{$tok} };
                die "Directions come after the coordinates\n" unless $unit;
                $lalo[$ll++] *= $sign{$tok};
                ${ $seen{$tok} } = $ll;  # after the increment so that it's nonzero.
                $unit=0;
            } else {
                if ($unit>$#dms) {  $ll++; $unit=0; }
                die "Too many coordinates in '$_'\n" if $ll>1;
                $lalo[$ll] += $tok / $dms[$unit++];
            }
        }
        @lalo=reverse @lalo if (!$seenNS && !$seenEW && $lonlat or $seenNS==1 or $seenEW==2);
        print STDERR "lat/lon @lalo" if $verbose;
        my ($x, $y) = @{ $wgs84->transform($proj, [ @lalo ]) or die $wgs84->error };
        print xy2grid($x, $y, $precision);
    }
}
exit 0;

sub grid2xy(_) {
    local $_=shift;
    my ($letters, $numbers) = /$do{parse}/ or die "Malformed ref '$_'\n";
    my $len = length $numbers;
    die "Odd gridref length - '$_' ($len)\n" if $len % 2;
    $len /= 2;
    my @sq = map { $squares{$_} or die "Non-grid square '$_'\n" } split '', $letters;
    my ($x,$y) = $do{os2xy}(@sq);
    $x = 100000 * ("$x.".substr($numbers, 0, $len).'5');
    $y = 100000 * ("$y.".substr($numbers, $len).'5');
    return [$x, $y];
}

sub xy2grid($$$) { 
   my ($x, $y, $prec) = @_;
   local $,=","; # for the die()
   my $sq = join '', map { $tosquare[ $_->[0] ][ $_->[1] ] or die "No square for @$_\n" } $do{xy2os}($x,$y); 
   return sprintf('%s%s%s', $sq, map { substr(100000 + $_%100000, 1, $prec) } ($x, $y));
}

sub getnext() {
    if (@ARGV and $ARGV[0] eq '-') {
        if ($_ = <STDIN>) { chomp; return $_ }
        else { shift @ARGV }
    }
    return $_=shift @ARGV;
}