diff --git a/wetb/utils/geometry.py b/wetb/utils/geometry.py index e82594ab3f77f7c8e43963006ec834553e223c84..3a7761cbd2bdb34add9c223b0df269a1fba9f18f 100644 --- a/wetb/utils/geometry.py +++ b/wetb/utils/geometry.py @@ -29,13 +29,33 @@ def mean_deg(dir, axis=0): ---------- dir : array_like Angles in degrees + axis : int + if dir is 2d array_like, axis defines which axis to take the mean of Returns ------- mean_deg : float Mean angle """ - return deg(np.arctan2(np.mean(sind(dir[:]), axis), np.mean(cosd(dir[:]), axis))) + return deg(mean_rad(rad(dir), axis)) + +def mean_rad(dir, axis=0): + """Mean of angles in radians + + Parameters + ---------- + dir : array_like + Angles in radians + axis : int + if dir is 2d array_like, axis defines which axis to take the mean of + + Returns + ------- + mean_rad : float + Mean angle + """ + return np.arctan2(np.nanmean(np.sin(dir[:]), axis), np.nanmean(np.cos(dir[:]), axis)) + def std_deg(dir): """Standard deviation of angles in degrees @@ -50,8 +70,22 @@ def std_deg(dir): std_deg : float standard deviation """ - return deg(np.sqrt(1 - (np.mean(sind(dir)) ** 2 + np.mean(cosd(dir)) ** 2))) + return deg(std_rad(rad(dir))) + +def std_rad(dir): + """Standard deviation of angles in radians + + Parameters + ---------- + dir : array_like + Angles in radians + Returns + ------- + std_rad : float + standard deviation + """ + return np.sqrt(1 - (np.nanmean(np.sin(dir)) ** 2 + np.nanmean(np.cos(dir)) ** 2)) def wsp_dir2uv(wsp, dir, dir_ref=None): """Convert horizontal wind speed and direction to u,v diff --git a/wetb/utils/tests/test_geometry.py b/wetb/utils/tests/test_geometry.py index e2fbf1f1faa5d951d611b4386cf0a95ca800ab42..1ae82832e36ab0d269a26a4aea16d03e818a2520 100644 --- a/wetb/utils/tests/test_geometry.py +++ b/wetb/utils/tests/test_geometry.py @@ -57,11 +57,17 @@ class TestGeometry(unittest.TestCase): np.testing.assert_array_almost_equal(mean_deg(a, 1), [45, 0, -45]) np.testing.assert_array_almost_equal(mean_deg(a.T, 0), [45, 0, -45]) + def test_mean_deg_nan(self): + self.assertEqual(mean_deg(np.array([0., 90, np.nan])), 45) + def test_std_deg(self): self.assertEqual(std_deg(np.array([0, 0, 0])), 0) self.assertAlmostEqual(std_deg(np.array([0, 90, 180, 270])), 57.296, 2) + def test_std_deg_nan(self): + self.assertAlmostEqual(std_deg(np.array([0, 90, 180, 270, np.nan])), 57.296, 2) + def test_wspdir2uv(self): u, v = wsp_dir2uv(np.array([1, 1, 1]), np.array([30, 0, 330])) np.testing.assert_array_almost_equal(u, [0.8660, 1, 0.8660], 3)