Clarify documentation of orbital_parameters metadata

doc/source/conf.py

@@ -249,16 +249,17 @@ def __getattr__(cls, name):
 
 # Example configuration for intersphinx: refer to the Python standard library.
 intersphinx_mapping = {
-    'python': ('https://docs.python.org/3', None),
-    'numpy': ('https://docs.scipy.org/doc/numpy', None),
-    'scipy': ('https://docs.scipy.org/doc/scipy/reference', None),
-    'xarray': ('https://xarray.pydata.org/en/stable', None),
     'dask': ('https://docs.dask.org/en/latest', None),
+    'geoviews': ('http://geoviews.org', None),
     'jobqueue': ('https://jobqueue.dask.org/en/latest', None),
+    'numpy': ('https://docs.scipy.org/doc/numpy', None),
+    'pydecorate': ('https://pydecorate.readthedocs.io/en/stable', None),
+    'pyorbital': ('https://pyorbital.readthedocs.io/en/stable', None),
+    'pyproj': ('https://pyproj4.github.io/pyproj/dev', None),
     'pyresample': ('https://pyresample.readthedocs.io/en/stable', None),
-    'trollsift': ('https://trollsift.readthedocs.io/en/stable', None),
+    'python': ('https://docs.python.org/3', None),
+    'scipy': ('https://docs.scipy.org/doc/scipy/reference', None),
     'trollimage': ('https://trollimage.readthedocs.io/en/stable', None),
-    'pydecorate': ('https://pydecorate.readthedocs.io/en/stable', None),
-    'geoviews': ('http://geoviews.org', None),
-    'pyproj': ('https://pyproj4.github.io/pyproj/dev', None)
+    'trollsift': ('https://trollsift.readthedocs.io/en/stable', None),
+    'xarray': ('https://xarray.pydata.org/en/stable', None)
 }

doc/source/dev_guide/custom_reader.rst

@@ -391,8 +391,8 @@ needs to implement a few methods:
    - the dataset info that is the description of the channel in the YAML file
 
   This method has to return an xarray.DataArray instance if the loading is
-  successful, containing the data and metadata of the loaded dataset, or
-  return None if the loading was unsuccessful.
+  successful, containing the data and :ref:`metadata <dataset_metadata>` of the
+  loaded dataset, or return None if the loading was unsuccessful.
 
  - the ``get_area_def`` method, that takes as single argument the dataset ID for which we want
    the area. For the data that cannot be geolocated with an area

doc/source/readers.rst

@@ -112,37 +112,46 @@ more information on the possible parameters.
 Metadata
 ========
 
-The datasets held by a scene also provide vital metadata such as dataset name, units, observation time etc. The
-following attributes are standardized across all readers:
+.. _dataset_metadata:
 
-* ``name``, ``wavelength``, ``resolution``, ``polarization``, ``calibration``, ``level``, ``modifiers``: See
-  :class:`satpy.dataset.DatasetID`.
+The datasets held by a scene also provide vital metadata such as dataset name, units, observation
+time etc. The following attributes are standardized across all readers:
+
+* ``name``, ``wavelength``, ``resolution``, ``polarization``, ``calibration``, ``level``,
+  ``modifiers``: See :class:`satpy.dataset.DatasetID`.
 * ``start_time``: Left boundary of the time interval covered by the dataset.
 * ``end_time``: Right boundary of the time interval covered by the dataset.
-* ``area``: :class:`~pyresample.geometry.AreaDefinition` or :class:`~pyresample.geometry.SwathDefinition` if
-  if data is geolocated. Areas are used for gridded projected data and Swaths when data must be
-  described by individual longitude/latitude coordinates. See the Coordinates section below.
+* ``area``: :class:`~pyresample.geometry.AreaDefinition` or
+  :class:`~pyresample.geometry.SwathDefinition` if data is geolocated. Areas are used for gridded
+  projected data and Swaths when data must be described by individual longitude/latitude
+  coordinates. See the Coordinates section below.
 * ``orbital_parameters``: Dictionary of orbital parameters describing the satellite's position.
 
   * For *geostationary* satellites it is described using the following scalar attributes:
 
