igc_lib.py

"""A simple library for parsing IGC files.

Main abstraction defined in this file is the Flight class, which
represents a parsed IGC file. A Flight is a collection of:
    - GNNSFix objects, one per B record in the original file,
    - IGC metadata, extracted from A/I/H records
    - a list of detected Thermals,
    - a list of detected Glides.

A Flight is assumed to have one actual flight,
from takeoff to landing. If it has zero then it will
be invalid, i.e. `Flight.valid` will be False.

If there's more than one actual flight in the IGC file then
the `which_flight_to_pick` option in FlightParsingConfig
will determine behavior.

For example usage see the attached igc_lib_demo.py file. Please note
that after creating a Flight instance you should always check for its
validity via the `Flight.valid` attribute prior to using it, as many
IGC records are broken. See `Flight.notes` for details on why a file
was considered broken.
"""

from __future__ import print_function

import pandas as pd
import collections
import datetime
import math
import re
import xml.dom.minidom
from pathlib2 import Path

from collections import defaultdict

import lib.viterbi as viterbi
import lib.geo as geo


def _strip_non_printable_chars(string):
    """Filters a string removing non-printable characters.

    Args:
        string: A string to be filtered.

    Returns:
        A string, where non-printable characters are removed.
    """
    printable = set("0123456789abcdefghijklmnopqrstuvwxyzABCDEFGHIJKL"
                    "MNOPQRSTUVWXYZ!\"#$%&\'()*+,-./:;<=>?@[\\]^_`{|}~ ")

    printable_string = [x for x in string if x in printable]
    return ''.join(printable_string)


def _rawtime_float_to_hms(timef):
    """Converts time from floating point seconds to hours/minutes/seconds.

    Args:
        timef: A floating point time in seconds to be converted

    Returns:
        A namedtuple with hours, minutes and seconds elements
    """
    time = int(round(timef))
    hms = collections.namedtuple('hms', ['hours', 'minutes', 'seconds'])

    return hms((time/3600), (time % 3600)/60, time % 60)


class Turnpoint:
    """A single turnpoint in a Task.

    Attributes:
        lat: a float, latitude in degrees
        lon: a float, longitude in degrees
        radius: a float, radius of cylinder or line in km
        kind: type of turnpoint; "start_exit", "start_enter", "cylinder",
        "End_of_speed_section", "goal_cylinder", "goal_line"
    """

    def __init__(self, lat, lon, radius, kind):
        self.lat = lat
        self.lon = lon
        self.radius = radius
        self.kind = kind
        assert kind in ["start_exit", "start_enter", "cylinder",
                        "End_of_speed_section", "goal_cylinder",
                        "goal_line"], \
            "turnpoint type is not valid: %r" % kind

    def in_radius(self, fix):
        """Checks whether the provided GNSSFix is within the radius"""
        distance = geo.earth_distance(self.lat, self.lon, fix.lat, fix.lon)
        return distance < self.radius


