File Coverage

blib/lib/Geo/Coordinates/OSGB/Grid.pm

Criterion	Covered	Total	%
statement	179	192	93.2
branch	75	96	78.1
condition	23	39	58.9
subroutine	26	26	100.0
pod	7	7	100.0
total	310	360	86.1

package Geo::Coordinates::OSGB::Grid;

use Geo::Coordinates::OSGB::Maps qw{%maps %name_for_map_series};
use Geo::Coordinates::OSGB qw{get_ostn02_shift_pair};

use base qw(Exporter);
use strict;
use warnings;
use Carp;
use 5.008; # At least Perl 5.08 please, be sure to change POD below if you update

our $VERSION = '2.18';

our %EXPORT_TAGS = (all => [qw(
        parse_grid
        format_grid

        parse_trad_grid
        parse_GPS_grid
        parse_landranger_grid
        parse_map_grid

        format_grid_trad
        format_grid_GPS
        format_grid_landranger
        format_grid_map

        random_grid
        )]);

our @EXPORT_OK = ( @{ $EXPORT_TAGS{all} } );

use constant GRID_SQ_LETTERS => 'VWXYZQRSTULMNOPFGHJKABCDE';
use constant GRID_SIZE => sqrt length GRID_SQ_LETTERS;
use constant MINOR_GRID_SQ_SIZE => 100_000;
use constant MAJOR_GRID_SQ_SIZE => GRID_SIZE * MINOR_GRID_SQ_SIZE;
use constant MAJOR_GRID_SQ_EASTING_OFFSET  => 2 * MAJOR_GRID_SQ_SIZE;
use constant MAJOR_GRID_SQ_NORTHING_OFFSET => 1 * MAJOR_GRID_SQ_SIZE;
use constant MAX_GRID_SIZE => MINOR_GRID_SQ_SIZE * length GRID_SQ_LETTERS;

# Produce a random GR
# A simple approach would pick 0 < E < 700000 and 0 < N < 1250000 but that way
# many GRs produced would be in the sea, so pick a random map, and then find a
# random GR within its bbox and finally check that the resulting pair is
# inside the OSTN02 boundary and actually on some map.

sub random_grid {
    my @preferred_sheets = @_;
    # assume plain sheet numbers are Landranger
    for my $s (@preferred_sheets) {
        if ($s =~ m{\A\d+\Z}xsmio) {
            $s =~ s/\A/A:/xsmio
        }
    }
    my @sheets;
    if (@preferred_sheets) {
        @sheets = grep { exists $maps{$_} } @preferred_sheets;
    }
    else {
        @sheets = keys %maps;
    }
    my $margin = 5000;
    my ($map, $lle, $lln, $ure, $urn);
    my ($easting, $northing);
    while (1) {
        $map = $maps{$sheets[int rand @sheets]};
        ($lle, $lln) = @{$map->{bbox}->[0]};
        ($ure, $urn) = @{$map->{bbox}->[1]};
        $easting  = sprintf '%.3f', $lle + rand ($ure-$lle);
        $northing = sprintf '%.3f', $lln + rand ($urn-$lln);
        # check we are well inside OSTN02 and actually on the chosen map
        last if get_ostn02_shift_pair($easting, $northing)
             && get_ostn02_shift_pair($easting, $northing+$margin)
             && get_ostn02_shift_pair($easting, $northing-$margin)
             && get_ostn02_shift_pair($easting+$margin, $northing)
             && get_ostn02_shift_pair($easting-$margin, $northing)
             && 0 != _winding_number($easting, $northing, $map->{polygon});
    }
    return ($easting, $northing);
}

