File Coverage

blib/lib/Geo/Coordinates/KKJ.pm

Criterion	Covered	Total	%
statement	107	107	100.0
branch	29	32	90.6
condition	3	6	50.0
subroutine	13	13	100.0
pod	8	8	100.0
total	160	166	96.3

package Geo::Coordinates::KKJ;

use 5.008006;
use strict;
use warnings;
use Carp;
use Math::Trig;

require Exporter;

our $VERSION = '0.01';

our @ISA = qw(Exporter);

our @EXPORT = qw(
    KKJxy_to_WGS84lalo
    WGS84lalo_to_KKJxy
    KKJxy_to_KKJlalo
    KKJlalo_to_KKJxy
    KKJlalo_to_WGS84lalo
    WGS84lalo_to_KKJlalo
    KKJ_Zone_I
    KKJ_Zone_Lo
);

my $KKJ_ZONE_INFO = {
    '0' => [ '18.0', '500_000.0' ],
    '1' => [ '21.0', '1_500_000.0' ],
    '2' => [ '24.0', '2_500_000.0' ],
    '3' => [ '27.0', '3_500_000.0' ],
    '4' => [ '30.0', '4_500_000.0' ],
    '5' => [ '33.0', '5_500_000.0' ],
};

sub KKJxy_to_WGS84lalo {
    croak 'Geo::Coordinates::KKJ::KKJxy_to_WGS84lalo needs two arguments'
        if @_ != 2;

    my ( $x,    $y )    = @_;
    my ( $lat,  $lon )  = KKJxy_to_KKJlalo( $x, $y );
    my ( $wlat, $wlon ) = KKJlalo_to_WGS84lalo( $lat, $lon );

    return ( $wlat, $wlon );
}

sub WGS84lalo_to_KKJxy {
    croak 'Geo::Coordinates::KKJ::WGS84lalo_to_KKJxy needs two arguments'
        if @_ != 2;

    my ( $lat, $lon ) = @_;
    my ( $kkjlat, $kkjlon ) = WGS84lalo_to_KKJlalo( $lat, $lon );

    my $zone = KKJ_Zone_Lo($kkjlon);

    my $LALO = {
        'La' => $kkjlat,
        'Lo' => $kkjlon,
    };

    my $foo = KKJlalo_to_KKJxy( $LALO, $zone );

    return ( $foo->{'P'}, $foo->{'I'} );

}

sub KKJxy_to_KKJlalo {
    croak 'Geo::Coordinates::KKJ::KKJxy_to_KKJlalo needs two arguments'
        if @_ != 2;

    my ( $p, $i ) = @_;

    croak "Wrong coordenates" if ( $p < $i );

    my $zone = KKJ_Zone_I($i);
    my $LALO = {};

    my $min_la = deg2rad('59.0');
    my $max_la = deg2rad('70.5');
    my $min_lo = deg2rad('18.5');
    my $max_lo = deg2rad('32.0');

    #Scan iteratively the target area, until find matching
    #KKJ coordinate value.  Area is defined with Hayford Ellipsoid.
    for ( my $count = 1 ; $count < 35 ; $count++ ) {
        my $delta_la = $max_la - $min_la;
        my $delta_lo = $max_lo - $min_lo;

        $LALO->{'La'} = rad2deg( $min_la + 0.5 * $delta_la );
        $LALO->{'Lo'} = rad2deg( $min_lo + 0.5 * $delta_lo );

        my $KKJt = KKJlalo_to_KKJxy( $LALO, $zone );

        if ( $KKJt->{'P'} < $p ) {
            $min_la = $min_la + 0.45 * $delta_la;
        }
        else {
            $max_la = $min_la + 0.55 * $delta_la;
        }

        if ( $KKJt->{'I'} < $i ) {
            $min_lo = $min_lo + 0.45 * $delta_lo;
        }
        else {
            $max_lo = $min_lo + 0.55 * $delta_lo;
        }
    }

    return ( $LALO->{'La'}, $LALO->{'Lo'} );
}