class Task:
    """Stores a single flight task definition

    Checks if a Flight has achieved the turnpoints in the Task.

    Attributes:
        turnpoints: A list of Turnpoint objects.
        start_time: Raw time (seconds past midnight). The time the race starts.
                    The pilots must start at or after this time.
        end_time: Raw time (seconds past midnight). The time the race ends.
                  The pilots must finish the race at or before this time.
                  No credit is given for distance covered after this time.
    """

    @staticmethod
    def create_from_lkt_file(filename):
        """ Creates Task from LK8000 task file, which is in xml format.
            LK8000 does not have End of Speed Section or task finish time.
            For the goal, at the moment, Turnpoints can't handle goal cones or
            lines, for this reason we default to goal_cylinder.
        """

        # Open XML document using minidom parser
        DOMTree = xml.dom.minidom.parse(filename)
        task = DOMTree.documentElement

        # Get the taskpoints, waypoints and time gate
        # TODO: add code to handle if these tags are missing.
        taskpoints = task.getElementsByTagName("taskpoints")[0]
        waypoints = task.getElementsByTagName("waypoints")[0]
        gate = task.getElementsByTagName("time-gate")[0]
        tpoints = taskpoints.getElementsByTagName("point")
        wpoints = waypoints.getElementsByTagName("point")
        start_time = gate.getAttribute("open-time")

        start_hours, start_minutes = start_time.split(':')
        start_time = int(start_hours) * 3600 + int(start_minutes) * 60
        end_time = 23*3600 + 59*60 + 59  # default end_time of 23:59:59

        # Create a dictionary of names and a list of longitudes and latitudes
        # as the waypoints co-ordinates are stored separate to turnpoint
        # details.
        coords = defaultdict(list)

        for point in wpoints:
            name = point.getAttribute("name")
            longitude = float(point.getAttribute("longitude"))
            latitude = float(point.getAttribute("latitude"))
            coords[name].append(longitude)
            coords[name].append(latitude)

        # Create list of turnpoints
        turnpoints = []
        for point in tpoints:
            lat = coords[point.getAttribute("name")][1]
            lon = coords[point.getAttribute("name")][0]
            radius = float(point.getAttribute("radius"))/1000

            if point == tpoints[0]:
                # It is the first turnpoint, the start
                if point.getAttribute("Exit") == "true":
                    kind = "start_exit"
                else:
                    kind = "start_enter"
            else:
                if point == tpoints[-1]:
                    # It is the last turnpoint, i.e. the goal
                    if point.getAttribute("type") == "line":
                        # TODO(kuaka): change to 'line' once we can process it
                        kind = "goal_cylinder"
                    else:
                        kind = "goal_cylinder"
                else:
                    # All turnpoints other than the 1st and the last are
                    # "cylinders". In theory they could be
                    # "End_of_speed_section" but this is not supported by
                    # LK8000. For paragliders it would be safe to assume
                    # that the 2nd to last is always "End_of_speed_section".
                    kind = "cylinder"

            turnpoint = Turnpoint(lat, lon, radius, kind)
            turnpoints.append(turnpoint)
        task = Task(turnpoints, start_time, end_time)
        return task

    def __init__(self, turnpoints, start_time, end_time):
        self.turnpoints = turnpoints
        self.start_time = start_time
        self.end_time = end_time

    def check_flight(self, flight):
        """ Checks a Flight object against the task.

            Args:
                flight: a Flight object

            Returns:
                a list of GNSSFixes of when turnpoints were achieved.
        """
        reached_turnpoints = []
        proceed_to_start = False
        t = 0
        for fix in flight.fixes:
            if t >= len(self.turnpoints):
                # Pilot has arrived in goal (last turnpoint) so we can stop.
                break

            if self.end_time < fix.rawtime:
                # Task has ended
                break

            # Pilot must have at least 1 fix inside the start after the start
            # time, then exit.
            if self.turnpoints[t].kind == "start_exit":
                if proceed_to_start:
                    if not self.turnpoints[t].in_radius(fix):
                        reached_turnpoints.append(fix)  # pilot has started
                        t += 1
                if fix.rawtime > self.start_time and not proceed_to_start:
                    if self.turnpoints[t].in_radius(fix):
                        # Pilot is inside start after the start time.
                        proceed_to_start = True

            # Pilot must have at least 1 fix outside the start after
            # the start time, then enter.
            elif self.turnpoints[t].kind == "start_enter":
                if proceed_to_start:
                    if self.turnpoints[t].in_radius(fix):
                        # Pilot has started
                        reached_turnpoints.append(fix)
                        t += 1
                if fix.rawtime > self.start_time and not proceed_to_start:
                    if not self.turnpoints[t].in_radius(fix):
                        # Pilot is outside start after the start time.
                        proceed_to_start = True

            elif self.turnpoints[t].kind in ["cylinder",
                                             "End_of_speed_section",
                                             "goal_cylinder"]:
                if self.turnpoints[t].in_radius(fix):
                    # pilot has achieved turnpoint
                    reached_turnpoints.append(fix)
                    t += 1
            else:
                assert False, (
                    "Unknown turnpoint kind: %s" % self.turnpoints[t].kind)

        return reached_turnpoints