-    * ``satellite_actual_longitude/latitude/altitude``: Current position of the satellite at the time of observation in
-      geodetic coordinates (i.e. altitude is normal to the surface).
-    * ``satellite_nominal_longitude/latitude/altitude``: Centre of the station keeping box (a confined area in which
-      the satellite is actively maintained in using maneuvres). Inbetween major maneuvres, when the satellite
-      is permanently moved, the nominal position is constant.
-    * ``nadir_longitude/latitude``: Intersection of the instrument's Nadir with the surface of the earth. May differ
-      from the actual satellite position, if the instrument is poiting slightly off the axis (satellite, earth-centre).
-      If available, this should be used to compute viewing angles etc. Otherwise, use the actual satellite position.
-    * ``projection_longitude/latitude/altitude``: Projection centre of the re-projected data. This should be used to
-      compute lat/lon coordinates. Note that the projection centre can differ considerably from the actual satellite
-      position. For example MSG-1 was at times positioned at 3.4 degrees west, while the image data was re-projected
-      to 0 degrees.
-    * [DEPRECATED] ``satellite_longitude/latitude/altitude``: Current position of the satellite at the time of observation
-      in geodetic coordinates.
-
-  * For *polar orbiting* satellites the readers usually provide coordinates and viewing angles of the swath as
-    ancillary datasets. Additional metadata related to the satellite position include:
+    * ``satellite_actual_longitude/latitude/altitude``: Current position of the satellite at the
+      time of observation in geodetic coordinates (i.e. altitude is relative and normal to the
+      surface of the ellipsoid).
+    * ``satellite_nominal_longitude/latitude/altitude``: Center of the station keeping box (a
+      confined area in which the satellite is actively maintained in using maneuvres). Inbetween
+      major maneuvres, when the satellite is permanently moved, the nominal position is constant.
+    * ``nadir_longitude/latitude``: Intersection of the instrument's Nadir with the surface of the
+      earth. May differ from the actual satellite position, if the instrument is pointing slightly
+      off the axis (satellite, earth-center). If available, this should be used to compute viewing
+      angles etc. Otherwise, use the actual satellite position.
+    * ``projection_longitude/latitude/altitude``: Projection center of the re-projected data. This
+      should be used to compute lat/lon coordinates. Note that the projection center can differ
+      considerably from the actual satellite position. For example MSG-1 was at times positioned
+      at 3.4 degrees west, while the image data was re-projected to 0 degrees.
+    * [DEPRECATED] ``satellite_longitude/latitude/altitude``: Current position of the satellite at
+      the time of observation in geodetic coordinates.
+
+    .. note:: Longitudes and latitudes are given in degrees, altitude in meters. For use in
+              pyorbital, the altitude has to be converted to kilometers, see for example
+              :func:`pyorbital.orbital.get_observer_look`.
+
+  * For *polar orbiting* satellites the readers usually provide coordinates and viewing angles of
+    the swath as ancillary datasets. Additional metadata related to the satellite position include:
 
       * ``tle``: Two-Line Element (TLE) set used to compute the satellite's orbit
 

satpy/composites/__init__.py

@@ -508,7 +508,7 @@ class PSPRayleighReflectance(CompositeBase):
     _rayleigh_cache = WeakValueDictionary()
 
     def get_angles(self, vis):
-        """Get the sun and satellite angles fro the current dataarray."""
+        """Get the sun and satellite angles from the current dataarray."""
         from pyorbital.astronomy import get_alt_az, sun_zenith_angle
         from pyorbital.orbital import get_observer_look
 
@@ -521,7 +521,7 @@ def get_angles(self, vis):
         sata, satel = get_observer_look(
             sat_lon,
             sat_lat,
-            sat_alt,
+            sat_alt / 1000.0,  # km
             vis.attrs['start_time'],
             lons, lats, 0)
         satz = 90 - satel
@@ -685,7 +685,7 @@ def __call__(self, projectables, optional_datasets=None, **info):
             try:
                 dummy, satel = get_observer_look(sat_lon,
                                                  sat_lat,
-                                                 sat_alt,
+                                                 sat_alt / 1000.0,  # km
                                                  band.attrs['start_time'],
                                                  lons, lats, 0)
             except KeyError:

