File Coverage

blib/lib/Geo/Coordinates/RDNAP.pm

Criterion	Covered	Total	%
statement	86	86	100.0
branch	6	12	50.0
condition	8	22	36.3
subroutine	16	16	100.0
pod	4	4	100.0
total	120	140	85.7

package Geo::Coordinates::RDNAP;

use strict;
use warnings;

use vars qw($VERSION);
$VERSION = '0.11';

use Carp;
use Params::Validate qw/validate BOOLEAN SCALAR/;

use Exporter;
use vars qw/@ISA @EXPORT_OK/;
@ISA = qw/Exporter/;
@EXPORT_OK = qw/from_rd to_rd deg dms/;

sub deg {
    my (@in) = @_;
    my @out;

    while (my($d, $m, $s) = splice (@in, 0, 3)) {
        push @out, $d + ($m||0)/60 + ($s||0)/3600;
    }
    return @out;
}

sub dms {
    return map {int($_), int($_*60)%60, ($_-int($_*60)/60)*3600} @_;
}

my %a = (
    '01' => 3236.0331637,
    20 => -32.5915821,
    '02' => -0.2472814,
    21 => -0.8501341,
    '03' => -0.0655238,
    22 => -0.0171137,
    40 =>  0.0052771,
    23 => -0.0003859,
    41 =>  0.0003314,
    '04' =>  0.0000371,
    42 =>  0.0000143,
    24 => -0.0000090,
);

my %b = (
    10 => 5261.3028966,
    11 => 105.9780241,
    12 =>  2.4576469,
    30 => -0.8192156,
    31 => -0.0560092,
    13 =>  0.0560089,
    32 => -0.0025614,
    14 =>  0.0012770,
    50 =>  0.0002574,
    33 => -0.0000973,
    51 =>  0.0000293,
    15 =>  0.0000291,
);

my %c = (
    '01' => 190066.98903,
    11 => -11830.85831,
    21 => -114.19754,
    '03' => -32.38360,
    31 => -2.34078,
    13 => -0.60639,
    23 => 0.15774,
    41 => -0.04158,
    '05' => -0.00661,
);

my %d = (
    10 => 309020.31810,
    '02' => 3638.36193,
    12 => -157.95222,
    20 => 72.97141,
    30 => 59.79734,
    22 => -6.43481,
    '04' => 0.09351,
    32 => -0.07379,
    14 => -0.05419,
    40 => -0.03444,
);

my %bessel = (
    a   => 6377397.155,
    e2  => 6674372e-9,
    f_i => 299.1528128,
);

my %etrs89 = (
    a   => 6378137,
    e2  => 6694380e-9,
    f_i => 298.257222101,
);

# Transformation parameters from Bessel to ETRS89 with respect to
# Amersfoort.

my %b2e = (
    tx  => 593.032,
    ty  => 26,
    tz  => 478.741,
    a   => 1.9848e-6,
    b   => -1.7439e-6,
    c   => 9.0587e-6,
    d   => 4.0772e-6,
);

my %e2b = map {$_ => -$b2e{$_}} keys %b2e;

my @amersfoort_b = ( 3903_453.148, 368_135.313, 5012_970.306 );
my @amersfoort_e = ( 3904_046.180, 368_161.313, 5013_449.047 );

sub from_rd {
    croak 'Geo::Coordinates::RDNAP::from_rd needs two or three arguments'
        if (@_ !=2 && @_ != 3);

    my ($x, $y, $h) = (@_, 0);

    croak "Geo::Coordinates::RDNAP::from_rd: X out of bounds: $x"
        if ($x < -7_000 or $x > 300_000);
    croak "Geo::Coordinates::RDNAP::from_rd: Y out of bounds: $y"
        if ($y < 289_000 or $y > 629_000);

    # Use the approximated transformation.
    # Step 1: RD -> Bessel (spherical coords)

    $x = ($x/100_000) - 1.55;
    $y = ($y/100_000) - 4.63;

    my $lat = (52*60*60) + (9*60) + 22.178;
    my $lon = (5 *60*60) + (23*60) + 15.5;

    foreach my $i (keys %a) {
        my ($m, $n) = split //, $i;
        $lat += $a{$i} * ($x**$m) * ($y**$n);
    }

    foreach my $i (keys %b) {
        my ($m, $n) = split //, $i;
        $lon += $b{$i} * ($x**$m) * ($y**$n);
    }

    # Step 2: spherical coords -> X, Y, Z
    my @coords = _ellipsoid_to_cartesian($lat/3600, $lon/3600, $h, \%bessel);

    # Step 3: Bessel -> ETRS89
    @coords = _transform_datum( @coords, \%b2e, \@amersfoort_b );

    # Step 4: X, Y, Z -> spherical coords
    return _cartesian_to_ellipsoid(@coords, \%etrs89);
}