class GNSSFix:
    """Stores single GNSS flight recorder fix (a B-record).

    Raw attributes (i.e. attributes read directly from the B record):
        rawtime: a float, time since last midnight, UTC, seconds
        lat: a float, latitude in degrees
        lon: a float, longitude in degrees
        validity: a string, GPS validity information from flight recorder
        press_alt: a float, pressure altitude, meters
        gnss_alt: a float, GNSS altitude, meters
        extras: a string, B record extensions

    Derived attributes:
        index: an integer, the position of the fix in the IGC file
        timestamp: a float, true timestamp (since epoch), UTC, seconds
        alt: a float, either press_alt or gnss_alt
        gsp: a float, current ground speed, km/h
        bearing: a float, aircraft bearing, in degrees
        bearing_change_rate: a float, bearing change rate, degrees/second
        flying: a bool, whether this fix is during a flight
        circling: a bool, whether this fix is inside a thermal
    """

    @staticmethod
    def build_from_B_record(B_record_line, index):
        """Creates GNSSFix object from IGC B-record line.

        Args:
            B_record_line: a string, B record line from an IGC file
            index: the zero-based position of the fix in the parent IGC file

        Returns:
            The created GNSSFix object
        """
        match = re.match(
            '^B' + '(\d\d)(\d\d)(\d\d)'
            + '(\d\d)(\d\d)(\d\d\d)([NS])'
            + '(\d\d\d)(\d\d)(\d\d\d)([EW])'
            + '([AV])' + '([-\d]\d\d\d\d)' + '([-\d]\d\d\d\d)'
            + '([0-9a-zA-Z\-]*).*$', B_record_line)
        if match is None:
            return None
        (hours, minutes, seconds,
         lat_deg, lat_min, lat_min_dec, lat_sign,
         lon_deg, lon_min, lon_min_dec, lon_sign,
         validity, press_alt, gnss_alt,
         extras) = match.groups()

        rawtime = (float(hours)*60.0 + float(minutes))*60.0 + float(seconds)

        lat = float(lat_deg)
        lat += float(lat_min) / 60.0
        lat += float(lat_min_dec) / 1000.0 / 60.0
        if lat_sign == 'S':
            lat = -lat

        lon = float(lon_deg)
        lon += float(lon_min) / 60.0
        lon += float(lon_min_dec) / 1000.0 / 60.0
        if lon_sign == 'W':
            lon = -lon

        press_alt = float(press_alt)
        gnss_alt = float(gnss_alt)

        return GNSSFix(rawtime, lat, lon, validity, press_alt, gnss_alt,
                       index, extras)

    def __init__(self, rawtime, lat, lon, validity, press_alt, gnss_alt,
                 index, extras):
        """Initializer of GNSSFix. Not meant to be used directly."""
        self.rawtime = rawtime
        self.lat = lat
        self.lon = lon
        self.validity = validity
        self.press_alt = press_alt
        self.gnss_alt = gnss_alt
        self.index = index
        self.extras = extras
        self.flight = None

    def set_flight(self, flight):
        """Sets parent Flight object."""
        self.flight = flight
        if self.flight.alt_source == "PRESS":
            self.alt = self.press_alt
        elif self.flight.alt_source == "GNSS":
            self.alt = self.gnss_alt
        else:
            assert(False)
        self.timestamp = self.rawtime + flight.date_timestamp

    def __repr__(self):
        return self.__str__()

    def __str__(self):
        return (
            "GNSSFix(rawtime=%02d:%02d:%02d, lat=%f, lon=%f, press_alt=%.1f, gnss_alt=%.1f)" %
            (_rawtime_float_to_hms(self.rawtime) +
             (self.lat, self.lon, self.press_alt, self.gnss_alt)))

    def bearing_to(self, other):
        """Computes bearing in degrees to another GNSSFix."""
        return geo.bearing_to(self.lat, self.lon, other.lat, other.lon)

    def distance_to(self, other):
        """Computes great circle distance in kilometers to another GNSSFix."""
        return geo.earth_distance(self.lat, self.lon, other.lat, other.lon)

    def to_B_record(self):
        """Reconstructs an IGC B-record."""
        rawtime = int(self.rawtime)
        hours = rawtime / 3600
        minutes = (rawtime % 3600) / 60
        seconds = rawtime % 60

        if self.lat < 0.0:
            lat = -self.lat
            lat_sign = 'S'
        else:
            lat = self.lat
            lat_sign = 'N'
        lat = int(round(lat*60000.0))
        lat_deg = lat / 60000
        lat_min = (lat % 60000) / 1000
        lat_min_dec = lat % 1000

        if self.lon < 0.0:
            lon = -self.lon
            lon_sign = 'W'
        else:
            lon = self.lon
            lon_sign = 'E'
        lon = int(round(lon*60000.0))
        lon_deg = lon / 60000
        lon_min = (lon % 60000) / 1000
        lon_min_dec = lon % 1000

        validity = self.validity
        press_alt = int(self.press_alt)
        gnss_alt = int(self.gnss_alt)
        extras = self.extras

        return (
            "B" +
            "%02d%02d%02d" % (hours, minutes, seconds) +
            "%02d%02d%03d%s" % (lat_deg, lat_min, lat_min_dec, lat_sign) +
            "%03d%02d%03d%s" % (lon_deg, lon_min, lon_min_dec, lon_sign) +
            validity +
            "%05d%05d" % (press_alt, gnss_alt) +
            extras)


class Thermal:
    """Represents a single thermal detected in a flight.

    Attributes:
        enter_fix: a GNSSFix, entry point of the thermal
        exit_fix: a GNSSFix, exit point of the thermal
    """
    def __init__(self, enter_fix, exit_fix):
        self.enter_fix = enter_fix
        self.exit_fix = exit_fix

    def time_change(self):
        """Returns the time spent in the thermal, seconds."""
        return self.exit_fix.rawtime - self.enter_fix.rawtime

    def alt_change(self):
        """Returns the altitude gained/lost in the thermal, meters."""
        return self.exit_fix.alt - self.enter_fix.alt

    def vertical_velocity(self):
        """Returns average vertical velocity in the thermal, m/s."""
        if math.fabs(self.time_change()) < 1e-7:
            return 0.0
        return self.alt_change() / self.time_change()

    def __repr__(self):
        return self.__str__()

    def __str__(self):
        hms = _rawtime_float_to_hms(self.time_change())
        return ("Thermal(vertical_velocity=%.2f m/s, duration=%dm %ds)" %
                (self.vertical_velocity(), hms.minutes, hms.seconds))


class Glide:
    """Represents a single glide detected in a flight.

    Glides are portions of the recorded track between thermals.

    Attributes:
        enter_fix: a GNSSFix, entry point of the glide
        exit_fix: a GNSSFix, exit point of the glide
        track_length: a float, the total length, in kilometers, of the recorded
        track, between the entry point and the exit point; note that this is
        not the same as the distance between these points
    """

    def __init__(self, enter_fix, exit_fix, track_length):
        self.enter_fix = enter_fix
        self.exit_fix = exit_fix
        self.track_length = track_length

    def time_change(self):
        """Returns the time spent in the glide, seconds."""
        return self.exit_fix.timestamp - self.enter_fix.timestamp
    
### homebrew start ###

# We need something, that returns singel elements to avoid str cutting