satpy/composites/viirs.py

@@ -29,6 +29,7 @@
 from satpy.composites import CompositeBase, GenericCompositor
 from satpy.config import get_environ_ancpath
 from satpy.dataset import combine_metadata
+from satpy.utils import get_satpos
 
 LOG = logging.getLogger(__name__)
 
@@ -156,10 +157,11 @@ def get_angles(self, vis):
         suna = get_alt_az(vis.attrs['start_time'], lons, lats)[1]
         suna = np.rad2deg(suna)
         sunz = sun_zenith_angle(vis.attrs['start_time'], lons, lats)
+        sat_lon, sat_lat, sat_alt = get_satpos(vis)
         sata, satel = get_observer_look(
-            vis.attrs['satellite_longitude'],
-            vis.attrs['satellite_latitude'],
-            vis.attrs['satellite_altitude'],
+            sat_lon,
+            sat_lat,
+            sat_alt / 1000.0,  # km
             vis.attrs['start_time'],
             lons, lats, 0)
         satz = 90 - satel

satpy/readers/abi_l1b.py

@@ -60,7 +60,7 @@ def get_dataset(self, key, info):
                           'sensor': self.sensor,
                           'satellite_latitude': float(self['nominal_satellite_subpoint_lat']),
                           'satellite_longitude': float(self['nominal_satellite_subpoint_lon']),
-                          'satellite_altitude': float(self['nominal_satellite_height'])})
+                          'satellite_altitude': float(self['nominal_satellite_height']) * 1000.})
 
         # Add orbital parameters
         projection = self.nc["goes_imager_projection"]
@@ -70,7 +70,7 @@ def get_dataset(self, key, info):
             'projection_altitude': float(projection.attrs['perspective_point_height']),
             'satellite_nominal_latitude': float(self['nominal_satellite_subpoint_lat']),
             'satellite_nominal_longitude': float(self['nominal_satellite_subpoint_lon']),
-            'satellite_nominal_altitude': float(self['nominal_satellite_height']),
+            'satellite_nominal_altitude': float(self['nominal_satellite_height']) * 1000.,
             'yaw_flip': bool(self['yaw_flip_flag']),
         }
 

satpy/readers/ami_l1b.py

@@ -137,7 +137,7 @@ def get_orbital_parameters(self):
             'projection_altitude': h,
             'satellite_actual_longitude': sc_position[0],
             'satellite_actual_latitude': sc_position[1],
-            'satellite_actual_altitude': sc_position[2] / 1000.0,  # km
+            'satellite_actual_altitude': sc_position[2],  # meters
         }
         return orbital_parameters
 

satpy/tests/compositor_tests/__init__.py

@@ -1066,6 +1066,88 @@ def test_multiple_sensors(self):
         self.assertEqual(res.attrs['sensor'], {'abi', 'glm'})
 
 
