File Coverage

blib/lib/Geo/Lookup/ByTime.pm

Criterion	Covered	Total	%
statement	124	129	96.1
branch	44	54	81.4
condition	20	36	55.5
subroutine	19	20	95.0
pod	6	6	100.0
total	213	245	86.9

package Geo::Lookup::ByTime;

use warnings;
use strict;
use Carp;
use Scalar::Util qw(blessed);
use base qw(Exporter);

our @EXPORT_OK = qw(hav_distance);

our $VERSION = '0.10';

use constant EARTH_RADIUS => 6_378_137.0;
use constant PI           => 4 * atan2( 1, 1 );
use constant DEG_TO_RAD   => PI / 180.0;
use constant RAD_TO_DEG   => 180.0 / PI;

sub new {
    my $class = shift;
    my $self  = {
        points    => [],
        need_sort => 0
    };

    bless( $self, $class );

    if ( @_ ) {
        $self->add_points( @_ );
    }

    return $self;
}

sub add_points {
    my $self = shift;

    $self->{need_sort}++ if @_;

    for my $pt ( @_ ) {
        if (   blessed( $pt )
            && $pt->can( 'latitude' )
            && $pt->can( 'longitude' )
            && $pt->can( 'time' ) ) {
            push @{ $self->{points} },
              {
                lat  => $pt->latitude(),
                lon  => $pt->longitude(),
                time => $pt->time(),
                orig => $pt
              };
        }
        elsif ( ref( $pt ) eq 'CODE' ) {
            my @pts = ();
            while ( my $ipt = $pt->() ) {
                push @pts, $ipt;
                if ( @pts >= 100 ) {
                    # Add points 100 at a time.
                    $self->add_points( @pts );
                    @pts = ();
                }
            }
            $self->add_points( @pts );
        }
        elsif ( ref( $pt ) eq 'ARRAY' ) {
            $self->add_points( @{$pt} );
        }
        elsif ( ref( $pt ) eq 'HASH' ) {
            croak(
                "Point hashes must have the following keys: lat, lon, time\n"
              )
              unless exists( $pt->{lat} )
                  && exists( $pt->{lon} )
                  && exists( $pt->{time} );
            push @{ $self->{points} }, $pt;
        }
        else {
            croak( "Don't know how to add "
                  . ( defined( $pt ) ? $pt : '(undef)' ) );
        }
    }

    return;
}

sub get_points {
    my $self = shift;

    if ( $self->{need_sort} ) {
        my $np = [
            sort { $a->{time} <=> $b->{time} }
              grep {
                     defined( $_->{lat} )
                  && defined( $_->{lon} )
                  && defined( $_->{time} )
              } @{ $self->{points} }
        ];
        $self->{points}    = $np;
        $self->{need_sort} = 0;
    }

    return $self->{points};
}

# Returns the index of the first point with time >= the supplied time
sub _search {
    my $pts  = shift;
    my $time = shift;

    my $max = scalar( @{$pts} );
    my ( $lo, $mid, $hi ) = ( 0, 0, $max - 1 );

    TRY:
    while ( $lo <= $hi ) {
        $mid = int( ( $lo + $hi ) / 2 );
        my $cmp = $pts->[$mid]->{time} <=> $time;
        if ( $cmp < 0 ) {
            $lo = $mid + 1;
        }
        elsif ( $cmp > 0 ) {
            $hi = $mid - 1;
        }
        else {
            last TRY;
        }
    }

    while ( $mid < $max && $pts->[$mid]->{time} < $time ) {
        $mid++;
    }

    return ( $mid < $max ) ? $mid : undef;
}

sub _interp {
    my ( $lo, $mid, $hi, $val1, $val2 ) = @_;
    confess "$lo <= $mid <= $hi !"
      unless $lo <= $mid && $mid <= $hi;
    my $scale = $hi - $lo;
    my $posn  = $mid - $lo;
    return ( $val1 * ( $scale - $posn ) + $val2 * $posn ) / $scale;
}

sub nearest {
    my $self     = shift;
    my $time     = shift;
    my $max_dist = shift;

    my $pts = $self->get_points();
    my $pos = _search( $pts, $time );

    return unless defined( $pos );

    if ( $pts->[$pos]->{time} == $time ) {
        # Exact match - just return the point
        my $pt = {
            lat  => $pts->[$pos]->{lat},
            lon  => $pts->[$pos]->{lon},
            time => $time
        };

        return
          wantarray
          ? ( $pt, $pts->[$pos]->{orig} || $pts->[$pos], 0 )
          : $pt;
    }

    # If we're at the first point we can't
    # interpolate with anything.
    return if $pos == 0;

    my ( $p1, $p2 ) = @$pts[ $pos - 1, $pos ];

    # Linear interpolation between nearest points
    my $lat = _interp( $p1->{time}, $time, $p2->{time}, $p1->{lat},
        $p2->{lat} );
    my $lon = _interp( $p1->{time}, $time, $p2->{time}, $p1->{lon},
        $p2->{lon} );

    my $pt = {
        lat  => $lat,
        lon  => $lon,
        time => $time
    };

    my $best_dist = 0;
    my $best      = undef;

    # Compute nearest if we need to return it or check proximity
    if ( wantarray || defined( $max_dist ) ) {
        my $d1 = abs( $pt->{time} - $p1->{time} );
        my $d2 = abs( $pt->{time} - $p2->{time} );

        $best = ( $d1 < $d2 ) ? $p1 : $p2;
        $best_dist = hav_distance( $pt, $best );

        # Nearest point out of range?
        return if defined( $max_dist ) && $best_dist > $max_dist;
    }

    # Return a synthetic point
    return
      wantarray ? ( $pt, $best->{orig} || $best, $best_dist ) : $pt;
}

sub _deg {
    return map { $_ * RAD_TO_DEG } @_;
}

sub _rad {
    return map { $_ * DEG_TO_RAD } @_;
}

# From
#  http://perldoc.perl.org/functions/sin.html
sub _asin {
    return atan2( $_[0], sqrt( 1 - $_[0] * $_[0] ) );
}

# Not a method
sub hav_distance {
    my $dist = 0;
    my ( $lat1, $lon1 );
    while ( my $pt = shift ) {
        my ( $lat2, $lon2 ) = _rad( $pt->{lat}, $pt->{lon} );
        if ( defined( $lat1 ) ) {
            my $sdlat = sin( ( $lat1 - $lat2 ) / 2.0 );
            my $sdlon = sin( ( $lon1 - $lon2 ) / 2.0 );
            my $res = sqrt( $sdlat * $sdlat
                  + cos( $lat1 ) * cos( $lat2 ) * $sdlon * $sdlon );
            if ( $res > 1.0 ) {
                $res = 1.0;
            }
            elsif ( $res < -1.0 ) {
                $res = -1.0;
            }
            $dist += 2.0 * _asin( $res );
        }
        ( $lat1, $lon1 ) = ( $lat2, $lon2 );
    }

    return $dist * EARTH_RADIUS;
}

sub time_range {
    my $self = shift;
    my $pts  = $self->get_points();
    return unless @{$pts};
    return ( $pts->[0]->{time}, $pts->[-1]->{time} );
}

1;
__END__