sub to_rd {
    croak 'Geo::Coordinates::RDNAP::to_rd needs two or three arguments'
        if (@_ !=2 && @_ != 3);

    my ($lat, $lon, $h) = (@_, 0);

    # Use the approximated transformation.
    # Step 1: spherical coords -> X, Y, Z
    my @coords = _ellipsoid_to_cartesian($lat, $lon, $h, \%etrs89);

    # Step 2: ETRS89 -> Bessel
    @coords = _transform_datum( @coords, \%e2b, \@amersfoort_e );

    # Step 3: X, Y, Z -> spherical coords
    ($lat, $lon, $h) = _cartesian_to_ellipsoid(@coords, \%bessel);

    # Step 4: Bessel -> RD'

    # Convert to units of 10_000 secs; as deltas from Amersfoort.
    $lat = ($lat * 3600 - ((52*60*60) + (9*60) + 22.178))/10_000;
    $lon = ($lon * 3600 - ((5 *60*60) + (23*60) + 15.5))/10_000;

    my $x = 155e3;
    my $y = 463e3;

    foreach my $i (keys %c) {
        my ($m, $n) = split //, $i;
        $x += $c{$i} * ($lat**$m) * ($lon**$n);
    }

    foreach my $i (keys %d) {
        my ($m, $n) = split //, $i;
        $y += $d{$i} * ($lat**$m) * ($lon**$n);
    }

    croak "Geo::Coordinates::RDNAP::to_rd: X out of bounds: $x"
        if ($x < -7_000 or $x > 300_000);
    croak "Geo::Coordinates::RDNAP::to_rd: Y out of bounds: $y"
        if ($y < 289_000 or $y > 629_000);

    return ($x, $y, $h);
}

sub _to_rads {
    return $_[0] * 2*3.14159_26535_89793 /360;
}

sub _from_rads {
    return $_[0] / (2*3.14159_26535_89793) *360;
}

sub _ellipsoid_to_cartesian {
    my ($lat, $lon, $h, $para) = @_;

    my $sinphi = sin(_to_rads($lat));
    my $cosphi = cos(_to_rads($lat));
    my $n = $para->{a}/sqrt(1 - $para->{e2}*$sinphi*$sinphi);

    return (($n+$h)*$cosphi*cos(_to_rads($lon)),
            ($n+$h)*$cosphi*sin(_to_rads($lon)),
            ($n*(1-$para->{e2})+$h)*$sinphi );
}

# Returns (lat, lon, h) in degrees.

sub _cartesian_to_ellipsoid {
    my ($x, $y, $z, $para) = @_;

    my $lon = atan2($y, $x);

    my $r = sqrt($x*$x+$y*$y);
    my $phi = 0;
    my $n_sinphi = $z;
    my $n;
    my $oldphi;

    do {
        $oldphi = $phi;
        $phi = atan2($z + $para->{e2}*$n_sinphi, $r);
        my $sinphi = sin($phi);
        $n = $para->{a}/sqrt(1-$para->{e2}*$sinphi*$sinphi);
        $n_sinphi = $n*$sinphi;
    } while (abs($oldphi-$phi) > 1e-8);

    my $h = $r/cos($phi) - $n;

    return (_from_rads($phi), _from_rads($lon), $h);
}

sub _transform_datum {
    my ($x, $y, $z, $t, $centre) = @_;

    return (
        $x + $t->{d}*($x-$centre->[0]) + $t->{c}*($y-$centre->[1])
            - $t->{b}*($z-$centre->[2]) + $t->{tx},
        $y - $t->{c}*($x-$centre->[0]) + $t->{d}*($y-$centre->[1])
            + $t->{a}*($z-$centre->[2]) + $t->{ty},
        $z + $t->{b}*($x-$centre->[0]) - $t->{a}*($y-$centre->[1])
            + $t->{d}*($z-$centre->[2]) + $t->{tz}
    );
}

1;
__END__