diff --git a/doc/source/api.rst b/doc/source/api.rst index c00804b2..04a803ac 100644 --- a/doc/source/api.rst +++ b/doc/source/api.rst @@ -73,7 +73,7 @@ Computation of angles :members: :undoc-members: -The azimuth angles are calculated using get_alt_az and get_observer_look from pyorbital. In pyorbital documentation there is a link describing how calculations are done http://celestrak.com/columns/v02n02/. The azimuth described in the link is measured as clockwise from North instead of counter-clockwise from South. Counter clockwise from south would be the standard for a right-handed orthogonal coordinate system. Pygac was updated to use the same definition for angles as pyorbital (20190926, version 1.1.0). Previous versions used azimuth +/-180 degrees, which correspond to degrees clockwise from south. All angles are converted to degrees. +The azimuth angles are calculated using get_alt_az and get_observer_look from pyorbital. In pyorbital documentation there is a link describing how calculations are done http://celestrak.com/columns/v02n02/. The azimuth described in the link is measured as clockwise from North instead of counter-clockwise from South. Counter clockwise from south would be the standard for a right-handed orthogonal coordinate system. Pygac was updated to use the same definition for angles as pyorbital (2019, September, version > 1.1.0). Previous versions used azimuth +/-180 degrees, which correspond to degrees clockwise from south. All angles are converted to degrees. All azimuth angles are converted to range ]-180, 180] (2019 October version > 1.1.0). Note that ]-180, 180] is an open interval. GAC calibration/inter-calibration diff --git a/pygac/__init__.py b/pygac/__init__.py index 815dad89..92259812 100644 --- a/pygac/__init__.py +++ b/pygac/__init__.py @@ -45,3 +45,10 @@ def get_absolute_azimuth_angle_diff(sat_azi, sun_azi): # Not using np.where to avoid copying array rel_azi[rel_azi > 180] = 360.0 - rel_azi[rel_azi > 180] return rel_azi + + +def centered_modulus(array, divisor): + """Transform array to half open range ]-divisor/2, divisor/2].""" + arr = array % divisor + arr[arr > divisor / 2] -= divisor + return arr diff --git a/pygac/gac_reader.py b/pygac/gac_reader.py index 43ac81d0..a64cbe75 100644 --- a/pygac/gac_reader.py +++ b/pygac/gac_reader.py @@ -30,7 +30,8 @@ import six import types -from pygac import CONFIG_FILE, get_absolute_azimuth_angle_diff +from pygac import (CONFIG_FILE, centered_modulus, + get_absolute_azimuth_angle_diff) try: import ConfigParser except ImportError: @@ -452,6 +453,23 @@ def get_tle_lines(self): return tle1, tle2 def get_angles(self): + """Get azimuth and zenith angles. + + Azimuth angle definition is the same as in pyorbital, but with + different units (degrees not radians for sun azimuth angles) + and different ranges. + + Returns: + sat_azi: satellite azimuth angle + degree clockwise from north in range ]-180, 180], + sat_zentih: satellite zenith angles in degrees in range [0,90], + sun_azi: sun azimuth angle + degree clockwise from north in range ]-180, 180], + sun_zentih: sun zenith angles in degrees in range [0,90], + rel_azi: absolute azimuth angle difference in degrees between sun + and sensor in range [0, 180] + + """ self.get_times() self.get_lonlat() tle1, tle2 = self.get_tle_lines() @@ -472,6 +490,10 @@ def get_angles(self): sun_azi = np.rad2deg(sun_azi) rel_azi = get_absolute_azimuth_angle_diff(sun_azi, sat_azi) + # Scale angles range to half open interval ]-180, 180] + sat_azi = centered_modulus(sat_azi, 360.0) + sun_azi = centered_modulus(sun_azi, 360.0) + # Mask corrupt scanlines for arr in (sat_azi, sat_zenith, sun_azi, sun_zenith, rel_azi): arr[self.mask] = np.nan diff --git a/pygac/tests/test_angles.py b/pygac/tests/test_angles.py index 3fc0b1c0..ec1a0a95 100644 --- a/pygac/tests/test_angles.py +++ b/pygac/tests/test_angles.py @@ -27,7 +27,8 @@ import numpy as np -from pygac import get_absolute_azimuth_angle_diff +from pygac import (get_absolute_azimuth_angle_diff, + centered_modulus) class TestAngles(unittest.TestCase): @@ -53,10 +54,24 @@ def test_azidiff_angle(self): np.testing.assert_allclose(rel_azi, res) + def test_centered_modulus(self): + """Test centered_modulus.""" + angles = np.ma.array( + [[180.0, -180.0, -179.9, 201.0], + [80.0, 360.0, -360.0, 604.0], + [-80.0, -88.0, -796.0, -104.0], + [-3.0, -188.0, -196.0, -204.0]], mask=False) + expected = np.ma.array( + [[180.0, 180.0, -179.9, -159.0], + [80.0, 0.0, 0.0, -116.0], + [-80.0, -88.0, -76.0, -104.0], + [-3.0, 172.0, 164.0, 156.0]], mask=False) + transformed = centered_modulus(angles, 360.0) + np.testing.assert_allclose(transformed, expected) + def suite(): - """The suite for test_slerp - """ + """Test suite for test_angles.""" loader = unittest.TestLoader() mysuite = unittest.TestSuite() mysuite.addTest(loader.loadTestsFromTestCase(TestAngles)) diff --git a/pygac/tests/test_reader.py b/pygac/tests/test_reader.py index 51bf3d30..6fea7750 100644 --- a/pygac/tests/test_reader.py +++ b/pygac/tests/test_reader.py @@ -65,7 +65,6 @@ def test_get_lonlat(self, interpolator, adjust_clockdrift, lat_i = np.array([1, 2, 3, np.nan, -90.1, 90.1]) get_lonlat.return_value = lon_i, lat_i interpolator.return_value = lon_i, lat_i - get_corrupt_mask.return_value = np.array( [0, 0, 1, 0, 0, 0], dtype=bool) @@ -192,9 +191,50 @@ def test_get_tle_lines(self, read_tle_file, *mocks): self.assertEqual(tle1, tle_data[tle_idx]) self.assertEqual(tle2, tle_data[tle_idx + 1]) + @mock.patch('pygac.gac_reader.GACReader._get_corrupt_mask') + def test_get_angles(self, get_corrupt_mask): + """Test get_angles function of the reader.""" + # Line: 1, 649, 6198 and 12658 from Tiros-N file (1980-01-03 11:47) + lon_i = np.array( + [69.41555, 152.10587, 164.3131, 67.23855, np.nan])[:, np.newaxis] + lat_i = np.array( + [71.6283, 85.24265, -62.076958, 82.72296, np.nan])[:, np.newaxis] + get_corrupt_mask.return_value = np.isnan(lon_i) + self.reader.lons = lon_i + self.reader.lats = lat_i + self.reader.tle_lines = [ + '1 11060U 78096A 80003.54792075 .00000937 00000-0 52481-3 0 2588\r\n', # noqa + '2 11060 98.9783 332.1605 0012789 88.8047 271.4583 14.11682873 63073\r\n'] # noqa + self.reader.utcs = np.array( + [315748035469, 315748359969, + 315751135469, 315754371969, + 315754371969]).astype('datetime64[ms]') + self.reader.spacecrafts_orbital = {25: 'tiros n'} + self.reader.spacecraft_id = 25 + self.reader.times = self.reader.to_datetime(self.reader.utcs) + expected_sat_azi = np.array( + [-76.90, 11.08, 145.33, -50.01, np.nan])[:, np.newaxis] + expected_sun_azi = np.array( + [-120.36, -31.94, -173.51, -93.67, np.nan])[:, np.newaxis] + expected_sat_zenith = np.array( + [69.05, 69.04, 69.55, 69.07, np.nan])[:, np.newaxis] + expected_sun_zenith = np.array( + [104.30, 116.94, 94.86, 112.60, np.nan])[:, np.newaxis] + expected_rel_azi = np.array( + [43.45, 43.01, 41.16, 43.65, np.nan])[:, np.newaxis] + + retv = self.reader.get_angles() + (sat_azi, sat_zenith, sun_azi, sun_zenith, rel_azi) = retv + np.testing.assert_allclose(sat_azi, expected_sat_azi, atol=0.01) + np.testing.assert_allclose(sun_azi, expected_sun_azi, atol=0.01) + np.testing.assert_allclose(sat_zenith, expected_sat_zenith, atol=0.01) + np.testing.assert_allclose(sun_zenith, expected_sun_zenith, atol=0.01) + np.testing.assert_allclose(rel_azi, expected_rel_azi, atol=0.01) + @mock.patch('pygac.gac_reader.ConfigParser.ConfigParser.read') @mock.patch('pygac.gac_reader.ConfigParser.ConfigParser.items') def test_get_tle_file(self, items, *mocks): + """Test get_tle_file.""" # Use TLE name/dir from config file items.return_value = [('tledir', 'a'), ('tlename', 'b')] tle_file = self.reader.get_tle_file()