Pygeodesy Testing

triangulation.py

# pip install pygeodesy
import pygeodesy.sphericalNvector


def intersection_point(lat_long_1, bearing_1, lat_long_2, bearing_2):
    vector1 = pygeodesy.sphericalNvector.LatLon(*lat_long_1)
    vector2 = pygeodesy.sphericalNvector.LatLon(*lat_long_2)
    return vector1.intersection(bearing_1, vector2, bearing_2)


# tests that make sense for intersection_point()
def test_north_pole():
    # any two gps locations with bearing 0 should intersect at north pole, latitude 90.
    gps_1 = (52.520, 13.405)  # berlin
    gps_2 = (49.246, -123.116)  # vancouver
    gps_3 = (31.224, 121.469)  # shanghai
    bearing = 0
    int1 = intersection_point(gps_1, bearing, gps_2, bearing)
    int2 = intersection_point(gps_1, bearing, gps_3, bearing)
    int3 = intersection_point(gps_2, bearing, gps_3, bearing)
    print(int1, int2, int3)
    assert int1.lat == 90.0
    assert int2.lat == 90.0
    assert int3.lat == 90.0
    print(
        f"longitude: {int1.lon, int2.lon, int3.lon}"
    )  # long of north pole is 0 for this library


def test_small_field():
    # tests some points on thunderbird field. Most similar to our use case.
    # https://imgur.com/a/KcbBFlJ
    # Fix B, choose A1 and A2 and record bearing for each and see if A1, A2, bearings give back B or close
    B = (49.25475, -123.24557)

    A1 = (49.25433, -123.24585)
    A1_bearing = 23.52

    A2 = (49.25437, -123.24548)
    A2_bearing = 351.21

    inter = intersection_point(A1, A1_bearing, A2, A2_bearing)
    print(f"Actual: {B}")
    print(f"Predicted: {inter}")
    print(f"Lat difference: {inter.lat - B[0]}")
    print(f"Long difference: {inter.lon - B[1]}")


def test_country():
    B = (43.86622, -99.9699)

    A1 = (41.5085800, -115.5716600)
    A1_bearing = 78.38

    A2 = (39.57182, -94.92699)
    A2_bearing = 318.78

    inter = intersection_point(A1, A1_bearing, A2, A2_bearing)
    print(f"Actual: {B}")
    print(f"Predicted: {inter}")
    # demonstrates that this method is not very accurate for long distances, which is fine because
    # we are only flying over a small field.


if __name__ == "__main__":
    # test_north_pole()
    # test_small_field()
    test_country()