From 1efe6769df76f819de2639d02f696b7f2320496e Mon Sep 17 00:00:00 2001 From: "Nina.Hakansson" Date: Mon, 12 Apr 2021 21:04:49 +0200 Subject: [PATCH 1/9] Fix get_observer_look for satellite exactly at nadir. When the satellite is exactly above the observer the satellite elevation angle is 90 degrees the azimuth angle is undefined. Set asimuth angle to 0 when satellite elevation is 90 degrees. Also handle the case when top_z / rg_ is larger than 1.0 due to rounding. --- pyorbital/orbital.py | 36 +++++++++++++++++++++++++++++---- pyorbital/tests/test_orbital.py | 8 ++++---- 2 files changed, 36 insertions(+), 8 deletions(-) diff --git a/pyorbital/orbital.py b/pyorbital/orbital.py index 9852277e..9fac1d3e 100644 --- a/pyorbital/orbital.py +++ b/pyorbital/orbital.py @@ -106,6 +106,9 @@ def get_observer_look(sat_lon, sat_lat, sat_alt, utc_time, lon, lat, alt): (opos_x, opos_y, opos_z), (ovel_x, ovel_y, ovel_z) = \ astronomy.observer_position(utc_time, lon, lat, alt) + # When satellit is at nadir, elevation is 90 degrees and azimuth is undefined. + sat_at_nadir = np.logical_and(lon == sat_lon, lat == sat_lat) + lon = np.deg2rad(lon) lat = np.deg2rad(lat) @@ -140,14 +143,39 @@ def get_observer_look(sat_lon, sat_lat, sat_alt, utc_time, lon, lat, alt): az_data[np.where(top_s > 0)] += np.pi az_data[np.where(az_data < 0)] += 2 * np.pi + rg_ = np.sqrt(rx * rx + ry * ry + rz * rz) + + top_z_divided_by_rg_ = top_z / rg_ + el_ = np.arcsin(top_z_divided_by_rg_) + if has_xarray and isinstance(az_, xr.DataArray): - az_.data = az_data + el_data = el_.data else: - az_ = az_data + el_data = el_ - rg_ = np.sqrt(rx * rx + ry * ry + rz * rz) - el_ = np.arcsin(top_z / rg_) + # Due to rounding top_z can be larger than rg_ (when el_ ~ 90). + # And azimuth undefined when elevation is 90 degrees + if has_dask and isinstance(az_data, da.Array): + el_data = da.where(top_z_divided_by_rg_ > 1.0, np.pi/2, el_data) + az_data = da.where(el_data == np.pi / 2, 0, az_data) + else: + el_data[np.where(top_z_divided_by_rg_ > 1.0)] = np.pi/2 + az_data[np.where(el_data == np.pi/2)] = 0 + + # When lat/lon is equal to sat_lat/sat_lon el_ is 90 degrees and az_ is undefined. + if has_dask and isinstance(az_data, da.Array): + el_data = da.where(sat_at_nadir, np.pi/2, el_data) + az_data = da.where(sat_at_nadir, 0, az_data) + else: + el_data[np.where(sat_at_nadir)] = np.pi / 2 + az_data[np.where(sat_at_nadir)] = 0 + if has_xarray and isinstance(az_, xr.DataArray): + az_.data = az_data + el_.data = el_data + else: + el_ = el_data + az_ = az_data return np.rad2deg(az_), np.rad2deg(el_) diff --git a/pyorbital/tests/test_orbital.py b/pyorbital/tests/test_orbital.py index 472e964a..0abee9a1 100644 --- a/pyorbital/tests/test_orbital.py +++ b/pyorbital/tests/test_orbital.py @@ -220,15 +220,15 @@ class TestGetObserverLook(unittest.TestCase): def setUp(self): self.t = datetime(2018, 1, 1, 0, 0, 0) self.sat_lon = np.array([[-89.5, -89.4], [-89.3, -89.2]]) - self.sat_lat = np.array([[45.5, 45.4], [45.3, 45.2]]) + self.sat_lat = np.array([[45.5, 45.4], [45.3, 40.2]]) self.sat_alt = np.array([[35786, 35786], [35786, 35786]]) - self.lon = np.array([[-85.5, -85.4], [-85.3, -85.2]]) + self.lon = np.array([[-85.5, -85.4], [-85.3, -89.2]]) self.lat = np.array([[40.5, 40.4], [40.3, 40.2]]) self.alt = np.zeros((2, 2)) self.exp_azi = np.array([[331.00275902, 330.95954165], - [330.91642994, 330.87342384]]) + [330.91642994, 0]]) self.exp_elev = np.array([[83.18070976, 83.17788976], - [83.17507167, 83.1722555]]) + [83.17507167, 90]]) def test_basic_numpy(self): """Test with numpy array inputs""" From 1bf54ef762ebe1ed9cfd06f38383ed56bdfa7652 Mon Sep 17 00:00:00 2001 From: "Nina.Hakansson" Date: Thu, 20 May 2021 17:42:08 +0200 Subject: [PATCH 2/9] Handle problematic az_data for zero divided by zero This happens when satellite is exactly at observer nadir. The problem depends on rounding and varies between array types: numpy, dask, and xarray. Removed the check for sat_at_nadir as is not needed any longer. --- pyorbital/orbital.py | 16 +++++----------- 1 file changed, 5 insertions(+), 11 deletions(-) diff --git a/pyorbital/orbital.py b/pyorbital/orbital.py index 9fac1d3e..f174d550 100644 --- a/pyorbital/orbital.py +++ b/pyorbital/orbital.py @@ -106,9 +106,6 @@ def get_observer_look(sat_lon, sat_lat, sat_alt, utc_time, lon, lat, alt): (opos_x, opos_y, opos_z), (ovel_x, ovel_y, ovel_z) = \ astronomy.observer_position(utc_time, lon, lat, alt) - # When satellit is at nadir, elevation is 90 degrees and azimuth is undefined. - sat_at_nadir = np.logical_and(lon == sat_lon, lat == sat_lat) - lon = np.deg2rad(lon) lat = np.deg2rad(lat) @@ -136,10 +133,15 @@ def get_observer_look(sat_lon, sat_lat, sat_alt, utc_time, lon, lat, alt): else: az_data = az_ + # When elevation is 90 degrees both top_s and top_e are close to zero + # and azimuth angle is undefined. + undet_azi = np.logical_and(top_s == 0, top_e == 0) if has_dask and isinstance(az_data, da.Array): + az_data = da.where(undet_azi, 0, az_data) az_data = da.where(top_s > 0, az_data + np.pi, az_data) az_data = da.where(az_data < 0, az_data + 2 * np.pi, az_data) else: + az_data[np.where(undet_azi)] = 0 az_data[np.where(top_s > 0)] += np.pi az_data[np.where(az_data < 0)] += 2 * np.pi @@ -162,14 +164,6 @@ def get_observer_look(sat_lon, sat_lat, sat_alt, utc_time, lon, lat, alt): el_data[np.where(top_z_divided_by_rg_ > 1.0)] = np.pi/2 az_data[np.where(el_data == np.pi/2)] = 0 - # When lat/lon is equal to sat_lat/sat_lon el_ is 90 degrees and az_ is undefined. - if has_dask and isinstance(az_data, da.Array): - el_data = da.where(sat_at_nadir, np.pi/2, el_data) - az_data = da.where(sat_at_nadir, 0, az_data) - else: - el_data[np.where(sat_at_nadir)] = np.pi / 2 - az_data[np.where(sat_at_nadir)] = 0 - if has_xarray and isinstance(az_, xr.DataArray): az_.data = az_data el_.data = el_data From 8023923f9fd3642025d3840ef0f7d53113c76280 Mon Sep 17 00:00:00 2001 From: "Nina.Hakansson" Date: Fri, 21 May 2021 13:37:22 +0200 Subject: [PATCH 3/9] Extend test of satellite azimuth differense angle --- pyorbital/tests/test_orbital.py | 116 +++++++++++++++++++++++++++++--- 1 file changed, 106 insertions(+), 10 deletions(-) diff --git a/pyorbital/tests/test_orbital.py b/pyorbital/tests/test_orbital.py index 0abee9a1..5495128b 100644 --- a/pyorbital/tests/test_orbital.py +++ b/pyorbital/tests/test_orbital.py @@ -219,16 +219,20 @@ class TestGetObserverLook(unittest.TestCase): def setUp(self): self.t = datetime(2018, 1, 1, 0, 0, 0) - self.sat_lon = np.array([[-89.5, -89.4], [-89.3, -89.2]]) - self.sat_lat = np.array([[45.5, 45.4], [45.3, 40.2]]) - self.sat_alt = np.array([[35786, 35786], [35786, 35786]]) - self.lon = np.array([[-85.5, -85.4], [-85.3, -89.2]]) - self.lat = np.array([[40.5, 40.4], [40.3, 40.2]]) - self.alt = np.zeros((2, 2)) - self.exp_azi = np.array([[331.00275902, 330.95954165], - [330.91642994, 0]]) - self.exp_elev = np.array([[83.18070976, 83.17788976], - [83.17507167, 90]]) + self.sat_lon = np.array([[-89.5, -89.4, -89.5, -89.4], + [-89.3, -89.2, -89.3, -89.2]]) + self.sat_lat = np.array([[45.5, 45.4, 45.5, 45.4], + [45.3, 40.2, 45.3, 40.2]]) + self.sat_alt = 35786 * np.ones((2, 4)) + self.lon = np.array([[-85.5, -85.4, -89.5, -99.4], + [-85.3, -89.2, -89.3, -79.2]]) + self.lat = np.array([[40.5, 40.4, 65.5, 45.4], + [40.3, 40.2, 25.3, 40.2]]) + self.alt = np.zeros((2, 4)) + self.exp_azi = np.array([[331.00275902, 330.95954165, 180, 86.435411], + [330.91642994, 0, 0, 273.232073]]) + self.exp_elev = np.array([[83.18070976, 83.17788976, 66.548467, 81.735221], + [83.17507167, 90, 66.559906, 81.010018]]) def test_basic_numpy(self): """Test with numpy array inputs""" @@ -295,6 +299,97 @@ def _xarr_conv(input): np.testing.assert_allclose(elev.data.compute(), self.exp_elev) +class TestGetObserverLookNadir(unittest.TestCase): + """Test the get_observer_look function when satellite is at nadir.""" + + def setUp(self): + """Setup for test observer at nadir. + Note that rounding error differs between array types. + With 1000 elements a test gives: + 1 error for basic numpy + 41 errors for basic dask + 63 errors for xarray with dask + 2 error for xarray with numpy + """ + + self.t = datetime(2018, 1, 1, 0, 0, 0) + self.sat_lon = 360 * np.random.rand(1000) - 180 + self.sat_lat = 180 * np.random.rand(1000) - 90 + self.sat_alt = np.random.rand(1000) + 850 + self.lon = self.sat_lon # + 10E-17 + self.lat = self.sat_lat # + 10E-17 + self.alt = np.zeros((1000)) + self.exp_elev = np.zeros((1000)) + 90 + + def test_basic_numpy(self): + """Test with numpy array inputs""" + from pyorbital import orbital + azi, elev = orbital.get_observer_look(self.sat_lon, self.sat_lat, + self.sat_alt, self.t, + self.lon, self.lat, self.alt) + self.assertEqual(np.sum(np.isnan(azi)), 0) + self.assertFalse(np.isnan(azi).any()) + np.testing.assert_allclose(elev, self.exp_elev) + + def test_basic_dask(self): + """Test with dask array inputs""" + from pyorbital import orbital + import dask.array as da + sat_lon = da.from_array(self.sat_lon, chunks=2) + sat_lat = da.from_array(self.sat_lat, chunks=2) + sat_alt = da.from_array(self.sat_alt, chunks=2) + lon = da.from_array(self.lon, chunks=2) + lat = da.from_array(self.lat, chunks=2) + alt = da.from_array(self.alt, chunks=2) + azi, elev = orbital.get_observer_look(sat_lon, sat_lat, + sat_alt, self.t, + lon, lat, alt) + self.assertEqual(np.sum(np.isnan(azi)), 0) + self.assertFalse(np.isnan(azi).any()) + np.testing.assert_allclose(elev.compute(), self.exp_elev) + + def test_xarray_with_numpy(self): + """Test with xarray DataArray with numpy array as inputs""" + from pyorbital import orbital + import xarray as xr + + def _xarr_conv(input): + return xr.DataArray(input) + sat_lon = _xarr_conv(self.sat_lon) + sat_lat = _xarr_conv(self.sat_lat) + sat_alt = _xarr_conv(self.sat_alt) + lon = _xarr_conv(self.lon) + lat = _xarr_conv(self.lat) + alt = _xarr_conv(self.alt) + azi, elev = orbital.get_observer_look(sat_lon, sat_lat, + sat_alt, self.t, + lon, lat, alt) + self.assertEqual(np.sum(np.isnan(azi)), 0) + self.assertFalse(np.isnan(azi).any()) + np.testing.assert_allclose(elev.data, self.exp_elev) + + def test_xarray_with_dask(self): + """Test with xarray DataArray with dask array as inputs""" + from pyorbital import orbital + import dask.array as da + import xarray as xr + + def _xarr_conv(input): + return xr.DataArray(da.from_array(input, chunks=2)) + sat_lon = _xarr_conv(self.sat_lon) + sat_lat = _xarr_conv(self.sat_lat) + sat_alt = _xarr_conv(self.sat_alt) + lon = _xarr_conv(self.lon) + lat = _xarr_conv(self.lat) + alt = _xarr_conv(self.alt) + azi, elev = orbital.get_observer_look(sat_lon, sat_lat, + sat_alt, self.t, + lon, lat, alt) + self.assertEqual(np.sum(np.isnan(azi)), 0) + self.assertFalse(np.isnan(azi).any()) + np.testing.assert_allclose(elev.data.compute(), self.exp_elev) + + class TestRegressions(unittest.TestCase): """Test regressions.""" @@ -318,6 +413,7 @@ def suite(): mysuite = unittest.TestSuite() mysuite.addTest(loader.loadTestsFromTestCase(Test)) mysuite.addTest(loader.loadTestsFromTestCase(TestGetObserverLook)) + mysuite.addTest(loader.loadTestsFromTestCase(TestGetObserverLookNadir)) mysuite.addTest(loader.loadTestsFromTestCase(TestRegressions)) return mysuite From 89beb6dbf844e8174660f11e1c5e49f9e9a4be0a Mon Sep 17 00:00:00 2001 From: "Nina.Hakansson" Date: Fri, 21 May 2021 13:38:02 +0200 Subject: [PATCH 4/9] Simplyfy get_observer_look --- pyorbital/orbital.py | 42 +++++------------------------------------- 1 file changed, 5 insertions(+), 37 deletions(-) diff --git a/pyorbital/orbital.py b/pyorbital/orbital.py index f174d550..d1923a49 100644 --- a/pyorbital/orbital.py +++ b/pyorbital/orbital.py @@ -126,50 +126,18 @@ def get_observer_look(sat_lon, sat_lat, sat_alt, utc_time, lon, lat, alt): top_z = cos_lat * cos_theta * rx + \ cos_lat * sin_theta * ry + sin_lat * rz - az_ = np.arctan(-top_e / top_s) - - if has_xarray and isinstance(az_, xr.DataArray): - az_data = az_.data - else: - az_data = az_ - - # When elevation is 90 degrees both top_s and top_e are close to zero - # and azimuth angle is undefined. - undet_azi = np.logical_and(top_s == 0, top_e == 0) - if has_dask and isinstance(az_data, da.Array): - az_data = da.where(undet_azi, 0, az_data) - az_data = da.where(top_s > 0, az_data + np.pi, az_data) - az_data = da.where(az_data < 0, az_data + 2 * np.pi, az_data) - else: - az_data[np.where(undet_azi)] = 0 - az_data[np.where(top_s > 0)] += np.pi - az_data[np.where(az_data < 0)] += 2 * np.pi + az_ = np.arctan2(-top_e, top_s) + np.pi rg_ = np.sqrt(rx * rx + ry * ry + rz * rz) top_z_divided_by_rg_ = top_z / rg_ - el_ = np.arcsin(top_z_divided_by_rg_) - - if has_xarray and isinstance(az_, xr.DataArray): - el_data = el_.data - else: - el_data = el_ # Due to rounding top_z can be larger than rg_ (when el_ ~ 90). # And azimuth undefined when elevation is 90 degrees - if has_dask and isinstance(az_data, da.Array): - el_data = da.where(top_z_divided_by_rg_ > 1.0, np.pi/2, el_data) - az_data = da.where(el_data == np.pi / 2, 0, az_data) - else: - el_data[np.where(top_z_divided_by_rg_ > 1.0)] = np.pi/2 - az_data[np.where(el_data == np.pi/2)] = 0 - - if has_xarray and isinstance(az_, xr.DataArray): - az_.data = az_data - el_.data = el_data - else: - el_ = el_data - az_ = az_data + top_z_divided_by_rg_ = top_z_divided_by_rg_.clip(max=1) + el_ = np.arcsin(top_z_divided_by_rg_) + + az_ -= np.pi * (top_z_divided_by_rg_ == 1) return np.rad2deg(az_), np.rad2deg(el_) From 1efbd7594a7147905bdadacd2b273b2f6b4c9202 Mon Sep 17 00:00:00 2001 From: "Nina.Hakansson" Date: Fri, 21 May 2021 14:03:04 +0200 Subject: [PATCH 5/9] Change undefined azimuth value to 180 --- pyorbital/orbital.py | 1 - pyorbital/tests/test_orbital.py | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/pyorbital/orbital.py b/pyorbital/orbital.py index d1923a49..ba1a46de 100644 --- a/pyorbital/orbital.py +++ b/pyorbital/orbital.py @@ -137,7 +137,6 @@ def get_observer_look(sat_lon, sat_lat, sat_alt, utc_time, lon, lat, alt): top_z_divided_by_rg_ = top_z_divided_by_rg_.clip(max=1) el_ = np.arcsin(top_z_divided_by_rg_) - az_ -= np.pi * (top_z_divided_by_rg_ == 1) return np.rad2deg(az_), np.rad2deg(el_) diff --git a/pyorbital/tests/test_orbital.py b/pyorbital/tests/test_orbital.py index 5495128b..622d7bc8 100644 --- a/pyorbital/tests/test_orbital.py +++ b/pyorbital/tests/test_orbital.py @@ -230,7 +230,7 @@ def setUp(self): [40.3, 40.2, 25.3, 40.2]]) self.alt = np.zeros((2, 4)) self.exp_azi = np.array([[331.00275902, 330.95954165, 180, 86.435411], - [330.91642994, 0, 0, 273.232073]]) + [330.91642994, 180, 0, 273.232073]]) self.exp_elev = np.array([[83.18070976, 83.17788976, 66.548467, 81.735221], [83.17507167, 90, 66.559906, 81.010018]]) From 5e715109a45dd58272b1f9a25a9899bb6fcbbe4e Mon Sep 17 00:00:00 2001 From: "Nina.Hakansson" Date: Fri, 21 May 2021 14:43:40 +0200 Subject: [PATCH 6/9] Fix problem for macos (hopefully) --- pyorbital/orbital.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyorbital/orbital.py b/pyorbital/orbital.py index ba1a46de..9b9d0f73 100644 --- a/pyorbital/orbital.py +++ b/pyorbital/orbital.py @@ -126,14 +126,15 @@ def get_observer_look(sat_lon, sat_lat, sat_alt, utc_time, lon, lat, alt): top_z = cos_lat * cos_theta * rx + \ cos_lat * sin_theta * ry + sin_lat * rz + # Azimuth is undefined when elevation is 90 degrees, 180 (pi) will be returned. az_ = np.arctan2(-top_e, top_s) + np.pi + az_ = np.mod(az_, 2 * np.pi) # Needed on some platforms rg_ = np.sqrt(rx * rx + ry * ry + rz * rz) top_z_divided_by_rg_ = top_z / rg_ # Due to rounding top_z can be larger than rg_ (when el_ ~ 90). - # And azimuth undefined when elevation is 90 degrees top_z_divided_by_rg_ = top_z_divided_by_rg_.clip(max=1) el_ = np.arcsin(top_z_divided_by_rg_) From a885dcea07239d66f44c968c6fe52aebcce2491b Mon Sep 17 00:00:00 2001 From: "Nina.Hakansson" Date: Fri, 21 May 2021 15:29:38 +0200 Subject: [PATCH 7/9] Fix deepcode issue and reduce number of elements for test --- pyorbital/tests/test_orbital.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/pyorbital/tests/test_orbital.py b/pyorbital/tests/test_orbital.py index 622d7bc8..b08701bb 100644 --- a/pyorbital/tests/test_orbital.py +++ b/pyorbital/tests/test_orbital.py @@ -312,14 +312,15 @@ def setUp(self): 2 error for xarray with numpy """ + np.random.seed(125) self.t = datetime(2018, 1, 1, 0, 0, 0) - self.sat_lon = 360 * np.random.rand(1000) - 180 - self.sat_lat = 180 * np.random.rand(1000) - 90 - self.sat_alt = np.random.rand(1000) + 850 + self.sat_lon = 360 * np.random.rand(100) - 180 + self.sat_lat = 180 * np.random.rand(100) - 90 + self.sat_alt = np.random.rand(100) + 850 self.lon = self.sat_lon # + 10E-17 self.lat = self.sat_lat # + 10E-17 - self.alt = np.zeros((1000)) - self.exp_elev = np.zeros((1000)) + 90 + self.alt = np.zeros((100)) + self.exp_elev = np.zeros((100)) + 90 def test_basic_numpy(self): """Test with numpy array inputs""" From 7099f8199f372fbcb6c2d2fda00a6d3f135a2f83 Mon Sep 17 00:00:00 2001 From: "Nina.Hakansson" Date: Mon, 24 May 2021 12:12:27 +0200 Subject: [PATCH 8/9] Attempt two to fix deepcode issue with random numbers --- pyorbital/tests/test_orbital.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/pyorbital/tests/test_orbital.py b/pyorbital/tests/test_orbital.py index b08701bb..511bbbf0 100644 --- a/pyorbital/tests/test_orbital.py +++ b/pyorbital/tests/test_orbital.py @@ -311,11 +311,10 @@ def setUp(self): 63 errors for xarray with dask 2 error for xarray with numpy """ - - np.random.seed(125) + rng = np.random.RandomState(125) self.t = datetime(2018, 1, 1, 0, 0, 0) - self.sat_lon = 360 * np.random.rand(100) - 180 - self.sat_lat = 180 * np.random.rand(100) - 90 + self.sat_lon = 360 * rng.rand(100) - 180 + self.sat_lat = 180 * rng.rand(100) - 90 self.sat_alt = np.random.rand(100) + 850 self.lon = self.sat_lon # + 10E-17 self.lat = self.sat_lat # + 10E-17 From ca1d4d20ebc554174b7d74e943f988a225c04cd9 Mon Sep 17 00:00:00 2001 From: "Nina.Hakansson" Date: Mon, 24 May 2021 13:08:38 +0200 Subject: [PATCH 9/9] use the random number generator for all calls (forgot one) --- pyorbital/tests/test_orbital.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyorbital/tests/test_orbital.py b/pyorbital/tests/test_orbital.py index 511bbbf0..af974cff 100644 --- a/pyorbital/tests/test_orbital.py +++ b/pyorbital/tests/test_orbital.py @@ -315,7 +315,7 @@ def setUp(self): self.t = datetime(2018, 1, 1, 0, 0, 0) self.sat_lon = 360 * rng.rand(100) - 180 self.sat_lat = 180 * rng.rand(100) - 90 - self.sat_alt = np.random.rand(100) + 850 + self.sat_alt = rng.rand(100) + 850 self.lon = self.sat_lon # + 10E-17 self.lat = self.sat_lat # + 10E-17 self.alt = np.zeros((100))