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";
}
}
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