sub KKJlalo_to_KKJxy {
    croak 'Geo::Coordinates::KKJ::KKJlalo_to_KKJxy needs two arguments'
        if @_ != 2;

    my ( $INP, $zone ) = @_;

    my $Lo = deg2rad( $INP->{'Lo'} ) - deg2rad( $KKJ_ZONE_INFO->{$zone}->[0] );

    # Hayford ellipsoid
    my $a = '6378388.0';
    my $f = 1 / 297.0;

    my $b  = ( 1.0 - $f ) * $a;
    my $bb = $b * $b;
    my $c  = ( $a / $b ) * $a;
    my $ee = ( $a * $a - $bb ) / $bb;
    my $n  = ( $a - $b ) / ( $a + $b );
    my $nn = $n * $n;

    my $cosLa = cos( deg2rad( $INP->{'La'} ) );

    my $NN = $ee * $cosLa * $cosLa;

    my $LaF =
      Math::Trig::atan( Math::Trig::tan( deg2rad( $INP->{'La'} ) ) /
          cos( $Lo * sqrt( 1 + $NN ) ) );

    my $cosLaF = cos($LaF);

    my $t =
      ( Math::Trig::tan($Lo) * $cosLaF ) / sqrt( 1 + $ee * $cosLaF * $cosLaF );

    my $A  = $a / ( 1 + $n );
    my $A1 = $A * ( 1 + $nn / 4 + $nn * $nn / 64 );
    my $A2 = $A * 1.5 * $n *     ( 1 - $nn / 8 );
    my $A3 = $A * 0.9375 * $nn * ( 1 - $nn / 4 );
    my $A4 = $A * 35 / 48.0 * $nn * $n;

    my $OUT = {};

    $OUT->{'P'} =
      $A1 * $LaF -
      $A2 * sin( 2 * $LaF ) +
      $A3 * sin( 4 * $LaF ) -
      $A4 * sin( 6 * $LaF );

    $OUT->{'I'} =
      $c * log( $t + sqrt( 1 + $t * $t ) ) + 500_000.0 + $zone * 1_000_000.0;

    return $OUT;
}

sub KKJlalo_to_WGS84lalo {
    croak 'Geo::Coordinates::KKJ::KKJlalo_to_WGS84lalo needs two arguments'
        if @_ != 2;

    my ( $La, $Lo ) = @_;

    my $d_la =
      deg2rad( 0.124867E+01 +
          -0.269982E+00 * $La +
          0.191330E+00 * $Lo +
          0.356119E-02 * $La * $La +
          -0.122312E-02 * $La * $Lo +
          -0.335514E-03 * $Lo * $Lo ) / 3600.0;

    my $d_lo =
      deg2rad( -0.286111E+02 +
          0.114183E+01 * $La +
          -0.581428E+00 * $Lo +
          -0.152421E-01 * $La * $La +
          0.118177E-01 * $La * $Lo +
          0.826646E-03 * $Lo * $Lo ) / 3600.0;

    my $WGS = {};
    $WGS->{'La'} = rad2deg( deg2rad($La) + $d_la );
    $WGS->{'Lo'} = rad2deg( deg2rad($Lo) + $d_lo );

    return ( $WGS->{'La'}, $WGS->{'Lo'} );
}

sub WGS84lalo_to_KKJlalo {
    croak 'Geo::Coordinates::KKJ::WGS84lalo_to_KKJlalo needs two arguments'
        if @_ != 2;

    my ( $La, $Lo ) = @_;
    croak "Wrong parameters"
        if ( $La < $Lo );

    my $d_la =
      deg2rad( -0.124766E+01 +
          0.269941E+00 * $La +
          -0.191342E+00 * $Lo +
          -0.356086E-02 * $La * $La +
          0.122353E-02 * $La * $Lo +
          0.335456E-03 * $Lo * $Lo ) / 3600.0;

    my $d_lo =
      deg2rad( 0.286008E+02 +
          -0.114139E+01 * $La +
          0.581329E+00 * $Lo +
          0.152376E-01 * $La * $La +
          -0.118166E-01 * $La * $Lo +
          -0.826201E-03 * $Lo * $Lo ) / 3600.0;

    my $KKJ = {};
    $KKJ->{'La'} = rad2deg( deg2rad($La) + $d_la );
    $KKJ->{'Lo'} = rad2deg( deg2rad($Lo) + $d_lo );

    return ( $KKJ->{'La'}, $KKJ->{'Lo'} );
}

sub KKJ_Zone_I {
    croak 'Geo::Coordinates::KKJ::KKJ_Zone_I needs one argument'
        if @_ != 1;

    my $i    = shift;
    my $zone = int( $i / 1_000_000.0 );
    croak "Zone value ($zone) invalid."
        if ( $zone < 0 || $zone > 5 );
    return $zone;
}

sub KKJ_Zone_Lo {
    croak 'Geo::Coordinates::KKJ::KKJ_Zone_Lo needs one argument'
        if @_ != 1;

    my $Lo   = shift;

    croak 'Longitude is undefined' unless defined($Lo);

    my $zone = 5;
    while ( $zone >= 0 ) {
        if ( abs( $Lo - $KKJ_ZONE_INFO->{$zone}->[0] ) <= 1.5 ) {
            last;
        }
        $zone--;
    }

    if ( $zone >= 0 && $zone <= 5 ) {
        return $zone;
    }
    else {
        croak "Zone outside range";
    }
}

1;

__END__