+class TestPSPAtmosphericalCorrection(unittest.TestCase):
+    def setUp(self):
+        """Patch in-class imports."""
+        self.orbital = mock.MagicMock()
+        modules = {
+            'pyspectral.atm_correction_ir': mock.MagicMock(),
+            'pyorbital.orbital': self.orbital,
+        }
+        self.module_patcher = mock.patch.dict('sys.modules', modules)
+        self.module_patcher.start()
+
+    def tearDown(self):
+        """Unpatch in-class imports."""
+        self.module_patcher.stop()
+
+    @mock.patch('satpy.composites.PSPAtmosphericalCorrection.apply_modifier_info')
+    @mock.patch('satpy.composites.get_satpos')
+    def test_call(self, get_satpos, *mocks):
+        """Test atmospherical correction."""
+        from satpy.composites import PSPAtmosphericalCorrection
+
+        # Patch methods
+        get_satpos.return_value = 'sat_lon', 'sat_lat', 12345678
+        self.orbital.get_observer_look.return_value = 0, 0
+        area = mock.MagicMock()
+        area.get_lonlats.return_value = 'lons', 'lats'
+        band = mock.MagicMock(attrs={'area': area,
+                                     'start_time': 'start_time',
+                                     'name': 'name',
+                                     'platform_name': 'platform',
+                                     'sensor': 'sensor'})
+
+        # Perform atmospherical correction
+        psp = PSPAtmosphericalCorrection(name='dummy')
+        psp(projectables=[band])
+
+        # Check arguments of get_orbserver_look() call, especially the altitude
+        # unit conversion from meters to kilometers
+        self.orbital.get_observer_look.assert_called_with(
+            'sat_lon', 'sat_lat', 12345.678, 'start_time', 'lons', 'lats', 0)
+
+
+class TestPSPRayleighReflectance(unittest.TestCase):
+    def setUp(self):
+        """Patch in-class imports."""
+        self.astronomy = mock.MagicMock()
+        self.orbital = mock.MagicMock()
+        modules = {
+            'pyorbital.astronomy': self.astronomy,
+            'pyorbital.orbital': self.orbital,
+        }
+        self.module_patcher = mock.patch.dict('sys.modules', modules)
+        self.module_patcher.start()
+
+    def tearDown(self):
+        """Unpatch in-class imports."""
+        self.module_patcher.stop()
+
+    @mock.patch('satpy.composites.get_satpos')
+    def test_get_angles(self, get_satpos):
+        """Test sun and satellite angle calculation."""
+        from satpy.composites import PSPRayleighReflectance
+
+        # Patch methods
+        get_satpos.return_value = 'sat_lon', 'sat_lat', 12345678
+        self.orbital.get_observer_look.return_value = 0, 0
+        self.astronomy.get_alt_az.return_value = 0, 0
+        area = mock.MagicMock()
+        area.get_lonlats.return_value = 'lons', 'lats'
+        vis = mock.MagicMock(attrs={'area': area,
+                                    'start_time': 'start_time'})
+
+        # Compute angles
+        psp = PSPRayleighReflectance(name='dummy')
+        psp.get_angles(vis)
+
+        # Check arguments of get_orbserver_look() call, especially the altitude
+        # unit conversion from meters to kilometers
+        self.orbital.get_observer_look.assert_called_with(
+            'sat_lon', 'sat_lat', 12345.678, 'start_time', 'lons', 'lats', 0)
+
+
 def suite():
     """Test suite for all reader tests."""
     loader = unittest.TestLoader()
@@ -1092,6 +1174,8 @@ def suite():
     mysuite.addTest(loader.loadTestsFromTestCase(TestAddBands))
     mysuite.addTest(loader.loadTestsFromTestCase(TestBackgroundCompositor))
     mysuite.addTest(loader.loadTestsFromTestCase(TestStaticImageCompositor))
+    mysuite.addTest(loader.loadTestsFromTestCase(TestPSPAtmosphericalCorrection))
+    mysuite.addTest(loader.loadTestsFromTestCase(TestPSPRayleighReflectance))
 
     return mysuite
 

satpy/tests/compositor_tests/test_viirs.py

@@ -23,6 +23,10 @@
     import unittest2 as unittest
 else:
     import unittest
+try:
+    from unittest import mock
+except ImportError:
+    import mock
 
 
 class TestVIIRSComposites(unittest.TestCase):