### homebrew end ###

    def speed(self):
        """Returns the average speed in the glide, km/h."""
        return self.track_length / (self.time_change() / 3600.0)

    def alt_change(self):
        """Return the overall altitude change in the glide, meters."""
        return self.enter_fix.alt - self.exit_fix.alt

    def glide_ratio(self):
        """Returns the L/D of the glide."""
        if math.fabs(self.alt_change()) < 1e-7:
            return 0.0
        return (self.track_length * 1000.0) / self.alt_change()
    
    def __repr__(self):
        return self.__str__()

    def __str__(self):
        hms = _rawtime_float_to_hms(self.time_change())
        
        return (
            ("Glide(dist=%.2f km, avg_speed=%.2f kph, "
             "avg L/D=%.2f duration=%dm %ds)") % (
                self.track_length, self.speed(), self.glide_ratio(),
                hms.minutes, hms.seconds))
                


class FlightParsingConfig(object):
    """Configuration for parsing an IGC file.

    Defines a set of parameters used to validate a file, and to detect
    thermals and flight mode. Details in comments.
    """

    #
    # Flight validation parameters.
    #

    # Minimum number of fixes in a file.
    min_fixes = 50

    # Maximum time between fixes, seconds.
    # Soft limit, some fixes are allowed to exceed.
    max_seconds_between_fixes = 50.0

    # Minimum time between fixes, seconds.
    # Soft limit, some fixes are allowed to exceed.
    min_seconds_between_fixes = 1.0

    # Maximum number of fixes exceeding time between fix constraints.
    max_time_violations = 10

    # Maximum number of times a file can cross the 0:00 UTC time.
    max_new_days_in_flight = 2

    # Minimum average of absolute values of altitude changes in a file.
    # This is needed to discover altitude sensors (either pressure or
    # gps) that report either always constant altitude, or almost
    # always constant altitude, and therefore are invalid. The unit
    # is meters/fix.
    min_avg_abs_alt_change = 0.01

    # Maximum altitude change per second between fixes, meters per second.
    # Soft limit, some fixes are allowed to exceed.
    max_alt_change_rate = 50.0

    # Maximum number of fixes that exceed the altitude change limit.
    max_alt_change_violations = 3

    # Absolute maximum altitude, meters.
    max_alt = 10000.0

    # Absolute minimum altitude, meters.
    min_alt = -600.0

    #
    # Flight detection parameters.
    #

    # Minimum ground speed to switch to flight mode, km/h.
    min_gsp_flight = 15.0

    # Minimum idle time (i.e. time with speed below min_gsp_flight) to switch
    # to landing, seconds. Exception: end of the file (tail fixes that
    # do not trigger the above condition), no limit is applied there.
    min_landing_time = 5.0 * 60.0

    # In case there are multiple continuous segments with ground
    # speed exceeding the limit, which one should be taken?
    # Available options:
    #   - "first": take the first segment, ignore the part after
    #     the first detected landing.
    #   - "concat": concatenate all segments; will include the down
    #     periods between segments (legacy behavior)
    which_flight_to_pick = "concat"

    #
    # Thermal detection parameters.
    #

    # Minimum bearing change to enter a thermal, deg/sec.
    min_bearing_change_circling = 6.0

    # Minimum time between fixes to calculate bearing change, seconds.
    # See the usage for a more detailed comment on why this is useful.
    min_time_for_bearing_change = 5.0

    # Minimum time to consider circling a thermal, seconds.
    min_time_for_thermal = 60.0


