along/cross track error calculation: for validation of storm track forecast

track_errors.py

from dataclasses import dataclass
import numpy as np
from numpy import sin, cos


R = EARTH_RADIUS = 6371000

@dataclass
class Point:

    x: float
    y: float

    @property
    def xr(self):
        return np.deg2rad(self.x)

    @property
    def yr(self):
        return np.deg2rad(self.y)

    def angular_distance(self, p):
        haversine = sin((p.yr - self.yr) / 2) ** 2 + \
            cos(self.yr) * cos(p.yr) * sin((p.xr - self.xr) /2 ) ** 2
        c = 2 * np.arctan2(np.sqrt(haversine), np.sqrt(1 - haversine))
        return c

    def distance(self, p):
        return self.angular_distance(p) * R

    def bearing(self, p):
        return np.arctan2(sin(p.xr - self.xr) * cos(p.yr),
            cos(self.yr) * sin(p.yr) - \
                sin(self.yr) * cos(p.yr) * cos(p.xr - self.xr))


def get_error(actual_p0_tup, actual_p_tup, forecast_p_tup):
    actual_p0 = Point(*actual_p0_tup)
    actual_p = Point(*actual_p_tup)
    forecast_p = Point(*forecast_p_tup)
    if actual_p == actual_p0:
        # Stationary
        along_track_error = cross_track_error = actual_p.distance(forecast_p)
        return along_track_error, cross_track_error
    # cross track
    angular_distance = actual_p0.angular_distance(forecast_p)
    cross_track_error = abs(np.arcsin(sin(angular_distance) \
        * sin(actual_p0.bearing(forecast_p) - actual_p0.bearing(actual_p))) * R)
    along_track_error = np.arccos(cos(angular_distance) / cos(cross_track_error / R))\
        * R - actual_p0.distance(actual_p)
    return along_track_error, cross_track_error

if __name__ == '__main__':
    actual_p0 = (120, 20)
    actual_p  = (121, 20)
    forecast_p_1 = (120.8, 20.2)
    forecast_p_2 = (121.2, 19.8)
    print(get_error(actual_p0, actual_p, forecast_p_1))
    print(get_error(actual_p0, actual_p, forecast_p_2))