10
\$\begingroup\$

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";
    }
}
\$\endgroup\$
  • \$\begingroup\$ You can probably tell that I learnt Perl 4 a couple of decades ago and haven't kept up. Perl 5 suggestions are of course welcome, though you might need to assume much lower understanding from me if that's what you choose to contribute. Thanks! \$\endgroup\$ – Toby Speight Nov 22 '18 at 21:32
  • \$\begingroup\$ Something seems to have got mangled in the GetOptions() call when I pasted this (as if tabs had been crushed to 4 spaces - but I didn't use tabs for indentation!). For the record, the options do all line up, and it's just a SE display bug there (assuming I'm not the only one seeing that!). \$\endgroup\$ – Toby Speight Nov 22 '18 at 21:35
  • \$\begingroup\$ I see the wrong indentation as well. However, in edit mode, the lines line up correctly. \$\endgroup\$ – Martin R Nov 23 '18 at 9:44
  • \$\begingroup\$ Thanks @Martin - I've tried everything I can think of, but no difference. The really annoying thing is that it looks fine in preview. If I can be bothered, perhaps I might ask about it on Meta.SE. \$\endgroup\$ – Toby Speight Nov 23 '18 at 9:56
  • \$\begingroup\$ Have you considered turning those long numbers into constants? \$\endgroup\$ – yuri Mar 10 at 13:03
4
\$\begingroup\$

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;
}
\$\endgroup\$
  • 1
    \$\begingroup\$ Excellent review; I've learnt a lot from this. Thank you! \$\endgroup\$ – Toby Speight Mar 11 at 8:47
  • 1
    \$\begingroup\$ I will be keeping a 5-line version of %squares - the newlines there do actually show the 5x5 structure embodied in the coordinates. \$\endgroup\$ – Toby Speight Mar 11 at 12:31
  • \$\begingroup\$ I think I remember why the error message includes the part-processed (aka mangled) input: I thought it would be useful to see what had been matched. Of course, that makes sense only to the author (I probably added that when debugging) and not to other users, and I need to learn to improve my users-eye-view of what I write! \$\endgroup\$ – Toby Speight Mar 13 at 8:18
  • \$\begingroup\$ yes, that kind of telemetry comes with the territory of parse-by-s///, which is a trap we've all fallen into, because it's seductively immediate and concise. That sentence in my review that begins "In general: …" was multiple, rambling paragraphs on this topic; they didn't make the final draft. \$\endgroup\$ – Oh My Goodness Mar 13 at 10:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.