From a489bfa5f9c5b01d1b8c79d11e875d4265dd5f64 Mon Sep 17 00:00:00 2001 From: Raphael <raphael.sturgis@lis-lab.fr> Date: Thu, 10 Feb 2022 11:26:15 +0100 Subject: [PATCH] Ais trajectory shifting --- .gitignore | 3 + skais/ais/ais_trajectory.py | 53 ++++++++++++++++++ skais/tests/ais/test_ais_trajectory.py | 40 ++++++++++++- skais/tests/utils/test_geography.py | 39 +++++++++++++ .../utils/__pycache__/__init__.cpython-38.pyc | Bin 167 -> 0 bytes .../__pycache__/geography.cpython-38.pyc | Bin 746 -> 0 bytes skais/utils/geography.py | 20 ++++++- 7 files changed, 151 insertions(+), 4 deletions(-) create mode 100644 skais/tests/utils/test_geography.py delete mode 100644 skais/utils/__pycache__/__init__.cpython-38.pyc delete mode 100644 skais/utils/__pycache__/geography.cpython-38.pyc diff --git a/.gitignore b/.gitignore index b139f3e..400aa52 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,6 @@ build/ skais.egg-info/ *.coverage *__pycache__* +venv/ + +local.* \ No newline at end of file diff --git a/skais/ais/ais_trajectory.py b/skais/ais/ais_trajectory.py index de9b6b0..5b2af98 100644 --- a/skais/ais/ais_trajectory.py +++ b/skais/ais/ais_trajectory.py @@ -1,8 +1,11 @@ +import random + import pandas as pd import numpy as np from numba import jit from scipy.interpolate import interp1d +from skais.utils.geography import great_circle, position_from_distance from skais.ais.ais_points import AISPoints @@ -154,3 +157,53 @@ class AISTrajectory(AISPoints): index += i return result + + def shift_trajectory_to_coordinates(self, target_coordinate=(0, 0), point_index=None, in_place=False): + if point_index is None: + point_index = random.randint(0, len(self.df.index) - 1) + + df = self.df.copy() + new_df = df.copy() + + new_df['latitude'].iat[point_index] = target_coordinate[0] + new_df['longitude'].iat[point_index] = target_coordinate[1] + + new_point = target_coordinate + for i in range(point_index, 0, -1): + current_point = (df.iloc[i]['latitude'], df.iloc[i]['longitude']) + lat_dist = great_circle(current_point[0], df.iloc[i - 1]['latitude'], current_point[1], current_point[1]) + long_dist = great_circle(current_point[0], current_point[0], current_point[1], df.iloc[i - 1]['longitude']) + + if current_point[0] > df.iloc[i - 1]['latitude']: + lat_dist *= -1 + + if current_point[1] > df.iloc[i - 1]['longitude']: + long_dist *= -1 + + new_point = position_from_distance(new_point, (lat_dist, long_dist)) + + new_df['latitude'].iat[i - 1] = new_point[0] + new_df['longitude'].iat[i - 1] = new_point[1] + + new_point = target_coordinate + for i in range(point_index, len(df.index) - 1): + current_point = (df.iloc[i]['latitude'], df.iloc[i]['longitude']) + lat_dist = great_circle(current_point[0], df.iloc[i + 1]['latitude'], current_point[1], current_point[1]) + long_dist = great_circle(current_point[0], current_point[0], current_point[1], df.iloc[i + 1]['longitude']) + + if current_point[0] > df.iloc[i + 1]['latitude']: + lat_dist *= -1 + + if current_point[1] > df.iloc[i + 1]['longitude']: + long_dist *= -1 + + new_point = position_from_distance(new_point, (lat_dist, long_dist)) + + new_df['latitude'].iat[i + 1] = new_point[0] + new_df['longitude'].iat[i + 1] = new_point[1] + + if in_place: + self.df = new_df + return self + else: + return AISTrajectory(new_df, mmsi=self.mmsi) diff --git a/skais/tests/ais/test_ais_trajectory.py b/skais/tests/ais/test_ais_trajectory.py index 799825c..8d46240 100644 --- a/skais/tests/ais/test_ais_trajectory.py +++ b/skais/tests/ais/test_ais_trajectory.py @@ -319,4 +319,42 @@ class TestAISTrajectory(unittest.TestCase): self.assertEqual(0, compute_trajectory.py_func(times, 800)) def test_apply_func_on_window(self): - self.assertRaises(ValueError, apply_func_on_window,np.arange(10), 0, 0, 'not valid string') \ No newline at end of file + self.assertRaises(ValueError, apply_func_on_window,np.arange(10), 0, 0, 'not valid string') + + def test_shift_trajectory_to_coordinates_no_change(self): + trajectory = AISTrajectory( + pd.DataFrame( + { + "latitude": [0, 90, 0, -90], + "longitude": [0, 90, 180, -90], + "ts_sec": [i for i in range(4)] + } + ) + ) + + result = trajectory.shift_trajectory_to_coordinates(point_index=0) + + pd.testing.assert_frame_equal(trajectory.df.reset_index(drop=True), result.df.reset_index(drop=True)) + + def test_shift_trajectory_to_coordinates_change_1(self): + trajectory = AISTrajectory( + pd.DataFrame( + { + "latitude": [1, 2, 3, 4], + "longitude": [1, 2, 3, 4], + "ts_sec": [i for i in range(4)] + } + ) + ) + + result = trajectory.shift_trajectory_to_coordinates(point_index=0).df + + expected = pd.DataFrame( + { + "latitude": [0, 2, 3, 4], # TBD + "longitude": [0, 1, 2, 3], + "ts_sec": [i for i in range(4)] + } + ) + + pd.testing.assert_frame_equal(expected.reset_index(drop=True), result.reset_index(drop=True)) \ No newline at end of file diff --git a/skais/tests/utils/test_geography.py b/skais/tests/utils/test_geography.py new file mode 100644 index 0000000..8316524 --- /dev/null +++ b/skais/tests/utils/test_geography.py @@ -0,0 +1,39 @@ +import unittest +import numpy as np + +from skais.utils.geography import position_from_distance + +tol = 0.005 + + +class TestGeography(unittest.TestCase): + def test_position_from_distance_1(self): + position = (0, 0) + distance = (10000, 10000) + + expected = (0.09, 0.09) + result = position_from_distance(position, distance) + + np.testing.assert_allclose(expected, result, tol) + + def test_position_from_distance_2(self): + position = (10, 10) + distance = (10000, 10000) + + expected = (10.0636111, 10.064722222222223) + result = position_from_distance(position, distance) + + np.testing.assert_allclose(expected, result, tol) + + def test_position_from_distance_3(self): + position = (75, 75) + distance = (10000, -10000) + + expected = (75.0633333, 74.75333333333333) + result = position_from_distance(position, distance) + + np.testing.assert_allclose(expected, result, tol) + + +if __name__ == '__main__': + unittest.main() diff --git a/skais/utils/__pycache__/__init__.cpython-38.pyc b/skais/utils/__pycache__/__init__.cpython-38.pyc deleted file mode 100644 index a9577174079f5508328c6cc2288ec7bc1d22e257..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 167 zcmWIL<>g{vU|_JY%SZsxk3j@7W@KPsaA06yC}v?`U`SyIX3%8xTggxa5=IceeDpK& zb5r$;5(_dCQ*-n~GE$2(i}edKQsN5}i*gg=i&GPe@=Hrni}jQ9Q&RN<DoZl*^Yn|e p6ElnTOG`3yiuL2;GxIV_;^XxSDsOSv<mRW8=A_zzZ21gw4FH}9DbxS} diff --git a/skais/utils/__pycache__/geography.cpython-38.pyc b/skais/utils/__pycache__/geography.cpython-38.pyc deleted file mode 100644 index 4c93677db8d3da9f08059537b93a7a25c1f32d7d..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 746 zcmWIL<>g{vU|`_ZN=}qvVqka-;vi#Y1_lNP1_p*=GX@5R6owSW9EM!RC`LvQn<<AW ziaCWLg*k^Mmo<u&ks*~eg(ZcxmobGcg}sHLnK6njg(;Xplj9}GL`}w9!WoHWsYS(^ zd8zR!nZ+fEdC958ewvK8n6olVLNpm~apdI}RF-7q=Owd%bU-m50|NsG0|SFI$QfH0 z7#K<z7BJQ@EM%-@0kfDEGS;$!S&R!AYuUgo=7o&4>?w@395rk;ED{X0oHgt<tP%{h zTqO((m})pvn4}pNG9t4@K&CR+uz*FGYuFYtf<>5WxFBlLdEyMU+$l^Y3=2SRV6EW> z>tjt}uHlwoSjf1LNra)6r-q}3M}(o4w}v;EL6gNVbB9?X!%HRx28Nf63=9mKY`2*5 z3U09%C8lI1<`v&!F3!xm#hjd9e2YD?D7hpt&*&CQabZ!3Cf_ZVoWv4CFlls)H77qW z-4H?<-D0XRyv0;$c#EmR=oV9@(JjWHTkI*RIVFkl6);BSEyl!KjLEkcQ&ut*Ni#4o z{PNY$$j?pHFG?)PNKDPq56MU^&Mek1$ViDVOf1Sxj4w`2EXpq}NiEh-&QD3z4**Ar zesOkUX0d*0NoG#5etK$tI!Jq^UP0w8p7f&B#FF^r%%bF+R8~+V$$?^yiHi{onOK-O z7&(}97<v9PG4lLpV^(0|U@Vdag#de$L~(XbYGP5IUT$hhQD$<nUO{4JQF&%@Y7r=$ zZn5T-<`z^Iu`@6*fZ0ijMW6%;)(!Rx5+MUp#9@=0pHiBWY6lA8Vo<8#VdP-sVd7u} E0ICYU>i_@% diff --git a/skais/utils/geography.py b/skais/utils/geography.py index 82c5fa4..a45ec84 100644 --- a/skais/utils/geography.py +++ b/skais/utils/geography.py @@ -2,6 +2,9 @@ from sklearn.metrics.pairwise import haversine_distances import numpy as np from numba import jit +R = 6371000 + + @jit(nopython=True) def great_circle(lat1, lat2, long1, long2): x1 = np.radians(lat1) @@ -9,10 +12,8 @@ def great_circle(lat1, lat2, long1, long2): x2 = np.radians(lat2) y2 = np.radians(long2) - R = 6371000 - delta_x = x2 - x1 - delta_y= y2 -y1 + delta_y = y2 - y1 a = np.sin(delta_x / 2) * np.sin(delta_x / 2) + np.cos(x1) * np.cos(x2) * \ np.sin(delta_y / 2) * np.sin(delta_y / 2) @@ -21,3 +22,16 @@ def great_circle(lat1, lat2, long1, long2): d = R * c return d + + +def position_from_distance(position, distances): + lat = np.arcsin( + np.sin(np.radians(position[0])) * np.cos(distances[0] / R) + np.cos(np.radians(position[0])) * np.sin( + distances[0] / R)) + + long = np.radians(position[1]) + np.arctan2( + np.sin(distances[1] / R) * np.cos(np.radians(position[0])), + np.cos(distances[1] / R) - (np.sin(np.radians(position[1])) * np.sin(lat)) + ) + + return np.degrees(lat), np.degrees(long) -- GitLab