Map my gps data

sailmap.py

#! /usr/bin/env python
# vim:fenc=utf-8
#
# Copyright © 2024 fred <[email protected]>
#
# Distributed under terms of the BSD 3-Clause license.

"""
Send my GPS traces on a map.
Examples: https://bsdworld.org/misc/
"""
import argparse
import math
import pathlib
import webbrowser
from zoneinfo import ZoneInfo

import folium
import folium.plugins
import gpxpy
import gpxpy.gpx
import numpy as np


class SlidingWindow:
  def __init__(self, iterable, window_size=2):
    self.iterable = iterable
    self.window_size = window_size
    self.start = 0
    self.end = len(iterable)

  def __iter__(self):
    return self

  def __next__(self):
    if self.start + self.window_size > self.end:
      raise StopIteration
    window = self.iterable[self.start:self.start + self.window_size]
    self.start += 1
    return window


def haversine(lat1, lon1, lat2, lon2):
  R = 6371.0
  lat1, lon1, lat2, lon2 = map(math.radians, [lat1, lon1, lat2, lon2])
  dlat = lat2 - lat1
  dlon = lon2 - lon1

  a = math.sin(dlat / 2)**2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon / 2)**2
  c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
  distance = R * c
  return distance


def find_farthest_point(ref_point, gps_trace):
  farthest_point = None
  max_distance = 0

  ref_lat, ref_lon = ref_point.latitude, ref_point.longitude

  for pts in gps_trace:
    distance = haversine(ref_lat, ref_lon, pts.latitude, pts.longitude)
    if distance > max_distance:
      max_distance = distance
      farthest_point = pts

  return farthest_point, max_distance


def every_nth_element(lst, n):
  return [lst[i] for i in range(0, len(lst), n)]


def read_gpx(filename):
  localtz = ZoneInfo('America/Los_Angeles')
  with filename.open('r') as gpx_fd:
    gpx_data = gpxpy.parse(gpx_fd)
    points = []
    for track in gpx_data.tracks:
      for segment in track.segments:
        for point in segment.points:
          point.time = point.time.astimezone(localtz)
          points.append(point)
    return points


def add_markers(filename, chart):
  try:
    with filename.open('r') as gpx_fd:
      gpx_data = gpxpy.parse(gpx_fd)
  except FileNotFoundError as err:
    print(err)
    return

  icon_html = """
  <div style="font-size: 14px; color: darkgreen;">
  <i class="fa fa-flag"></i>
  </div>"""

  for marker in gpx_data.waypoints:
    icon_only = folium.DivIcon(html=icon_html)
    folium.Marker(location=(marker.latitude, marker.longitude),
                  name=marker.name, popup=marker.comment,
                  icon=icon_only).add_to(chart)


def calculate_heading(lat1, lon1, lat2, lon2):
  delta_lon = lon2 - lon1
  x = math.sin(math.radians(delta_lon)) * math.cos(math.radians(lat2))
  y = (math.cos(math.radians(lat1)) * math.sin(math.radians(lat2)) -
       math.sin(math.radians(lat1)) * math.cos(math.radians(lat2)) *
       math.cos(math.radians(delta_lon)))
  heading = math.degrees(math.atan2(x, y))
  return (heading + 360) % 360


def detect_tacks(points, threshold_angle=90, sample_factor=10):
  last_heading = None
  cumulative_change = 0
  tack_in_progress = False

  tacks = [points[0]]
  for i in range(0, len(points) - sample_factor, sample_factor):
    pt1 = points[i]
    pt2 = points[i + sample_factor]
    heading = calculate_heading(pt1.latitude, pt1.longitude, pt2.latitude, pt2.longitude)

    if last_heading is not None:
      # Calculate the change in heading
      change = heading - last_heading
      if abs(change) > 180:
        change -= np.sign(change) * 360  # Normalize to -180 to 180
      cumulative_change += change

      # Check if the cumulative change indicates a tack
      if abs(cumulative_change) >= threshold_angle:
        if not tack_in_progress:
          tacks.append(pt1)
          tack_in_progress = True
          cumulative_change = 0
        else:
          tack_in_progress = False
    last_heading = heading

  tacks.append(points[-1])
  return tacks


def nbsp(var):
  return var.replace(' ', '&nbsp;') + ' '


def main():
  parser = argparse.ArgumentParser()
  parser.add_argument('--start', type=int, default=0,
                      help="number of points to cut at the beginning.")
  parser.add_argument('--end', type=int, default=0,
                      help="number of points to cut at the end.")
  parser.add_argument('-w', '--way-points', type=pathlib.Path, default=None,
                      help='Waypoints and markers')
  parser.add_argument('filename', type=pathlib.Path, nargs=1)
  parser.add_argument('-o', '--output', type=pathlib.Path, default=None,
                      help='Output filename')
  opts = parser.parse_args()

  try:
    gpx_file = opts.filename[0].absolute()
    out_file = opts.output.absolute() or gpx_file.with_suffix('.html')
    points = read_gpx(gpx_file)
  except FileNotFoundError as err:
    raise SystemExit(err) from None

  points = points[opts.start:-opts.end]
  tacks = detect_tacks(points, threshold_angle=80, sample_factor=8)

  farthest, distance = find_farthest_point(points[0], points)
  mid_point = points.index(farthest)

  chart = folium.Map()

  folium.TileLayer(tiles='USGS.USTopo', overlay=True, control=True).add_to(chart)
  folium.TileLayer(tiles='https://tiles.openseachart.org/seamark/{z}/{x}/{y}.png',
                   attr='OpenSeachart', name='OpenSeachart',
                   overlay=True, control=True).add_to(chart)

  trace = [(pts.latitude, pts.longitude) for pts in points]

  folium.PolyLine(trace[:mid_point]).add_to(chart)
  folium.PolyLine(trace[mid_point:], color='orange').add_to(chart)

  avg_speed = round(points[0].speed_between(farthest) * 1.994, 1)
  folium.Marker(location=(farthest.latitude, farthest.longitude),
                tooltip="Click me!",
                popup=(f"Distance:&nbsp;{round(distance, 1)}NM "
                       f'Average&nbsp;Speed:&nbsp;{avg_speed}'),
                icon=folium.Icon(color="red")).add_to(chart)

  for pt1, pt2 in SlidingWindow(tacks, 2):
    speed = pt1.speed_between(pt2)
    speed = round(speed * 1.994, 1)
    folium.Marker(location=(pt2.latitude, pt2.longitude),
                  tooltip=f"<b>Click Me</b> for more info!<br>Speed:&nbsp;{speed}&nbsp;Knots",
                  popup=nbsp((f"<b>Since the last marker</b>\n"
                              f"Speed: {speed} Knots\n"
                              f"Time: {pt2.time.strftime('%Y/%m/%d %H:%M:%S')}\n")),
                  icon=folium.Icon(color="blue", prefix='fa', icon="sailboat")).add_to(chart)

  if opts.way_points:
    add_markers(opts.way_points, chart)

  chart.fit_bounds(chart.get_bounds(), padding=(10, 10))
  chart.save(out_file)
  print(f'Writing: {out_file}')
  webbrowser.open(out_file.as_uri())


if __name__ == '__main__':
  main()