sub format_grid {
    my ($easting, $northing, $options) = @_;

    my $form      = exists $options->{form}   ? uc $options->{form}   : 'SS EEE NNN';
    my $with_maps = exists $options->{maps}   ?    $options->{maps}   : 0;
    my $map_keys  = exists $options->{series} ? uc $options->{series} : join '', sort keys %name_for_map_series;

    my $sq = _grid_to_sq($easting,$northing);
    if ( !$sq ) {
        croak "Too far off the grid: $easting $northing\n";
    }

    my ($e,$n) = map { int } map { $_ % MINOR_GRID_SQ_SIZE } ($easting, $northing);

    my @sheets = ();
    if ( $with_maps ) {
        while (my ($k,$m) = each %maps) {
            next if index($map_keys, substr($k,0,1)) == -1;
            if ($m->{bbox}->[0][0] <= $easting  && $easting  < $m->{bbox}->[1][0]
             && $m->{bbox}->[0][1] <= $northing && $northing < $m->{bbox}->[1][1]) {
                my $w = _winding_number($easting, $northing, $m->{polygon});
                if ($w != 0) {
                    push @sheets, $k;
                }
            }
        }
        @sheets = sort @sheets;
    }

    # special cases
    if ( $form eq 'TRAD' ) {
        $form = 'SS EEE NNN';
    }
    elsif ( $form eq 'GPS' ) {
        $form = 'SS EEEEE NNNNN';
    }
    elsif ( $form eq 'SS' ) {
        return $sq;
    }

    if ( my ($space_a, $e_spec, $space_b, $n_spec) = $form =~ m{ \A S{1,2}(\s*)(E{1,5})(\s*)(N{1,5}) \Z }iosxm ) {
        my $e_len = length $e_spec;
        my $n_len = length $n_spec;
        $e = int($e / 10**(5 - $e_len));
        $n = int($n / 10**(5 - $n_len));

        if ( wantarray ) {
            return ($sq, $e, $n, @sheets)
        }

        my $gr = sprintf '%s%s%0.*d%s%0.*d', $sq, $space_a, $e_len, $e, $space_b, $n_len, $n;
        $gr =~ s/\s+/ /g;

        if ( $with_maps ) {
            if ( @sheets ) {
                return sprintf '%s on %s', $gr, join ', ', @sheets;
            }
            else {
                return sprintf '%s is not on any maps in series %s', $gr, $map_keys;
            }
        }
        else {
            return $gr;
        }
    }

    croak "Format $form was not recognized";
}

sub format_grid_trad {
    my ($easting, $northing) = @_;
    return format_grid($easting, $northing, { form => 'SS EEE NNN' })
}

sub format_grid_GPS {
    my ($easting, $northing) = @_;
    return format_grid($easting, $northing, { form => 'SS EEEEE NNNNN' })
}

sub format_grid_map {
    my ($easting, $northing, $options) = @_;
    if ( defined $options ) {
        $options->{maps} = 1
    }
    else {
        $options = { maps => 1 }
    }
    return format_grid($easting, $northing, $options)
}