@@ -282,7 +286,7 @@ def test_reflectance_corrector_abi(self):
         c01 = xr.DataArray(dnb,
                            dims=('y', 'x'),
                            attrs={'satellite_longitude': -89.5, 'satellite_latitude': 0.0,
-                                  'satellite_altitude': 35786.0234375, 'platform_name': 'GOES-16',
+                                  'satellite_altitude': 35786023.4375, 'platform_name': 'GOES-16',
                                   'calibration': 'reflectance', 'units': '%', 'wavelength': (0.45, 0.47, 0.49),
                                   'name': 'C01', 'resolution': 1000, 'sensor': 'abi',
                                   'start_time': '2017-09-20 17:30:40.800000', 'end_time': '2017-09-20 17:41:17.500000',
@@ -293,7 +297,7 @@ def test_reflectance_corrector_abi(self):
         self.assertIsInstance(res.data, da.Array)
         self.assertEqual(res.attrs['satellite_longitude'], -89.5)
         self.assertEqual(res.attrs['satellite_latitude'], 0.0)
-        self.assertEqual(res.attrs['satellite_altitude'], 35786.0234375)
+        self.assertEqual(res.attrs['satellite_altitude'], 35786023.4375)
         self.assertEqual(res.attrs['modifiers'], ('sunz_corrected', 'rayleigh_corrected_crefl',))
         self.assertEqual(res.attrs['platform_name'], 'GOES-16')
         self.assertEqual(res.attrs['calibration'], 'reflectance')
@@ -474,9 +478,50 @@ def make_xarray(self, name, calibration, wavelength=None, modifiers=None, resolu
         np.testing.assert_allclose(unique, [24.641586, 50.431692, 69.315375])
 
 
+class ViirsReflectanceCorrectorTest(unittest.TestCase):
+    def setUp(self):
+        """Patch in-class imports."""
+        self.astronomy = mock.MagicMock()
+        self.orbital = mock.MagicMock()
+        modules = {
+            'pyorbital.astronomy': self.astronomy,
+            'pyorbital.orbital': self.orbital,
+        }
+        self.module_patcher = mock.patch.dict('sys.modules', modules)
+        self.module_patcher.start()
+
+    def tearDown(self):
+        """Unpatch in-class imports."""
+        self.module_patcher.stop()
+
+    @mock.patch('satpy.composites.viirs.get_satpos')
+    def test_get_angles(self, get_satpos):
+        """Test sun and satellite angle calculation."""
+        from satpy.composites.viirs import ReflectanceCorrector
+
+        # Patch methods
+        get_satpos.return_value = 'sat_lon', 'sat_lat', 12345678
+        self.orbital.get_observer_look.return_value = 0, 0
+        self.astronomy.get_alt_az.return_value = 0, 0
+        area = mock.MagicMock()
+        area.get_lonlats_dask.return_value = 'lons', 'lats'
+        vis = mock.MagicMock(attrs={'area': area,
+                                    'start_time': 'start_time'})
+
+        # Compute angles
+        psp = ReflectanceCorrector(name='dummy')
+        psp.get_angles(vis)
+
+        # Check arguments of get_orbserver_look() call, especially the altitude
+        # unit conversion from meters to kilometers
+        self.orbital.get_observer_look.assert_called_with(
+            'sat_lon', 'sat_lat', 12345.678, 'start_time', 'lons', 'lats', 0)
+
+
 def suite():
     """Create test suite for test_ahi."""
     loader = unittest.TestLoader()
     mysuite = unittest.TestSuite()
     mysuite.addTest(loader.loadTestsFromTestCase(TestVIIRSComposites))
+    mysuite.addTest(loader.loadTestsFromTestCase(ViirsReflectanceCorrectorTest))
     return mysuite

satpy/tests/reader_tests/test_abi_l1b.py

@@ -154,15 +154,15 @@ def test_get_dataset(self):
                'orbital_parameters': {'projection_altitude': 1.0,
                                       'projection_latitude': 0.0,
                                       'projection_longitude': -90.0,
-                                      'satellite_nominal_altitude': 35786.02,
+                                      'satellite_nominal_altitude': 35786020.,
                                       'satellite_nominal_latitude': 0.0,
                                       'satellite_nominal_longitude': -89.5,
                                       'yaw_flip': True},
                'orbital_slot': None,
                'platform_name': 'GOES-16',
                'platform_shortname': 'G16',
                'production_site': None,
-               'satellite_altitude': 35786.02,
+               'satellite_altitude': 35786020.,
                'satellite_latitude': 0.0,
                'satellite_longitude': -89.5,
                'scan_mode': 'M3',

satpy/tests/reader_tests/test_ami_l1b.py

@@ -141,7 +141,7 @@ def _check_orbital_parameters(self, orb_params):
             'projection_altitude': 35785863.0,
             'projection_latitude': 0.0,
             'projection_longitude': 128.2,
-            'satellite_actual_altitude': 35782.65456070405,
+            'satellite_actual_altitude': 35782654.56070405,
             'satellite_actual_latitude': 0.005364927,
             'satellite_actual_longitude': 128.2707,
         }