Bugfix for Sentinel-2 radiance calculation

satpy/readers/msi_safe.py

@@ -15,7 +15,7 @@
 #
 # You should have received a copy of the GNU General Public License along with
 # satpy.  If not, see <http://www.gnu.org/licenses/>.
-"""SAFE MSI L1C reader.
+"""SAFE MSI L1C/L2A reader.
 
 The MSI data has a special value for saturated pixels. By default, these
 pixels are set to np.inf, but for some applications it might be desirable
@@ -32,6 +32,10 @@
 
   https://sentinels.copernicus.eu/documents/247904/685211/S2-PDGS-TAS-DI-PSD-V14.9.pdf/3d3b6c9c-4334-dcc4-3aa7-f7c0deffbaf7?t=1643013091529
 
+NOTE: At present, L1B data is not supported. If the user needs radiance data instead of counts or reflectances, these
+are retrieved by first calculating the reflectance and then working back to the radiance. L1B radiance data support
+will be added once the data is published onto the Copernicus data ecosystem.
+
 """
 
 import logging
@@ -59,13 +63,16 @@
 class SAFEMSIL1C(BaseFileHandler):
     """File handler for SAFE MSI files (jp2)."""
 
-    def __init__(self, filename, filename_info, filetype_info, mda, tile_mda, mask_saturated=True):
+    def __init__(self, filename, filename_info, filetype_info, mda, tile_mda,
+                 mask_saturated=True):
         """Initialize the reader."""
         super(SAFEMSIL1C, self).__init__(filename, filename_info,
                                          filetype_info)
         del mask_saturated
         self._channel = filename_info["band_name"]
         self.process_level = filename_info["process_level"]
+        if self.process_level not in ["L1C", "L2A"]:
+            raise ValueError(f"Unsupported process level: {self.process_level}")
         self._tile_mda = tile_mda
         self._mda = mda
         self.platform_name = PLATFORMS[filename_info["fmission_id"]]
@@ -83,7 +90,6 @@
         if proj is None:
             return
         proj.attrs = info.copy()
-        proj.attrs["units"] = "%"
         proj.attrs["platform_name"] = self.platform_name
         return proj
 
@@ -93,7 +99,21 @@
         if key["calibration"] == "reflectance":
             return self._mda.calibrate_to_reflectances(proj, self._channel)
         if key["calibration"] == "radiance":
-            return self._mda.calibrate_to_radiances(proj, self._channel)
+            # The calibration procedure differs for L1B and L1C/L2A data!
+            if self.process_level in ["L1C", "L2A"]:
+                # For higher level data, radiances must be computed from the reflectance.
+                # By default, we use the mean solar angles so that the user does not need to resample,
+                # but the user can also choose to use the solar angles from the tile metadata.
+                # This is on a coarse grid so for most bands must be resampled before use.
+                dq = dict(name="solar_zenith_angle", resolution=key["resolution"])
+                zen = self._tile_mda.get_dataset(dq, dict(xml_tag="Sun_Angles_Grid/Zenith"))
+                tmp_refl = self._mda.calibrate_to_reflectances(proj, self._channel)
+                return self._mda.calibrate_to_radiances(tmp_refl, zen, self._channel)
+            #else:
+                # For L1B the radiances can be directly computed from the digital counts.
+                #return self._mda.calibrate_to_radiances_l1b(proj, self._channel)
+
+
         if key["calibration"] == "counts":
             return self._mda._sanitize_data(proj)
         if key["calibration"] in ["aerosol_thickness", "water_vapor"]:
@@ -149,15 +169,15 @@
 
     def calibrate_to_reflectances(self, data, band_name):
         """Calibrate *data* using the radiometric information for the metadata."""
-        quantification = int(self.root.find(".//QUANTIFICATION_VALUE").text) if self.process_level == "L1C" else \
+        quantification = int(self.root.find(".//QUANTIFICATION_VALUE").text) if self.process_level[:2] == "L1" else \
             int(self.root.find(".//BOA_QUANTIFICATION_VALUE").text)
         data = self._sanitize_data(data)
         return (data + self.band_offset(band_name)) / quantification * 100
 
     def calibrate_to_atmospheric(self, data, band_name):
         """Calibrate L2A AOT/WVP product."""
         atmospheric_bands = ["AOT", "WVP"]
-        if self.process_level == "L1C":
+        if self.process_level == "L1C" or self.process_level == "L1B":
             return
         elif self.process_level == "L2A" and band_name not in atmospheric_bands:
             return
@@ -194,14 +214,37 @@
     @cached_property
     def band_offsets(self):
         """Get the band offsets from the metadata."""
-        offsets = self.root.find(".//Radiometric_Offset_List") if self.process_level == "L1C" else \
+        offsets = self.root.find(".//Radiometric_Offset_List") if self.process_level[:2] == "L1" else \
             self.root.find(".//BOA_ADD_OFFSET_VALUES_LIST")
         if offsets is not None:
             band_offsets = {int(off.attrib["band_id"]): float(off.text) for off in offsets}
         else:
             band_offsets = {}
         return band_offsets
 
+    def solar_irradiance(self, band_name):
+        """Get the solar irradiance for a given *band_name*."""
+        band_index = self._band_index(band_name)
+        return self.solar_irradiances[band_index]
+
+    @cached_property
+    def solar_irradiances(self):
+        """Get the TOA solar irradiance values from the metadata."""
+        irrads = self.root.find(".//Solar_Irradiance_List")
+        if irrads is not None:
+            solar_irrad = {int(irr.attrib["bandId"]): float(irr.text) for irr in irrads}
+        else:
+            solar_irrad = {}
+        return solar_irrad
+
+    @cached_property
+    def sun_earth_dist(self):
+        """Get the sun-earth distance from the metadata."""
+        sed = self.root.find(".//U")
+        if sed is not None:
+            return float(sed.text)
+        return -1
+
     @cached_property
     def special_values(self):
         """Get the special values from the metadata."""
@@ -219,12 +262,23 @@
         """Get the saturated value from the metadata."""
         return self.special_values["SATURATED"]
 
-    def calibrate_to_radiances(self, data, band_name):
+    def calibrate_to_radiances_l1b(self, data, band_name):
         """Calibrate *data* to radiance using the radiometric information for the metadata."""
         physical_gain = self.physical_gain(band_name)
         data = self._sanitize_data(data)
         return (data + self.band_offset(band_name)) / physical_gain
 
+    def calibrate_to_radiances(self, data, solar_zenith, band_name):
+        """Calibrate *data* to radiance using the radiometric information for the metadata."""
+        sed = self.sun_earth_dist
+        if sed < 0.5 or sed > 1.5:
+            raise ValueError(f"Sun-Earth distance is incorrect in the metadata: {sed}")
+        solar_irrad_band = self.solar_irradiance(band_name)
+
+        solar_zenith = np.deg2rad(solar_zenith)
+
+        return (data / 100.) * solar_irrad_band * np.cos(solar_zenith) / (np.pi * sed * sed)
+
     def physical_gain(self, band_name):
         """Get the physical gain for a given *band_name*."""
         band_index = self._band_index(band_name)

satpy/tests/reader_tests/test_msi_safe.py

@@ -1435,18 +1435,25 @@
     return xml_fh, xml_tile_fh
 
 
-def jp2_builder(process_level, band_name, mask_saturated=True):
+def jp2_builder(process_level, band_name, mask_saturated=True, test_l1b=False):
     """Build fake SAFE jp2 image file."""
     from satpy.readers.msi_safe import SAFEMSIL1C, SAFEMSITileMDXML
     filename_info = dict(observation_time=fname_dt, dtile_number=None, band_name=band_name, fmission_id="S2A",
                          process_level=process_level.replace("old", ""))
+    if test_l1b:
+        filename_info["process_level"] = "L1B"
+
     xml_fh = xml_builder(process_level, mask_saturated, band_name)[0]
     tile_xml_fh = mock.create_autospec(SAFEMSITileMDXML)(BytesIO(TILE_XMLS[PROCESS_LEVELS.index(process_level)]),
-                                   filename_info, mock.MagicMock())
+                                                         filename_info, mock.MagicMock())
     tile_xml_fh.start_time.return_value = tilemd_dt
+    tile_xml_fh.get_dataset.return_value = xr.DataArray([[22.5, 23.8],
+                                                          [22.5, 24.8]],
+                                                        dims=["x", "y"])
     jp2_fh = SAFEMSIL1C("somefile", filename_info, mock.MagicMock(), xml_fh, tile_xml_fh)
     return jp2_fh
 
+
 def make_alt_dataid(**items):
     """Make a DataID with modified keys."""
     from satpy.dataset.dataid import DataID, ModifierTuple, WavelengthRange
@@ -1578,26 +1585,26 @@
                              [
                                  ("L1C", True, "B01", ([[[np.nan, -9.99, -9.98, -9.97],
                                                          [-9.96, 0, 645.34, np.inf]]],
-                                                       [[[np.nan, -251.584265, -251.332429, -251.080593],
-                                                         [-250.828757, 0., 16251.99095, np.inf]]],
+                                                       [[[0.0, 5.60879825, 11.2175965, 16.8263948,],
+                                                         [22.435193, 5608.79825, 367566.985, 367572.593]]],
                                                        [[[np.nan, 1, 2, 3],
                                                          [4, 1000, 65534, np.inf]]])),
                                  ("L1C", False, "B10", ([[[np.nan, -19.99, -19.98, -19.97],
                                                           [-19.96, -10, 635.34, 635.35]]],
-                                                        [[[np.nan, -35.465976, -35.448234, -35.430493],
-                                                          [-35.412751, -17.741859, 1127.211275, 1127.229017]]],
+                                                        [[[0.0, 1.09348075, 2.1869615, 3.28044225],
+                                                          [4.373923, 1093.48075, 71660.1675, 71661.2609]]],
                                                         [[[np.nan, 1, 2, 3],
                                                           [4, 1000, 65534, 65535]]])),
                                  ("oldL1C", True, "B01", ([[[np.nan, 0.01, 0.02, 0.03],
                                                             [0.04, 10, 655.34, np.inf]]],
-                                                          [[[np.nan, 0.251836101, 0.503672202, 0.755508303],
-                                                            [1.00734440, 251.836101, 16503.8271, np.inf]]],
+                                                          [[[0.0, 5.60879825, 11.2175965, 16.8263948,],
+                                                            [22.435193, 5608.79825, 367566.985, 367572.593]]],
                                                           [[[np.nan, 1, 2, 3],
                                                             [4, 1000, 65534, np.inf]]])),
                                  ("L2A", False, "B03", ([[[np.nan, -9.99, -9.98, -9.97],
                                                           [-9.96, 0, 645.34, 645.35]]],
-                                                        [[[np.nan, -238.571863, -238.333052, -238.094241],
-                                                          [-237.855431, 0, 15411.407995, 15411.646806]]],
+                                                        [[[0.0, 5.25188783, 10.5037757, 15.7556635,],
+                                                          [21.0075513, 5251.88783, 344177.217, 344182.469]]],
                                                         [[[np.nan, 1, 2, 3],
                                                           [4, 1000, 65534, 65535]]])),
                              ])
@@ -1606,10 +1613,11 @@
         xml_fh = xml_builder(process_level, mask_saturated)[0]
 
         res1 = xml_fh.calibrate_to_reflectances(self.fake_data, band_name)
-        res2 = xml_fh.calibrate_to_radiances(self.fake_data, band_name)
+        res2 = xml_fh.calibrate_to_radiances(self.fake_data, 25.6, band_name)
         res3 = xml_fh._sanitize_data(self.fake_data)
 
         results = (res1, res2, res3)
+
         np.testing.assert_allclose(results, expected)
 
     @pytest.mark.parametrize(("process_level", "mask_saturated", "band_name", "expected"),
@@ -1640,22 +1648,24 @@
     def setup_method(self):
         """Set up the test."""
         self.fake_data = xr.Dataset({"band_data": xr.DataArray([[[0, 1], [65534, 65535]]], dims=["band", "x", "y"])})
+        self.fake_data = xr.Dataset({"band_data": xr.DataArray([[[0, 1], [65534, 65535]]], dims=["band", "x", "y"])})
 
-    @pytest.mark.parametrize(("mask_saturated", "dataset_name", "calibration", "expected"),
+    @pytest.mark.parametrize(("process_level", "mask_saturated", "dataset_name", "calibration", "expected"),
                              [
-                                 (False, "B01", "reflectance", [[np.nan, -9.99], [645.34, 645.35]]),
-                                 (True, "B02", "radiance", [[np.nan, -265.970568], [17181.325973, np.inf]]),
-                                 (True, "B03", "counts", [[np.nan, 1], [65534, np.inf]]),
-                                 (False, "AOT", "aerosol_thickness", [[np.nan, 0.001], [65.534, 65.535]]),
-                                 (True, "WVP", "water_vapor", [[np.nan, 0.001], [65.534, np.inf]]),
-                                 (True, "SNOW", "water_vapor", None),
+                                 ("L2A", False, "B01", "reflectance", [[np.nan, -9.99], [645.34, 645.35]]),
+                                 ("L1C", True, "B02", "radiance", [[np.nan, -59.439197], [3877.121602, np.inf]]),
+                                 ("L2A", True, "B03", "counts", [[np.nan, 1], [65534, np.inf]]),
+                                 ("L2A", False, "AOT", "aerosol_thickness", [[np.nan, 0.001], [65.534, 65.535]]),
+                                 ("L2A", True, "WVP", "water_vapor", [[np.nan, 0.001], [65.534, np.inf]]),
+                                 ("L2A", True, "SNOW", "water_vapor", None),
                              ])
-    def test_calibration_and_masking(self, mask_saturated, dataset_name, calibration, expected):
+    def test_calibration_and_masking(self, process_level, mask_saturated, dataset_name, calibration, expected):
         """Test that saturated is masked with inf when requested and that calibration is performed."""
-        jp2_fh = jp2_builder("L2A", dataset_name, mask_saturated)
+        jp2_fh = jp2_builder(process_level, dataset_name, mask_saturated)
 
         with mock.patch("xarray.open_dataset", return_value=self.fake_data):
-            res = jp2_fh.get_dataset(make_alt_dataid(name=dataset_name, calibration=calibration), info=dict())
+            res = jp2_fh.get_dataset(make_alt_dataid(name=dataset_name, calibration=calibration, resolution="20"),
+                                     info=dict())
             if res is not None:
                 np.testing.assert_allclose(res, expected)
             else:
@@ -1677,7 +1687,13 @@
             assert res1 is None
             assert res2 is None
 
-    def test_start_time(self):
+    def test_start_end_time(self):
         """Test that the correct start time is returned."""
         jp2_fh = jp2_builder("L1C", "B01")
         assert tilemd_dt == jp2_fh.start_time
+        assert tilemd_dt == jp2_fh.end_time
+
+    def test_l1b_error(self):
+        """We can't process L1B data yet, so check an error is raised."""
+        with pytest.raises(ValueError, match="Unsupported process level: L1B"):
+            jp2_builder("L1C", "B01", test_l1b=True)