class Flight:
    """Parses IGC file, detects thermals and checks for record anomalies.

    Before using an instance of Flight check the `valid` attribute. An
    invalid Flight instance is not usable. For an explaination why is
    a Flight invalid see the `notes` attribute.

    General attributes:
        valid: a bool, whether the supplied record is considered valid
        notes: a list of strings, warnings and errors encountered while
        parsing/validating the file
        fixes: a list of GNSSFix objects, one per each valid B record
        thermals: a list of Thermal objects, the detected thermals
        glides: a list of Glide objects, the glides between thermals
        takeoff_fix: a GNSSFix object, the fix at which takeoff was detected
        landing_fix: a GNSSFix object, the fix at which landing was detected

    IGC metadata attributes (some might be missing if the flight does not
    define them):
        glider_type: a string, the declared glider type
        competition_class: a string, the declared competition class
        fr_manuf_code: a string, the flight recorder manufaturer code
        fr_uniq_id: a string, the flight recorded unique id
        i_record: a string, the I record (describing B record extensions)
        fr_firmware_version: a string, the version of the recorder firmware
        fr_hardware_version: a string, the version of the recorder hardware
        fr_recorder_type: a string, the type of the recorder
        fr_gps_receiver: a string, the used GPS receiver
        fr_pressure_sensor: a string, the used pressure sensor

    Other attributes:
        alt_source: a string, the chosen altitude sensor,
        either "PRESS" or "GNSS"
        press_alt_valid: a bool, whether the pressure altitude sensor is OK
        gnss_alt_valid: a bool, whether the GNSS altitude sensor is OK
    """

    @staticmethod
    def create_from_file(filename, config_class=FlightParsingConfig):
        """Creates an instance of Flight from a given file.

        Args:
            filename: a string, the name of the input IGC file
            config_class: a class that implements FlightParsingConfig

        Returns:
            An instance of Flight built from the supplied IGC file.
        """
        config = config_class()
        fixes = []
        a_records = []
        i_records = []
        h_records = []
        abs_filename = Path(filename).expanduser().absolute()
        with abs_filename.open('r', encoding="ISO-8859-1") as flight_file:
            for line in flight_file:
                line = line.replace('\n', '').replace('\r', '')
                if not line:
                    continue
                if line[0] == 'A':
                    a_records.append(line)
                elif line[0] == 'B':
                    fix = GNSSFix.build_from_B_record(line, index=len(fixes))
                    if fix is not None:
                        if fixes and math.fabs(fix.rawtime - fixes[-1].rawtime) < 1e-5:
                            # The time did not change since the previous fix.
                            # Ignore this fix.
                            pass
                        else:
                            fixes.append(fix)
                elif line[0] == 'I':
                    i_records.append(line)
                elif line[0] == 'H':
                    h_records.append(line)
                else:
                    # Do not parse any other types of IGC records
                    pass
        flight = Flight(fixes, a_records, h_records, i_records, config)
        return flight

    def __init__(self, fixes, a_records, h_records, i_records, config):
        """Initializer of the Flight class. Do not use directly."""
        self._config = config
        self.fixes = fixes
        self.valid = True
        self.notes = []
        if len(fixes) < self._config.min_fixes:
            self.notes.append(
                "Error: This file has %d fixes, less than "
                "the minimum %d." % (len(fixes), self._config.min_fixes))
            self.valid = False
            return

        self._check_altitudes()
        if not self.valid:
            return

        self._check_fix_rawtime()
        if not self.valid:
            return

        if self.press_alt_valid:
            self.alt_source = "PRESS"
        elif self.gnss_alt_valid:
            self.alt_source = "GNSS"
        else:
            self.notes.append(
                "Error: neither pressure nor gnss altitude is valid.")
            self.valid = False
            return

        if a_records:
            self._parse_a_records(a_records)
        if i_records:
            self._parse_i_records(i_records)
        if h_records:
            self._parse_h_records(h_records)

        if not hasattr(self, 'date_timestamp'):
            self.notes.append("Error: no date record (HFDTE) in the file")
            self.valid = False
            return

        for fix in self.fixes:
            fix.set_flight(self)

        self._compute_ground_speeds()
        self._compute_flight()
        self._compute_takeoff_landing()
        if not hasattr(self, 'takeoff_fix'):
            self.notes.append("Error: did not detect takeoff.")
            self.valid = False
            return

        self._compute_bearings()
        self._compute_bearing_change_rates()
        self._compute_circling()
        self._find_thermals()

    def _parse_a_records(self, a_records):
        """Parses the IGC A record.

        A record contains the flight recorder manufacturer ID and
        device unique ID.
        """
        self.fr_manuf_code = _strip_non_printable_chars(a_records[0][1:4])
        self.fr_uniq_id = _strip_non_printable_chars(a_records[0][4:7])

    def _parse_i_records(self, i_records):
        """Parses the IGC I records.

        I records contain a description of extensions used in B records.
        """
        self.i_record = _strip_non_printable_chars(" ".join(i_records))

    def _parse_h_records(self, h_records):
        """Parses the IGC H records.

        H records (header records) contain a lot of interesting metadata
        about the file, such as the date of the flight, name of the pilot,
        glider type, competition class, recorder accuracy and more.
        Consult the IGC manual for details.
        """
        for record in h_records:
            self._parse_h_record(record)

    def _parse_h_record(self, record):
        if record[0:5] == 'HFDTE':
            match = re.match(
                '(?:HFDTE|HFDTEDATE:[ ]*)(\d\d)(\d\d)(\d\d)',
                record, flags=re.IGNORECASE)
            if match:
                dd, mm, yy = [_strip_non_printable_chars(group) for group in match.groups()]
                year = int(2000 + int(yy))
                month = int(mm)
                day = int(dd)
                if 1 <= month <= 12 and 1 <= day <= 31:
                    epoch = datetime.datetime(year=1970, month=1, day=1)
                    date = datetime.datetime(year=year, month=month, day=day)
                    self.date_timestamp = (date - epoch).total_seconds()
        elif record[0:5] == 'HFGTY':
            match = re.match(
                'HFGTY[ ]*GLIDER[ ]*TYPE[ ]*:[ ]*(.*)',
                record, flags=re.IGNORECASE)
            if match:
                (self.glider_type,) = map(
                    _strip_non_printable_chars, match.groups())
        elif record[0:5] == 'HFRFW' or record[0:5] == 'HFRHW':
            match = re.match(
                'HFR[FH]W[ ]*FIRMWARE[ ]*VERSION[ ]*:[ ]*(.*)',
                record, flags=re.IGNORECASE)
            if match:
                (self.fr_firmware_version,) = map(
                    _strip_non_printable_chars, match.groups())
            match = re.match(
                'HFR[FH]W[ ]*HARDWARE[ ]*VERSION[ ]*:[ ]*(.*)',
                record, flags=re.IGNORECASE)
            if match:
                (self.fr_hardware_version,) = map(
                    _strip_non_printable_chars, match.groups())
        elif record[0:5] == 'HFFTY':
            match = re.match(
                'HFFTY[ ]*FR[ ]*TYPE[ ]*:[ ]*(.*)',
                record, flags=re.IGNORECASE)
            if match:
                (self.fr_recorder_type,) = map(_strip_non_printable_chars,
                                               match.groups())
        elif record[0:5] == 'HFGPS':
            match = re.match(
                'HFGPS(?:[: ]|(?:GPS))*(.*)',
                record, flags=re.IGNORECASE)
            if match:
                (self.fr_gps_receiver,) = map(_strip_non_printable_chars,
                                              match.groups())
        elif record[0:5] == 'HFPRS':
            match = re.match(
                'HFPRS[ ]*PRESS[ ]*ALT[ ]*SENSOR[ ]*:[ ]*(.*)',
                record, flags=re.IGNORECASE)
            if match:
                (self.fr_pressure_sensor,) = map(_strip_non_printable_chars,
                                                 match.groups())
        elif record[0:5] == 'HFCCL':
            match = re.match(
                'HFCCL[ ]*COMPETITION[ ]*CLASS[ ]*:[ ]*(.*)',
                record, flags=re.IGNORECASE)
            if match:
                (self.competition_class,) = map(_strip_non_printable_chars,
                                                match.groups())

    def __str__(self):
        descr = "Flight(valid=%s, fixes: %d" % (
            str(self.valid), len(self.fixes))
        if hasattr(self, 'thermals'):
            descr += ", thermals: %d" % len(self.thermals)
        descr += ")"
        return descr

    def _check_altitudes(self):
        press_alt_violations_num = 0
        gnss_alt_violations_num = 0
        press_huge_changes_num = 0
        gnss_huge_changes_num = 0
        press_chgs_sum = 0.0
        gnss_chgs_sum = 0.0
        for i in range(len(self.fixes) - 1):
            press_alt_delta = math.fabs(
                self.fixes[i+1].press_alt - self.fixes[i].press_alt)
            gnss_alt_delta = math.fabs(
                self.fixes[i+1].gnss_alt - self.fixes[i].gnss_alt)
            rawtime_delta = math.fabs(
                self.fixes[i+1].rawtime - self.fixes[i].rawtime)
            if rawtime_delta > 0.5:
                if (press_alt_delta / rawtime_delta >
                        self._config.max_alt_change_rate):
                    press_huge_changes_num += 1
                else:
                    press_chgs_sum += press_alt_delta
                if (gnss_alt_delta / rawtime_delta >
                        self._config.max_alt_change_rate):
                    gnss_huge_changes_num += 1
                else:
                    gnss_chgs_sum += gnss_alt_delta
            if (self.fixes[i].press_alt > self._config.max_alt
                    or self.fixes[i].press_alt < self._config.min_alt):
                press_alt_violations_num += 1
            if (self.fixes[i].gnss_alt > self._config.max_alt or
                    self.fixes[i].gnss_alt < self._config.min_alt):
                gnss_alt_violations_num += 1
        press_chgs_avg = press_chgs_sum / float(len(self.fixes) - 1)
        gnss_chgs_avg = gnss_chgs_sum / float(len(self.fixes) - 1)

        press_alt_ok = True
        if press_chgs_avg < self._config.min_avg_abs_alt_change:
            self.notes.append(
                "Warning: average pressure altitude change between fixes "
                "is: %f. It is lower than the minimum: %f."
                % (press_chgs_avg, self._config.min_avg_abs_alt_change))
            press_alt_ok = False

        if press_huge_changes_num > self._config.max_alt_change_violations:
            self.notes.append(
                "Warning: too many high changes in pressure altitude: %d. "
                "Maximum allowed: %d."
                % (press_huge_changes_num,
                   self._config.max_alt_change_violations))
            press_alt_ok = False

        if press_alt_violations_num > 0:
            self.notes.append(
                "Warning: pressure altitude limits exceeded in %d fixes."
                % (press_alt_violations_num))
            press_alt_ok = False

        gnss_alt_ok = True
        if gnss_chgs_avg < self._config.min_avg_abs_alt_change:
            self.notes.append(
                "Warning: average gnss altitude change between fixes is: %f. "
                "It is lower than the minimum: %f."
                % (gnss_chgs_avg, self._config.min_avg_abs_alt_change))
            gnss_alt_ok = False

        if gnss_huge_changes_num > self._config.max_alt_change_violations:
            self.notes.append(
                "Warning: too many high changes in gnss altitude: %d. "
                "Maximum allowed: %d."
                % (gnss_huge_changes_num,
                   self._config.max_alt_change_violations))
            gnss_alt_ok = False

        if gnss_alt_violations_num > 0:
            self.notes.append(
                "Warning: gnss altitude limits exceeded in %d fixes." %
                gnss_alt_violations_num)
            gnss_alt_ok = False

        self.press_alt_valid = press_alt_ok
        self.gnss_alt_valid = gnss_alt_ok

    def _check_fix_rawtime(self):
        """Checks for rawtime anomalies, fixes 0:00 UTC crossing.

        The B records do not have fully qualified timestamps (just the current
        time in UTC), therefore flights that cross 0:00 UTC need special
        handling.
        """
        DAY = 24.0 * 60.0 * 60.0
        days_added = 0
        rawtime_to_add = 0.0
        rawtime_between_fix_exceeded = 0
        for i in range(1, len(self.fixes)):
            f0 = self.fixes[i-1]
            f1 = self.fixes[i]
            f1.rawtime += rawtime_to_add

            if (f0.rawtime > f1.rawtime and
                    f1.rawtime + DAY < f0.rawtime + 200.0):
                # Day switch
                days_added += 1
                rawtime_to_add += DAY
                f1.rawtime += DAY

            time_change = f1.rawtime - f0.rawtime
            if time_change < self._config.min_seconds_between_fixes - 1e-5:
                rawtime_between_fix_exceeded += 1
            if time_change > self._config.max_seconds_between_fixes + 1e-5:
                rawtime_between_fix_exceeded += 1

        if rawtime_between_fix_exceeded > self._config.max_time_violations:
            self.notes.append(
                "Error: too many fixes intervals exceed time between fixes "
                "constraints. Allowed %d fixes, found %d fixes."
                % (self._config.max_time_violations,
                   rawtime_between_fix_exceeded))
            self.valid = False
        if days_added > self._config.max_new_days_in_flight:
            self.notes.append(
                "Error: too many times did the flight cross the UTC 0:00 "
                "barrier. Allowed %d times, found %d times."
                % (self._config.max_new_days_in_flight, days_added))
            self.valid = False

    def _compute_ground_speeds(self):
        """Adds ground speed info (km/h) to self.fixes."""
        self.fixes[0].gsp = 0.0
        for i in range(1, len(self.fixes)):
            dist = self.fixes[i].distance_to(self.fixes[i-1])
            rawtime = self.fixes[i].rawtime - self.fixes[i-1].rawtime
            if math.fabs(rawtime) < 1e-5:
                self.fixes[i].gsp = 0.0
            else:
                self.fixes[i].gsp = dist/rawtime*3600.0

    def _flying_emissions(self):
        """Generates raw flying/not flying emissions from ground speed.

        Standing (i.e. not flying) is encoded as 0, flying is encoded as 1.
        Exported to a separate function to be used in Baum-Welch parameters
        learning.
        """
        emissions = []
        for fix in self.fixes:
            if fix.gsp > self._config.min_gsp_flight:
                emissions.append(1)
            else:
                emissions.append(0)

        return emissions

    def _compute_flight(self):
        """Adds boolean flag .flying to self.fixes.

        Two pass:
          1. Viterbi decoder
          2. Only emit landings (0) if the downtime is more than
             _config.min_landing_time (or it's the end of the log).
        """
        # Step 1: the Viterbi decoder
        emissions = self._flying_emissions()
        decoder = viterbi.SimpleViterbiDecoder(
            # More likely to start the log standing, i.e. not in flight
            init_probs=[0.80, 0.20],
            transition_probs=[
                [0.9995, 0.0005],  # transitions from standing
                [0.0005, 0.9995],  # transitions from flying
            ],
            emission_probs=[
                [0.8, 0.2],  # emissions from standing
                [0.2, 0.8],  # emissions from flying
            ])

        outputs = decoder.decode(emissions)

        # Step 2: apply _config.min_landing_time.
        ignore_next_downtime = False
        apply_next_downtime = False
        for i, (fix, output) in enumerate(zip(self.fixes, outputs)):
            if output == 1:
                fix.flying = True
                # We're in flying mode, therefore reset all expectations
                # about what's happening in the next down mode.
                ignore_next_downtime = False
                apply_next_downtime = False
            else:
                if apply_next_downtime or ignore_next_downtime:
                    if apply_next_downtime:
                        fix.flying = False
                    else:
                        fix.flying = True
                else:
                    # We need to determine whether to apply_next_downtime
                    # or to ignore_next_downtime. This requires a scan into
                    # upcoming fixes. Find the next fix on which
                    # the Viterbi decoder said "flying".
                    j = i + 1
                    while j < len(self.fixes):
                        upcoming_fix_decoded = outputs[j]
                        if upcoming_fix_decoded == 1:
                            break
                        j += 1

                    if j == len(self.fixes):
                        # No such fix, end of log. Then apply.
                        apply_next_downtime = True
                        fix.flying = False
                    else:
                        # Found next flying fix.
                        upcoming_fix = self.fixes[j]
                        upcoming_fix_time_ahead = upcoming_fix.rawtime - fix.rawtime
                        # If it's far enough into the future of then apply.
                        if upcoming_fix_time_ahead >= self._config.min_landing_time:
                            apply_next_downtime = True
                            fix.flying = False
                        else:
                            ignore_next_downtime = True
                            fix.flying = True

    def _compute_takeoff_landing(self):
        """Finds the takeoff and landing fixes in the log.

        Takeoff fix is the first fix in the flying mode. Landing fix
        is the next fix after the last fix in the flying mode or the
        last fix in the file.
        """
        takeoff_fix = None
        landing_fix = None
        was_flying = False
        for fix in self.fixes:
            if fix.flying and takeoff_fix is None:
                takeoff_fix = fix
            if not fix.flying and was_flying:
                landing_fix = fix
                if self._config.which_flight_to_pick == "first":
                    # User requested to select just the first flight in the log,
                    # terminate now.
                    break
            was_flying = fix.flying

        if takeoff_fix is None:
            # No takeoff found.
            return

        if landing_fix is None:
            # Landing on the last fix
            landing_fix = self.fixes[-1]

        self.takeoff_fix = takeoff_fix
        self.landing_fix = landing_fix

    def _compute_bearings(self):
        """Adds bearing info to self.fixes."""
        for i in range(len(self.fixes) - 1):
            self.fixes[i].bearing = self.fixes[i].bearing_to(self.fixes[i+1])
        self.fixes[-1].bearing = self.fixes[-2].bearing

    def _compute_bearing_change_rates(self):
        """Adds bearing change rate info to self.fixes.

        Computing bearing change rate between neighboring fixes proved
        itself to be noisy on tracks recorded with minimum interval (1 second).
        Therefore we compute rates between points that are at least
        min_time_for_bearing_change seconds apart.
        """
        def find_prev_fix(curr_fix):
            """Computes the previous fix to be used in bearing rate change."""
            prev_fix = None
            for i in range(curr_fix - 1, 0, -1):
                time_dist = math.fabs(self.fixes[curr_fix].timestamp -
                                      self.fixes[i].timestamp)
                if (time_dist >
                        self._config.min_time_for_bearing_change - 1e-7):
                    prev_fix = i
                    break
            return prev_fix

        for curr_fix in range(len(self.fixes)):
            prev_fix = find_prev_fix(curr_fix)

            if prev_fix is None:
                self.fixes[curr_fix].bearing_change_rate = 0.0
            else:
                bearing_change = (self.fixes[prev_fix].bearing -
                                  self.fixes[curr_fix].bearing)
                if math.fabs(bearing_change) > 180.0:
                    if bearing_change < 0.0:
                        bearing_change += 360.0
                    else:
                        bearing_change -= 360.0
                time_change = (self.fixes[prev_fix].timestamp -
                               self.fixes[curr_fix].timestamp)
                change_rate = bearing_change/time_change
                self.fixes[curr_fix].bearing_change_rate = change_rate

    def _circling_emissions(self):
        """Generates raw circling/straight emissions from bearing change.

        Staight flight is encoded as 0, circling is encoded as 1. Exported
        to a separate function to be used in Baum-Welch parameters learning.
        """
        emissions = []
        for fix in self.fixes:
            bearing_change = math.fabs(fix.bearing_change_rate)
            bearing_change_enough = (
                bearing_change > self._config.min_bearing_change_circling)
            if fix.flying and bearing_change_enough:
                emissions.append(1)
            else:
                emissions.append(0)
        return emissions

    def _compute_circling(self):
        """Adds .circling to self.fixes."""
        emissions = self._circling_emissions()
        decoder = viterbi.SimpleViterbiDecoder(
            # More likely to start in straight flight than in circling
            init_probs=[0.80, 0.20],
            transition_probs=[
                [0.982, 0.018],  # transitions from straight flight
                [0.030, 0.970],  # transitions from circling
            ],
            emission_probs=[
                [0.942, 0.058],  # emissions from straight flight
                [0.093, 0.907],  # emissions from circling
            ])

        output = decoder.decode(emissions)

        for i in range(len(self.fixes)):
            self.fixes[i].circling = (output[i] == 1)

    def _find_thermals(self):
        """Go through the fixes and find the thermals.

        Every point not in a thermal is put into a glide.If we get to end of
        the fixes and there is still an open glide (i.e. flight not finishing
        in a valid thermal) the glide will be closed.
        """
        takeoff_index = self.takeoff_fix.index
        landing_index = self.landing_fix.index
        flight_fixes = self.fixes[takeoff_index:landing_index + 1]

        self.thermals = []
        self.glides = []
        circling_now = False
        gliding_now = False
        first_fix = None
        first_glide_fix = None
        last_glide_fix = None
        distance = 0.0
        for fix in flight_fixes:
            if not circling_now and fix.circling:
                # Just started circling
                circling_now = True
                first_fix = fix
                distance_start_circling = distance
            elif circling_now and not fix.circling:
                # Just ended circling
                circling_now = False
                thermal = Thermal(first_fix, fix)
                if (thermal.time_change() >
                        self._config.min_time_for_thermal - 1e-5):
                    self.thermals.append(thermal)
                    # glide ends at start of thermal
                    glide = Glide(first_glide_fix, first_fix,
                                  distance_start_circling)
                    self.glides.append(glide)
                    gliding_now = False

            if gliding_now:
                distance = distance + fix.distance_to(last_glide_fix)
                last_glide_fix = fix
            else:
                # just started gliding
                first_glide_fix = fix
                last_glide_fix = fix
                gliding_now = True
                distance = 0.0

        if gliding_now:
            glide = Glide(first_glide_fix, last_glide_fix, distance)
            self.glides.append(glide)