sub format_grid_landranger {
    my ($easting,$northing) = @_;

    my ($sq, $e, $n, @sheets) = format_grid($easting, $northing, { form => 'SS EEE NNN', maps => 1, series => 'A' });

    for my $s (@sheets) {
        $s =~ s/\AA://xsm;
    }

    return ($sq, $e, $n, @sheets) if wantarray;

    if (!@sheets )    { return sprintf '%s %03d %03d is not on any Landranger sheet', $sq, $e, $n }
    if ( @sheets==1 ) { return sprintf '%s %03d %03d on Landranger sheet %s'        , $sq, $e, $n, $sheets[0] }
    if ( @sheets==2 ) { return sprintf '%s %03d %03d on Landranger sheets %s and %s', $sq, $e, $n, @sheets }

    my $phrase = join ', ', @sheets[0..($#sheets-1)], "and $sheets[-1]";
    return sprintf '%s %03d %03d on Landranger sheets %s', $sq, $e, $n, $phrase;

}

sub _grid_to_sq {
    my ($e, $n) = @_;

    $e += MAJOR_GRID_SQ_EASTING_OFFSET;
    $n += MAJOR_GRID_SQ_NORTHING_OFFSET;
    return if !(0 <= $e && $e < MAX_GRID_SIZE && 0 <= $n && $n < MAX_GRID_SIZE);

    my $major_index = int $e / MAJOR_GRID_SQ_SIZE + GRID_SIZE * int $n / MAJOR_GRID_SQ_SIZE;
    $e = $e % MAJOR_GRID_SQ_SIZE;
    $n = $n % MAJOR_GRID_SQ_SIZE;
    my $minor_index = int $e / MINOR_GRID_SQ_SIZE + GRID_SIZE * int $n / MINOR_GRID_SQ_SIZE;
    return
       substr(GRID_SQ_LETTERS, $major_index, 1) .
       substr(GRID_SQ_LETTERS, $minor_index, 1);
}

sub _get_grid_square_offsets {
    my $s = shift;
    return if length $s < 2;

    my $a = index GRID_SQ_LETTERS, uc substr $s, 0, 1;
    return if 0 > $a;

    my $b = index GRID_SQ_LETTERS, uc substr $s, 1, 1;
    return if 0 > $b;

    my ($X, $Y) = ($a % GRID_SIZE, int $a / GRID_SIZE);
    my ($x, $y) = ($b % GRID_SIZE, int $b / GRID_SIZE);

    return (
        MAJOR_GRID_SQ_SIZE * $X - MAJOR_GRID_SQ_EASTING_OFFSET  + MINOR_GRID_SQ_SIZE * $x,
        MAJOR_GRID_SQ_SIZE * $Y - MAJOR_GRID_SQ_NORTHING_OFFSET + MINOR_GRID_SQ_SIZE * $y
    );
}

sub _get_eastnorthings {
    my $s = shift;
    my $numbers = $s;
    $numbers =~ tr/0-9//cd; # avoid using "r" here as it requires perl >= 5.14
    my $len = length $numbers;
    croak 'No easting or northing found' if $len == 0;
    croak "Easting and northing have different lengths in $s" if $len % 2;
    croak "Too many digits in $s" if $len > 10;

    # this trick lets us pad with zeros on the right
    my $e = reverse sprintf '%05d', scalar reverse substr $numbers, 0, $len/2;
    my $n = reverse sprintf '%05d', scalar reverse substr $numbers,    $len/2;
    return ($e, $n)
}

sub parse_grid {

    my $options = 'HASH' eq ref $_[-1] ? pop @_ : { };

    my $figs = exists $options->{figs} ? $options->{figs} : 3;

    my @out;

    my $s = @_ < 3 ? "@_" : sprintf '%s %0.*d %0.*d', $_[0], $figs, $_[1], $figs, $_[2];

    # normal case : TQ 123 456 etc
    if ( my ($E, $N) = _get_grid_square_offsets($s) ) {
        my ($e, $n) = length $s > 2 ? _get_eastnorthings(substr $s, 2) : (0,0);
        @out = ($E+$e, $N+$n);
        return wantarray ? @out : "@out";
    }

    # sheet id instead of grid sq
    my $sheet_ref_pattern = qr{\A([A-Z]:)?([0-9NEWSOL/]+?)(\.[a-z]+)?(?:[ -/.]([ 0-9]+))?\Z};
    my ($prefix, $sheet, $suffix, $numbers) = $s =~ m/$sheet_ref_pattern/sioxm;

    if (defined $sheet) {

        if ( ! defined $prefix ) {
            $prefix = 'A:';
        }
        if ( ! defined $suffix ) {
            $suffix = q{};
        }
        $sheet = $prefix . $sheet . $suffix;

        if (exists $maps{$sheet}) {
            my ($E, $N)  = @{$maps{$sheet}->{bbox}[0]};  # NB we need the bbox corner so that it is left and below all points on the map

            if (defined $numbers) {
                my ($e, $n) = _get_eastnorthings($numbers);
                $E = $E + ($e-$E) % MINOR_GRID_SQ_SIZE;
                $N = $N + ($n-$N) % MINOR_GRID_SQ_SIZE;

                my $w = _winding_number($E, $N, $maps{$sheet}->{polygon});
                if ($w == 0) {
                    croak sprintf 'Grid reference %s = (%d, %d) is not on sheet %s', scalar format_grid($E,$N), $e, $n, $sheet;
                }
            }
            return wantarray ? ($E, $N) : "$E $N";
        }
    }

    # just a pair of numbers
    if ( @out = $s =~ m{\A (\d+(?:\.\d+)?) \s+ (\d+(?:\.\d+)?) \Z}xsmio ) { # eee nnn
        return wantarray ? @out : "@out";
    }

    croak "Failed to parse a grid reference from $s";
}

*parse_trad_grid       = \&parse_grid;
*parse_GPS_grid        = \&parse_grid;
*parse_landranger_grid = \&parse_grid;
*parse_map_grid        = \&parse_grid;

# is $pt left of $a--$b?
sub _is_left {
    my ($x, $y, $a, $b) = @_;
    return ( ($b->[0] - $a->[0]) * ($y - $a->[1]) - ($x - $a->[0]) * ($b->[1] - $a->[1]) );
}

# adapted from http://geomalgorithms.com/a03-_inclusion.html
sub _winding_number {
    my ($x, $y, $poly) = @_;
    my $w = 0;
    for (my $i=0; $i < $#$poly; $i++ ) {
        if ( $poly->[$i][1] <= $y ) {
            if ($poly->[$i+1][1]  > $y && _is_left($x, $y, $poly->[$i], $poly->[$i+1]) > 0 ) {
                $w++;
            }
        }
        else {
            if ($poly->[$i+1][1] <= $y && _is_left($x, $y, $poly->[$i], $poly->[$i+1]) < 0 ) {
                $w--;
            }
        }
    }
    return $w;
}

1;

__END__