From 836af09c0be569db7f3f399277d666dbb7541505 Mon Sep 17 00:00:00 2001 From: andream Date: Mon, 6 Apr 2020 17:23:09 +0000 Subject: [PATCH 01/14] implement fixes to calibration of visnir and ir channels --- satpy/readers/fci_l1c_fdhsi.py | 33 +++++++++++++++------------------ 1 file changed, 15 insertions(+), 18 deletions(-) diff --git a/satpy/readers/fci_l1c_fdhsi.py b/satpy/readers/fci_l1c_fdhsi.py index 448a37222b..87d5dbff9a 100644 --- a/satpy/readers/fci_l1c_fdhsi.py +++ b/satpy/readers/fci_l1c_fdhsi.py @@ -257,8 +257,8 @@ def calibrate(self, data, key, measured, root): def _ir_calibrate(self, radiance, measured, root): """IR channel calibration.""" - coef = self[measured + "/radiance_unit_conversion_coefficient"] - wl_c = self[root + "/central_wavelength_actual"] + # using the method from EUM/MET/TEN/11/0569 and PUG + vc = self[measured + "/radiance_to_bt_conversion_coefficient_wavenumber"] a = self[measured + "/radiance_to_bt_conversion_coefficient_a"] b = self[measured + "/radiance_to_bt_conversion_coefficient_b"] @@ -266,7 +266,7 @@ def _ir_calibrate(self, radiance, measured, root): c1 = self[measured + "/radiance_to_bt_conversion_constant_c1"] c2 = self[measured + "/radiance_to_bt_conversion_constant_c2"] - for v in (coef, wl_c, a, b, c1, c2): + for v in (vc, a, b, c1, c2): if v == v.attrs.get("FillValue", default_fillvals.get(v.dtype.str[1:])): logger.error( @@ -283,10 +283,8 @@ def _ir_calibrate(self, radiance, measured, root): coords=radiance.coords, attrs=radiance.attrs) - Lv = radiance * coef - vc = 1e6/wl_c # from wl in um to wn in m^-1 nom = c2 * vc - denom = a * np.log(1 + (c1 * vc**3) / Lv) + denom = a * np.log(1 + (c1 * vc**3) / radiance) res = nom / denom - b / a res.attrs["units"] = "K" @@ -294,25 +292,24 @@ def _ir_calibrate(self, radiance, measured, root): def _vis_calibrate(self, radiance, measured): """VIS channel calibration.""" - # radiance to reflectance taken as in mipp/xrit/MSG.py - # again FCI User Guide is not clear on how to do this - cesilab = measured + "/channel_effective_solar_irradiance" cesi = self[cesilab] + if cesi == cesi.attrs.get( "FillValue", default_fillvals.get(cesi.dtype.str[1:])): logger.error( "channel effective solar irradiance set to fill value, " "cannot produce reflectance for {:s}.".format(measured)) return xr.DataArray( - da.full(shape=radiance.shape, - chunks=radiance.chunks, - fill_value=np.nan), - dims=radiance.dims, - coords=radiance.coords, - attrs=radiance.attrs) - - sirr = float(cesi) - res = radiance / sirr * 100 + da.full(shape=radiance.shape, + chunks=radiance.chunks, + fill_value=np.nan), + dims=radiance.dims, + coords=radiance.coords, + attrs=radiance.attrs) + + sun_earth_distance = np.mean(self["state/celestial/earth_sun_distance"]) / 149597870.7 # [AU] + + res = 100 * radiance * np.pi * sun_earth_distance**2 / cesi res.attrs["units"] = "%" return res From b22f4bf2e5ed6a48dc9c7dc8be60412f9594ffac Mon Sep 17 00:00:00 2001 From: andream Date: Wed, 8 Apr 2020 09:21:53 +0000 Subject: [PATCH 02/14] fix tests for reflectance and bt calculation --- satpy/readers/fci_l1c_fdhsi.py | 6 +++--- .../tests/reader_tests/test_fci_l1c_fdhsi.py | 19 +++++++++---------- 2 files changed, 12 insertions(+), 13 deletions(-) diff --git a/satpy/readers/fci_l1c_fdhsi.py b/satpy/readers/fci_l1c_fdhsi.py index 87d5dbff9a..72ca90f68a 100644 --- a/satpy/readers/fci_l1c_fdhsi.py +++ b/satpy/readers/fci_l1c_fdhsi.py @@ -200,7 +200,7 @@ def calc_area_extent(self, key): ext[c] = (min_c.item(), max_c.item()) area_extent = (ext["x"][1], ext["y"][1], ext["x"][0], ext["y"][0]) - return (area_extent, nlines, ncols) + return area_extent, nlines, ncols def get_area_def(self, key, info=None): """Calculate on-fly area definition for 0 degree geos-projection for a dataset.""" @@ -216,7 +216,7 @@ def get_area_def(self, key, info=None): lon_0 = float(self["data/mtg_geos_projection/attr/longitude_of_projection_origin"]) sweep = str(self["data/mtg_geos_projection"].sweep_angle_axis) # Channel dependent swath resolution - (area_extent, nlines, ncols) = self.calc_area_extent(key) + area_extent, nlines, ncols = self.calc_area_extent(key) logger.debug('Calculated area extent: {}' .format(''.join(str(area_extent)))) @@ -308,7 +308,7 @@ def _vis_calibrate(self, radiance, measured): coords=radiance.coords, attrs=radiance.attrs) - sun_earth_distance = np.mean(self["state/celestial/earth_sun_distance"]) / 149597870.7 # [AU] + sun_earth_distance = np.mean(self["state/celestial/earth_sun_distance"]) / 149597870.7 # [AU] res = 100 * radiance * np.pi * sun_earth_distance**2 / cesi res.attrs["units"] = "%" diff --git a/satpy/tests/reader_tests/test_fci_l1c_fdhsi.py b/satpy/tests/reader_tests/test_fci_l1c_fdhsi.py index 7753fbe959..fc5453a8d5 100644 --- a/satpy/tests/reader_tests/test_fci_l1c_fdhsi.py +++ b/satpy/tests/reader_tests/test_fci_l1c_fdhsi.py @@ -38,17 +38,17 @@ def _get_test_calib_for_channel_ir(self, chroot, meas): C_SPEED as c) xrda = xr.DataArray data = {} - data[meas + "/radiance_unit_conversion_coefficient"] = xrda(1) - data[chroot + "/central_wavelength_actual"] = xrda(10) - data[meas + "/radiance_to_bt_conversion_coefficient_a"] = xrda(1000) - data[meas + "/radiance_to_bt_conversion_coefficient_b"] = xrda(1) - data[meas + "/radiance_to_bt_conversion_constant_c1"] = xrda(2*h*c**2) - data[meas + "/radiance_to_bt_conversion_constant_c2"] = xrda(h*c/k) + data[meas + "/radiance_to_bt_conversion_coefficient_wavenumber"] = xrda(955) + data[meas + "/radiance_to_bt_conversion_coefficient_a"] = xrda(1) + data[meas + "/radiance_to_bt_conversion_coefficient_b"] = xrda(0.4) + data[meas + "/radiance_to_bt_conversion_constant_c1"] = xrda(1e11*2*h*c**2) + data[meas + "/radiance_to_bt_conversion_constant_c2"] = xrda(1e2*h*c/k) return data def _get_test_calib_for_channel_vis(self, chroot, meas): xrda = xr.DataArray data = {} + data["state/celestial/earth_sun_distance"] = xrda(149597870.7) data[meas + "/channel_effective_solar_irradiance"] = xrda(50) return data @@ -173,8 +173,7 @@ def _get_test_calib_for_channel_ir(self, chroot, meas): from netCDF4 import default_fillvals v = xr.DataArray(default_fillvals["f4"]) data = {} - data[meas + "/radiance_unit_conversion_coefficient"] = v - data[chroot + "/central_wavelength_actual"] = v + data[meas + "/radiance_to_bt_conversion_coefficient_wavenumber"] = v data[meas + "/radiance_to_bt_conversion_coefficient_a"] = v data[meas + "/radiance_to_bt_conversion_coefficient_b"] = v data[meas + "/radiance_to_bt_conversion_constant_c1"] = v @@ -337,7 +336,7 @@ def test_load_reflectance(self): self.assertEqual(res[ch].dtype, np.float64) self.assertEqual(res[ch].attrs["calibration"], "reflectance") self.assertEqual(res[ch].attrs["units"], "%") - numpy.testing.assert_array_equal(res[ch], 15 / 50 * 100) + numpy.testing.assert_array_equal(res[ch], 100 * 15 * 1 * np.pi / 50) def test_load_bt(self): """Test loading with bt @@ -366,7 +365,7 @@ def test_load_bt(self): self.assertEqual(res[ch].attrs["units"], "K") numpy.testing.assert_array_almost_equal( res[ch], - 181.917084) + 209.68274099) def test_load_composite(self): """Test that composites are loadable From af8d4488c321c3fb23b41b7ee6d0ea103113baae Mon Sep 17 00:00:00 2001 From: andream Date: Wed, 8 Apr 2020 11:57:19 +0000 Subject: [PATCH 03/14] implement warm scaling for ir_38 --- satpy/readers/fci_l1c_fdhsi.py | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/satpy/readers/fci_l1c_fdhsi.py b/satpy/readers/fci_l1c_fdhsi.py index 72ca90f68a..0f5475142e 100644 --- a/satpy/readers/fci_l1c_fdhsi.py +++ b/satpy/readers/fci_l1c_fdhsi.py @@ -117,6 +117,7 @@ def get_dataset(self, key, info=None): attrs = data.attrs.copy() info = info.copy() + fv = attrs.pop( "FillValue", default_fillvals.get(data.dtype.str[1:], np.nan)) @@ -128,6 +129,7 @@ def get_dataset(self, key, info=None): nfv = np.nan data = data.where(data >= vr[0], nfv) data = data.where(data <= vr[1], nfv) + if key.calibration == "counts": # from package description, this just means not applying add_offset # and scale_factor @@ -136,8 +138,17 @@ def get_dataset(self, key, info=None): data.attrs["units"] = "1" res = data else: - data = (data * attrs.pop("scale_factor", 1) + - attrs.pop("add_offset", 0)) + # counts to radiance scaling + if key.name == 'ir_38': + data = xr.where(((2**12-1 < data) & (data <= 2**13-1)), + (data * attrs.pop("warm_scale_factor", 1) + + attrs.pop("warm_add_offset", 0)), + (data * attrs.pop("scale_factor", 1) + + attrs.pop("add_offset", 0)) + ) + else: + data = (data * attrs.pop("scale_factor", 1) + + attrs.pop("add_offset", 0)) if key.calibration in ("brightness_temperature", "reflectance"): res = self.calibrate(data, key, measured, root) From 381411bb641228cbba21052e152f656b950ac6e6 Mon Sep 17 00:00:00 2001 From: andream Date: Wed, 8 Apr 2020 15:27:24 +0000 Subject: [PATCH 04/14] add testing of warm range and fix flake8 --- .../tests/reader_tests/test_fci_l1c_fdhsi.py | 93 ++++++++++++------- 1 file changed, 59 insertions(+), 34 deletions(-) diff --git a/satpy/tests/reader_tests/test_fci_l1c_fdhsi.py b/satpy/tests/reader_tests/test_fci_l1c_fdhsi.py index fc5453a8d5..d649cae780 100644 --- a/satpy/tests/reader_tests/test_fci_l1c_fdhsi.py +++ b/satpy/tests/reader_tests/test_fci_l1c_fdhsi.py @@ -31,6 +31,8 @@ class FakeNetCDF4FileHandler2(FakeNetCDF4FileHandler): + """Class for faking the NetCDF4 Filehandler.""" + def _get_test_calib_for_channel_ir(self, chroot, meas): from pyspectral.blackbody import ( H_PLANCK as h, @@ -66,7 +68,24 @@ def _get_test_content_for_channel(self, pat, ch): data = {} ch_str = pat.format(ch) ch_path = rad.format(ch_str) - d = xrda( + + if ch == 38: + fire_line = da.ones((1, ncols), dtype="uint16", chunks=1024) * 5000 + data_without_fires = da.ones((nrows-1, ncols), dtype="uint16", chunks=1024) + d = xrda( + da.concatenate([fire_line, data_without_fires], axis=0), + dims=("y", "x"), + attrs={ + "valid_range": [0, 8191], + "scale_factor": 5, + "add_offset": 10, + "warm_scale_factor": 2, + "warm_add_offset": -300, + "units": "mW.m-2.sr-1.(cm-1)-1", + } + ) + else: + d = xrda( da.ones((nrows, ncols), dtype="uint16", chunks=1024), dims=("y", "x"), attrs={ @@ -74,8 +93,9 @@ def _get_test_content_for_channel(self, pat, ch): "scale_factor": 5, "add_offset": 10, "units": "mW.m-2.sr-1.(cm-1)-1", - } - ) + } + ) + data[ch_path] = d data[x.format(ch_str)] = xrda( da.arange(1, ncols+1, dtype="uint16"), @@ -152,6 +172,7 @@ def _get_test_content_areadef(self): return data def get_test_content(self, filename, filename_info, filetype_info): + """Get the content of the test data.""" # mock global attributes # - root groups global # - other groups global @@ -167,8 +188,8 @@ def get_test_content(self, filename, filename_info, filetype_info): class FakeNetCDF4FileHandler3(FakeNetCDF4FileHandler2): - """Mock bad data - """ + """Mock bad data.""" + def _get_test_calib_for_channel_ir(self, chroot, meas): from netCDF4 import default_fillvals v = xr.DataArray(default_fillvals["f4"]) @@ -182,14 +203,14 @@ def _get_test_calib_for_channel_ir(self, chroot, meas): class TestFCIL1CFDHSIReader(unittest.TestCase): + """Initialize the unittest TestCase for the FCI L1C FDHSI Reader.""" + yaml_file = "fci_l1c_fdhsi.yaml" _alt_handler = FakeNetCDF4FileHandler2 def setUp(self): - """Wrap NetCDF4 FileHandler with our own fake handler - """ - + """Wrap NetCDF4 FileHandler with our own fake handler.""" # implementation strongly inspired by test_viirs_l1b.py from satpy.config import config_search_paths from satpy.readers.fci_l1c_fdhsi import FCIFDHSIFileHandler @@ -204,25 +225,21 @@ def setUp(self): self.p.is_local = True def tearDown(self): - """Stop wrapping the NetCDF4 file handler - """ + """Stop wrapping the NetCDF4 file handler.""" # implementation strongly inspired by test_viirs_l1b.py self.p.stop() class TestFCIL1CFDHSIReaderGoodData(TestFCIL1CFDHSIReader): - """Test FCI L1C FDHSI reader - """ + """Test FCI L1C FDHSI reader.""" # TODO: - # - test special case for extended range IR38 # - test geolocation _alt_handler = FakeNetCDF4FileHandler2 def test_file_pattern(self): - """Test file pattern matching - """ + """Test file pattern matching.""" from satpy.readers import load_reader filenames = [ @@ -257,8 +274,7 @@ def test_file_pattern(self): "ir_123", "ir_133"]} def test_load_counts(self): - """Test loading with counts - """ + """Test loading with counts.""" from satpy import DatasetID from satpy.readers import load_reader @@ -284,11 +300,15 @@ def test_load_counts(self): self.assertEqual(res[ch].dtype, np.uint16) self.assertEqual(res[ch].attrs["calibration"], "counts") self.assertEqual(res[ch].attrs["units"], "1") - numpy.testing.assert_array_equal(res[ch], 1) + + if ch == 'ir_38': + numpy.testing.assert_array_equal(res[ch][~0], 1) + numpy.testing.assert_array_equal(res[ch][0], 5000) + else: + numpy.testing.assert_array_equal(res[ch], 1) def test_load_radiance(self): - """Test loading with radiance - """ + """Test loading with radiance.""" from satpy import DatasetID from satpy.readers import load_reader @@ -310,11 +330,15 @@ def test_load_radiance(self): self.assertEqual(res[ch].dtype, np.float64) self.assertEqual(res[ch].attrs["calibration"], "radiance") self.assertEqual(res[ch].attrs["units"], 'mW.m-2.sr-1.(cm-1)-1') - numpy.testing.assert_array_equal(res[ch], 15) + + if ch == 'ir_38': + numpy.testing.assert_array_equal(res[ch][~0], 15) + numpy.testing.assert_array_equal(res[ch][0], 9700) + else: + numpy.testing.assert_array_equal(res[ch], 15) def test_load_reflectance(self): - """Test loading with reflectance - """ + """Test loading with reflectance.""" from satpy import DatasetID from satpy.readers import load_reader @@ -339,8 +363,7 @@ def test_load_reflectance(self): numpy.testing.assert_array_equal(res[ch], 100 * 15 * 1 * np.pi / 50) def test_load_bt(self): - """Test loading with bt - """ + """Test loading with bt.""" from satpy import DatasetID from satpy.readers import load_reader @@ -363,14 +386,15 @@ def test_load_bt(self): self.assertEqual(res[ch].attrs["calibration"], "brightness_temperature") self.assertEqual(res[ch].attrs["units"], "K") - numpy.testing.assert_array_almost_equal( - res[ch], - 209.68274099) - def test_load_composite(self): - """Test that composites are loadable - """ + if ch == 'ir_38': + numpy.testing.assert_array_almost_equal(res[ch][~0], 209.68274099) + numpy.testing.assert_array_almost_equal(res[ch][0], 1888.851296) + else: + numpy.testing.assert_array_almost_equal(res[ch], 209.68274099) + def test_load_composite(self): + """Test that composites are loadable.""" # when dedicated composites for FCI FDHSI are implemented in satpy, # this method should probably move to a dedicated class and module # in the tests.compositor_tests package @@ -383,11 +407,12 @@ def test_load_composite(self): class TestFCIL1CFDHSIReaderBadData(TestFCIL1CFDHSIReader): + """Test the FCI L1C FDHSI Reader for bad data input.""" + _alt_handler = FakeNetCDF4FileHandler3 def test_handling_bad_data_ir(self): - """Test handling of bad data - """ + """Test handling of bad data.""" from satpy import DatasetID from satpy.readers import load_reader @@ -406,4 +431,4 @@ def test_handling_bad_data_ir(self): reader.load([DatasetID( name="ir_123", calibration="brightness_temperature")]) - self.assertRegex(cm.output[0], "cannot produce brightness temperatur") + self.assertRegex(cm.output[0], "cannot produce brightness temperature") From 608dbe869cf6434d9e4ce83bdada7fd1373831f6 Mon Sep 17 00:00:00 2001 From: andream Date: Wed, 8 Apr 2020 16:32:41 +0000 Subject: [PATCH 05/14] add link to technical note and fix PUG link --- satpy/readers/fci_l1c_fdhsi.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/satpy/readers/fci_l1c_fdhsi.py b/satpy/readers/fci_l1c_fdhsi.py index 0f5475142e..548dd0c045 100644 --- a/satpy/readers/fci_l1c_fdhsi.py +++ b/satpy/readers/fci_l1c_fdhsi.py @@ -37,7 +37,7 @@ radius * From the attribute ``perspective_point_height`` on the same data variable, the geostationary altitude in the normalised geostationary - projection (see PUG §5.2) + projection (see `PUG`_ §5.2) * From the attribute ``longitude_of_projection_origin`` on the same data variable, the longitude of the projection origin * From the attribute ``inverse_flattening`` on the same data variable, the @@ -53,6 +53,10 @@ ``pyresample.geometry.AreaDefinition``, which then uses proj4 for the actual geolocation calculations. +The brightness temperature calculation is based on the formulas indicated in +`PUG`_ and `RADTOBR`_. + +.. _RADTOBR: https://www.eumetsat.int/website/wcm/idc/idcplg?IdcService=GET_FILE&dDocName=PDF_EFFECT_RAD_TO_BRIGHTNESS&RevisionSelectionMethod=LatestReleased&Rendition=Web .. _PUG: http://www.eumetsat.int/website/wcm/idc/idcplg?IdcService=GET_FILE&dDocName=PDF_DMT_719113&RevisionSelectionMethod=LatestReleased&Rendition=Web .. _EUMETSAT: https://www.eumetsat.int/website/home/Satellites/FutureSatellites/MeteosatThirdGeneration/MTGDesign/index.html#fci # noqa: E501 """ @@ -268,7 +272,7 @@ def calibrate(self, data, key, measured, root): def _ir_calibrate(self, radiance, measured, root): """IR channel calibration.""" - # using the method from EUM/MET/TEN/11/0569 and PUG + # using the method from RADTOBR and PUG vc = self[measured + "/radiance_to_bt_conversion_coefficient_wavenumber"] a = self[measured + "/radiance_to_bt_conversion_coefficient_a"] From 41830316a3f07bcd58527dcbaf9a832550c787ae Mon Sep 17 00:00:00 2001 From: andream Date: Mon, 4 May 2020 14:54:32 +0000 Subject: [PATCH 06/14] delete intermediate cesilab variable --- satpy/readers/fci_l1c_fdhsi.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/satpy/readers/fci_l1c_fdhsi.py b/satpy/readers/fci_l1c_fdhsi.py index 548dd0c045..0ac4c3868a 100644 --- a/satpy/readers/fci_l1c_fdhsi.py +++ b/satpy/readers/fci_l1c_fdhsi.py @@ -307,8 +307,7 @@ def _ir_calibrate(self, radiance, measured, root): def _vis_calibrate(self, radiance, measured): """VIS channel calibration.""" - cesilab = measured + "/channel_effective_solar_irradiance" - cesi = self[cesilab] + cesi = self[measured + "/channel_effective_solar_irradiance"] if cesi == cesi.attrs.get( "FillValue", default_fillvals.get(cesi.dtype.str[1:])): From 59cd75cae2ab2fa998c7fcbaf962151f4b8528d1 Mon Sep 17 00:00:00 2001 From: andream Date: Mon, 4 May 2020 17:24:27 +0000 Subject: [PATCH 07/14] rearrange calibration functions and fix popping of attributes --- satpy/readers/fci_l1c_fdhsi.py | 73 +++++++++++++++++++--------------- 1 file changed, 40 insertions(+), 33 deletions(-) diff --git a/satpy/readers/fci_l1c_fdhsi.py b/satpy/readers/fci_l1c_fdhsi.py index 0ac4c3868a..45238ac2fc 100644 --- a/satpy/readers/fci_l1c_fdhsi.py +++ b/satpy/readers/fci_l1c_fdhsi.py @@ -123,9 +123,9 @@ def get_dataset(self, key, info=None): info = info.copy() fv = attrs.pop( - "FillValue", - default_fillvals.get(data.dtype.str[1:], np.nan)) - vr = attrs.pop("valid_range", [-np.inf, np.inf]) + "FillValue", + default_fillvals.get(data.dtype.str[1:], np.nan)) + vr = attrs.get("valid_range", [-np.inf, np.inf]) if key.calibration == "counts": attrs["_FillValue"] = fv nfv = fv @@ -134,31 +134,8 @@ def get_dataset(self, key, info=None): data = data.where(data >= vr[0], nfv) data = data.where(data <= vr[1], nfv) - if key.calibration == "counts": - # from package description, this just means not applying add_offset - # and scale_factor - attrs.pop("scale_factor") - attrs.pop("add_offset") - data.attrs["units"] = "1" - res = data - else: - # counts to radiance scaling - if key.name == 'ir_38': - data = xr.where(((2**12-1 < data) & (data <= 2**13-1)), - (data * attrs.pop("warm_scale_factor", 1) + - attrs.pop("warm_add_offset", 0)), - (data * attrs.pop("scale_factor", 1) + - attrs.pop("add_offset", 0)) - ) - else: - data = (data * attrs.pop("scale_factor", 1) + - attrs.pop("add_offset", 0)) + res = self.calibrate(data, key, measured, root) - if key.calibration in ("brightness_temperature", "reflectance"): - res = self.calibrate(data, key, measured, root) - else: - res = data - data.attrs["units"] = attrs["units"] # pre-calibration units no longer apply info.pop("units") attrs.pop("units") @@ -166,6 +143,17 @@ def get_dataset(self, key, info=None): res.attrs.update(key.to_dict()) res.attrs.update(info) res.attrs.update(attrs) + + # remove unpacking parameters for calibrated data + if key.calibration in ['brightness_temperature', 'reflectance']: + res.attrs.pop("add_offset") + res.attrs.pop("warm_add_offset") + res.attrs.pop("scale_factor") + res.attrs.pop("warm_scale_factor") + + # remove attributes from original file which don't apply anymore + res.attrs.pop('long_name') + return res def get_channel_dataset(self, channel): @@ -258,15 +246,34 @@ def get_area_def(self, key, info=None): def calibrate(self, data, key, measured, root): """Calibrate data.""" - if key.calibration == 'brightness_temperature': - data = self._ir_calibrate(data, measured, root) - elif key.calibration == 'reflectance': - data = self._vis_calibrate(data, measured) + if key.calibration == "counts": + # from package description, this just means not applying add_offset + # and scale_factor + data.attrs["units"] = "1" else: - raise ValueError( + original_radiance_units = data.attrs.get("units") + # counts to radiance scaling + if key.name == 'ir_38': + data = xr.where(((2**12-1 < data) & (data <= 2**13-1)), + (data * data.attrs.get("warm_scale_factor", 1) + + data.attrs.get("warm_add_offset", 0)), + (data * data.attrs.get("scale_factor", 1) + + data.attrs.get("add_offset", 0)) + ) + else: + data = (data * data.attrs.get("scale_factor", 1) + + data.attrs.get("add_offset", 0)) + + if key.calibration == 'brightness_temperature': + data = self._ir_calibrate(data, measured, root) + elif key.calibration == 'reflectance': + data = self._vis_calibrate(data, measured) + else: + logger.warning( "Received unknown calibration key. Expected " "'brightness_temperature' or 'reflectance', got " - + key.calibration) + + key.calibration + ". Returning radiances.") + data.attrs["units"] = original_radiance_units return data From 997eba7b373f7d2a9e58b7fde9f4e67fa96656d8 Mon Sep 17 00:00:00 2001 From: andream Date: Mon, 4 May 2020 17:41:43 +0000 Subject: [PATCH 08/14] make radiance calibration explicit, update tests --- satpy/readers/fci_l1c_fdhsi.py | 9 +++++---- satpy/tests/reader_tests/test_fci_l1c_fdhsi.py | 4 ++++ 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/satpy/readers/fci_l1c_fdhsi.py b/satpy/readers/fci_l1c_fdhsi.py index 45238ac2fc..a3cb4c7ba2 100644 --- a/satpy/readers/fci_l1c_fdhsi.py +++ b/satpy/readers/fci_l1c_fdhsi.py @@ -268,12 +268,13 @@ def calibrate(self, data, key, measured, root): data = self._ir_calibrate(data, measured, root) elif key.calibration == 'reflectance': data = self._vis_calibrate(data, measured) + elif key.calibration == 'radiance': + data.attrs["units"] = original_radiance_units else: - logger.warning( + logger.error( "Received unknown calibration key. Expected " - "'brightness_temperature' or 'reflectance', got " - + key.calibration + ". Returning radiances.") - data.attrs["units"] = original_radiance_units + "'brightness_temperature', 'reflectance' or 'radiance, got " + + key.calibration + ".") return data diff --git a/satpy/tests/reader_tests/test_fci_l1c_fdhsi.py b/satpy/tests/reader_tests/test_fci_l1c_fdhsi.py index d649cae780..e1bbdcd253 100644 --- a/satpy/tests/reader_tests/test_fci_l1c_fdhsi.py +++ b/satpy/tests/reader_tests/test_fci_l1c_fdhsi.py @@ -81,6 +81,7 @@ def _get_test_content_for_channel(self, pat, ch): "add_offset": 10, "warm_scale_factor": 2, "warm_add_offset": -300, + "long_name": "Effective Radiance", "units": "mW.m-2.sr-1.(cm-1)-1", } ) @@ -92,6 +93,9 @@ def _get_test_content_for_channel(self, pat, ch): "valid_range": [0, 4095], "scale_factor": 5, "add_offset": 10, + "warm_scale_factor": 1, + "warm_add_offset": 0, + "long_name": "Effective Radiance", "units": "mW.m-2.sr-1.(cm-1)-1", } ) From 5fd08f0e839d320dfa45276de1a424bd9f093255 Mon Sep 17 00:00:00 2001 From: andream Date: Mon, 4 May 2020 17:57:10 +0000 Subject: [PATCH 09/14] split calibration further --- satpy/readers/fci_l1c_fdhsi.py | 54 +++++++++++++++++++--------------- 1 file changed, 30 insertions(+), 24 deletions(-) diff --git a/satpy/readers/fci_l1c_fdhsi.py b/satpy/readers/fci_l1c_fdhsi.py index db6133b2b4..24e5edc471 100644 --- a/satpy/readers/fci_l1c_fdhsi.py +++ b/satpy/readers/fci_l1c_fdhsi.py @@ -261,37 +261,43 @@ def get_area_def(self, key, info=None): self._cache[key.resolution] = area return area + def calibrate_counts(self, data, key, measured, root): + """Calibrate counts.""" + # counts to radiance scaling + original_radiance_units = data.attrs.get("units") + if key.name == 'ir_38': + data = xr.where(((2 ** 12 - 1 < data) & (data <= 2 ** 13 - 1)), + (data * data.attrs.get("warm_scale_factor", 1) + + data.attrs.get("warm_add_offset", 0)), + (data * data.attrs.get("scale_factor", 1) + + data.attrs.get("add_offset", 0)) + ) + else: + data = (data * data.attrs.get("scale_factor", 1) + + data.attrs.get("add_offset", 0)) + + data.attrs["units"] = original_radiance_units + + if key.calibration == 'brightness_temperature': + data = self._ir_calibrate(data, measured, root) + elif key.calibration == 'reflectance': + data = self._vis_calibrate(data, measured) + + return data + def calibrate(self, data, key, measured, root): """Calibrate data.""" if key.calibration == "counts": # from package description, this just means not applying add_offset # and scale_factor data.attrs["units"] = "1" + elif key.calibration in ['brightness_temperature', 'reflectance', 'radiance']: + data = self.calibrate_counts(data, key, measured, root) else: - original_radiance_units = data.attrs.get("units") - # counts to radiance scaling - if key.name == 'ir_38': - data = xr.where(((2**12-1 < data) & (data <= 2**13-1)), - (data * data.attrs.get("warm_scale_factor", 1) + - data.attrs.get("warm_add_offset", 0)), - (data * data.attrs.get("scale_factor", 1) + - data.attrs.get("add_offset", 0)) - ) - else: - data = (data * data.attrs.get("scale_factor", 1) + - data.attrs.get("add_offset", 0)) - - if key.calibration == 'brightness_temperature': - data = self._ir_calibrate(data, measured, root) - elif key.calibration == 'reflectance': - data = self._vis_calibrate(data, measured) - elif key.calibration == 'radiance': - data.attrs["units"] = original_radiance_units - else: - logger.error( - "Received unknown calibration key. Expected " - "'brightness_temperature', 'reflectance' or 'radiance, got " - + key.calibration + ".") + logger.error( + "Received unknown calibration key. Expected " + "'brightness_temperature', 'reflectance' or 'radiance, got " + + key.calibration + ".") return data From c34d8430a847f6b2bf7f486d75e74ac102db3a5b Mon Sep 17 00:00:00 2001 From: andream Date: Mon, 4 May 2020 18:59:45 +0000 Subject: [PATCH 10/14] split counts calibration even further --- satpy/readers/fci_l1c_fdhsi.py | 35 ++++++++++++++++++++++------------ 1 file changed, 23 insertions(+), 12 deletions(-) diff --git a/satpy/readers/fci_l1c_fdhsi.py b/satpy/readers/fci_l1c_fdhsi.py index 24e5edc471..a7668d5fc5 100644 --- a/satpy/readers/fci_l1c_fdhsi.py +++ b/satpy/readers/fci_l1c_fdhsi.py @@ -148,7 +148,7 @@ def get_dataset(self, key, info=None): data = data.where(data >= vr[0], nfv) data = data.where(data <= vr[1], nfv) - res = self.calibrate(data, key, measured, root) + res = self.calibrate(data, key) # pre-calibration units no longer apply info.pop("units") @@ -261,10 +261,9 @@ def get_area_def(self, key, info=None): self._cache[key.resolution] = area return area - def calibrate_counts(self, data, key, measured, root): - """Calibrate counts.""" - # counts to radiance scaling - original_radiance_units = data.attrs.get("units") + def calibrate_to_radiances(self, data, key): + """Calibrate counts to radiances.""" + radiance_units = data.attrs["units"] if key.name == 'ir_38': data = xr.where(((2 ** 12 - 1 < data) & (data <= 2 ** 13 - 1)), (data * data.attrs.get("warm_scale_factor", 1) + @@ -276,23 +275,31 @@ def calibrate_counts(self, data, key, measured, root): data = (data * data.attrs.get("scale_factor", 1) + data.attrs.get("add_offset", 0)) - data.attrs["units"] = original_radiance_units + data.attrs["units"] = radiance_units + + return data + + def calibrate_counts(self, data, key): + """Calibrate counts to radiances, brightness temperatures, or reflectances.""" + # counts to radiance scaling + + data = self.calibrate_to_radiances(data, key) if key.calibration == 'brightness_temperature': - data = self._ir_calibrate(data, measured, root) + data = self._ir_calibrate(data, key) elif key.calibration == 'reflectance': - data = self._vis_calibrate(data, measured) + data = self._vis_calibrate(data, key) return data - def calibrate(self, data, key, measured, root): + def calibrate(self, data, key): """Calibrate data.""" if key.calibration == "counts": # from package description, this just means not applying add_offset # and scale_factor data.attrs["units"] = "1" elif key.calibration in ['brightness_temperature', 'reflectance', 'radiance']: - data = self.calibrate_counts(data, key, measured, root) + data = self.calibrate_counts(data, key) else: logger.error( "Received unknown calibration key. Expected " @@ -301,8 +308,10 @@ def calibrate(self, data, key, measured, root): return data - def _ir_calibrate(self, radiance, measured, root): + def _ir_calibrate(self, radiance, key): """IR channel calibration.""" + measured, root = self.get_channel_dataset(key.name) + # using the method from RADTOBR and PUG vc = self[measured + "/radiance_to_bt_conversion_coefficient_wavenumber"] @@ -336,8 +345,10 @@ def _ir_calibrate(self, radiance, measured, root): res.attrs["units"] = "K" return res - def _vis_calibrate(self, radiance, measured): + def _vis_calibrate(self, radiance, key): """VIS channel calibration.""" + measured, _ = self.get_channel_dataset(key.name) + cesi = self[measured + "/channel_effective_solar_irradiance"] if cesi == cesi.attrs.get( From eb3674426bd5a79a7ff3d239f0065e98476f6f6d Mon Sep 17 00:00:00 2001 From: andream Date: Mon, 4 May 2020 19:38:51 +0000 Subject: [PATCH 11/14] simplify nan return if coefficients are invalid --- satpy/readers/fci_l1c_fdhsi.py | 17 ++--------------- 1 file changed, 2 insertions(+), 15 deletions(-) diff --git a/satpy/readers/fci_l1c_fdhsi.py b/satpy/readers/fci_l1c_fdhsi.py index a7668d5fc5..1182d6561f 100644 --- a/satpy/readers/fci_l1c_fdhsi.py +++ b/satpy/readers/fci_l1c_fdhsi.py @@ -66,7 +66,6 @@ import logging import numpy as np -import dask.array as da import xarray as xr from pyresample import geometry @@ -330,13 +329,7 @@ def _ir_calibrate(self, radiance, key): v.attrs.get("long_name", "at least one necessary coefficient"), root)) - return xr.DataArray( - da.full(shape=radiance.shape, - chunks=radiance.chunks, - fill_value=np.nan), - dims=radiance.dims, - coords=radiance.coords, - attrs=radiance.attrs) + return radiance*np.nan nom = c2 * vc denom = a * np.log(1 + (c1 * vc**3) / radiance) @@ -356,13 +349,7 @@ def _vis_calibrate(self, radiance, key): logger.error( "channel effective solar irradiance set to fill value, " "cannot produce reflectance for {:s}.".format(measured)) - return xr.DataArray( - da.full(shape=radiance.shape, - chunks=radiance.chunks, - fill_value=np.nan), - dims=radiance.dims, - coords=radiance.coords, - attrs=radiance.attrs) + return radiance*np.nan sun_earth_distance = np.mean(self["state/celestial/earth_sun_distance"]) / 149597870.7 # [AU] From 565c21d079126ddebfa0cefbe1754e51f91a0444 Mon Sep 17 00:00:00 2001 From: andream Date: Mon, 4 May 2020 19:42:42 +0000 Subject: [PATCH 12/14] add missing quote to 'radiance' --- satpy/readers/fci_l1c_fdhsi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/satpy/readers/fci_l1c_fdhsi.py b/satpy/readers/fci_l1c_fdhsi.py index 1182d6561f..493cc34c47 100644 --- a/satpy/readers/fci_l1c_fdhsi.py +++ b/satpy/readers/fci_l1c_fdhsi.py @@ -302,7 +302,7 @@ def calibrate(self, data, key): else: logger.error( "Received unknown calibration key. Expected " - "'brightness_temperature', 'reflectance' or 'radiance, got " + "'brightness_temperature', 'reflectance' or 'radiance', got " + key.calibration + ".") return data From f78dd3a290d8db9305ef7273609bec875bc3bb06 Mon Sep 17 00:00:00 2001 From: andream Date: Tue, 5 May 2020 07:57:36 +0000 Subject: [PATCH 13/14] rename calibration functions, rename get_channel_dataset function --- satpy/readers/fci_l1c_fdhsi.py | 86 +++++++++++++++++----------------- 1 file changed, 42 insertions(+), 44 deletions(-) diff --git a/satpy/readers/fci_l1c_fdhsi.py b/satpy/readers/fci_l1c_fdhsi.py index 493cc34c47..a8ef233e8c 100644 --- a/satpy/readers/fci_l1c_fdhsi.py +++ b/satpy/readers/fci_l1c_fdhsi.py @@ -128,9 +128,8 @@ def get_dataset(self, key, info=None): logger.debug('Reading {} from {}'.format(key.name, self.filename)) # Get the dataset # Get metadata for given dataset - measured, root = self.get_channel_dataset(key.name) - radlab = measured + "/effective_radiance" - data = self[radlab] + measured = self.get_channel_group_path(key.name) + data = self[measured + "/effective_radiance"] attrs = data.attrs.copy() info = info.copy() @@ -172,17 +171,16 @@ def get_dataset(self, key, info=None): return res - def get_channel_dataset(self, channel): - """Get channel dataset.""" - root_group = 'data/{}'.format(channel) - group = 'data/{}/measured'.format(channel) + def get_channel_group_path(self, channel): + """Get channel group path.""" + group_path = 'data/{}/measured'.format(channel) - return group, root_group + return group_path def calc_area_extent(self, key): """Calculate area extent for a dataset.""" # Get metadata for given dataset - measured, root = self.get_channel_dataset(key.name) + measured = self.get_channel_group_path(key.name) # Get start/end line and column of loaded swath. nlines, ncols = self[measured + "/effective_radiance/shape"] @@ -260,56 +258,56 @@ def get_area_def(self, key, info=None): self._cache[key.resolution] = area return area - def calibrate_to_radiances(self, data, key): - """Calibrate counts to radiances.""" - radiance_units = data.attrs["units"] - if key.name == 'ir_38': - data = xr.where(((2 ** 12 - 1 < data) & (data <= 2 ** 13 - 1)), - (data * data.attrs.get("warm_scale_factor", 1) + - data.attrs.get("warm_add_offset", 0)), - (data * data.attrs.get("scale_factor", 1) + - data.attrs.get("add_offset", 0)) - ) + def calibrate(self, data, key): + """Calibrate data.""" + if key.calibration == "counts": + # from package description, this just means not applying add_offset + # and scale_factor + data.attrs["units"] = "1" + elif key.calibration in ['brightness_temperature', 'reflectance', 'radiance']: + data = self.calibrate_counts_to_physical_quantity(data, key) else: - data = (data * data.attrs.get("scale_factor", 1) + - data.attrs.get("add_offset", 0)) - - data.attrs["units"] = radiance_units + logger.error( + "Received unknown calibration key. Expected " + "'brightness_temperature', 'reflectance' or 'radiance', got " + + key.calibration + ".") return data - def calibrate_counts(self, data, key): + def calibrate_counts_to_physical_quantity(self, data, key): """Calibrate counts to radiances, brightness temperatures, or reflectances.""" # counts to radiance scaling - data = self.calibrate_to_radiances(data, key) + data = self.calibrate_counts_to_rad(data, key) if key.calibration == 'brightness_temperature': - data = self._ir_calibrate(data, key) + data = self.calibrate_rad_to_bt(data, key) elif key.calibration == 'reflectance': - data = self._vis_calibrate(data, key) + data = self.calibrate_rad_to_refl(data, key) return data - def calibrate(self, data, key): - """Calibrate data.""" - if key.calibration == "counts": - # from package description, this just means not applying add_offset - # and scale_factor - data.attrs["units"] = "1" - elif key.calibration in ['brightness_temperature', 'reflectance', 'radiance']: - data = self.calibrate_counts(data, key) + def calibrate_counts_to_rad(self, data, key): + """Calibrate counts to radiances.""" + radiance_units = data.attrs["units"] + if key.name == 'ir_38': + data = xr.where(((2 ** 12 - 1 < data) & (data <= 2 ** 13 - 1)), + (data * data.attrs.get("warm_scale_factor", 1) + + data.attrs.get("warm_add_offset", 0)), + (data * data.attrs.get("scale_factor", 1) + + data.attrs.get("add_offset", 0)) + ) else: - logger.error( - "Received unknown calibration key. Expected " - "'brightness_temperature', 'reflectance' or 'radiance', got " - + key.calibration + ".") + data = (data * data.attrs.get("scale_factor", 1) + + data.attrs.get("add_offset", 0)) + + data.attrs["units"] = radiance_units return data - def _ir_calibrate(self, radiance, key): + def calibrate_rad_to_bt(self, radiance, key): """IR channel calibration.""" - measured, root = self.get_channel_dataset(key.name) + measured = self.get_channel_group_path(key.name) # using the method from RADTOBR and PUG vc = self[measured + "/radiance_to_bt_conversion_coefficient_wavenumber"] @@ -328,7 +326,7 @@ def _ir_calibrate(self, radiance, key): "brightness temperatures for {:s}.".format( v.attrs.get("long_name", "at least one necessary coefficient"), - root)) + measured)) return radiance*np.nan nom = c2 * vc @@ -338,9 +336,9 @@ def _ir_calibrate(self, radiance, key): res.attrs["units"] = "K" return res - def _vis_calibrate(self, radiance, key): + def calibrate_rad_to_refl(self, radiance, key): """VIS channel calibration.""" - measured, _ = self.get_channel_dataset(key.name) + measured = self.get_channel_group_path(key.name) cesi = self[measured + "/channel_effective_solar_irradiance"] From bbf93d909e0abc1e9e6b96c50b0afd2d0cd42eaf Mon Sep 17 00:00:00 2001 From: andream Date: Tue, 5 May 2020 08:07:52 +0000 Subject: [PATCH 14/14] rename get_channel_group_path to get_channel_measured_group_path --- satpy/readers/fci_l1c_fdhsi.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/satpy/readers/fci_l1c_fdhsi.py b/satpy/readers/fci_l1c_fdhsi.py index a8ef233e8c..6116b3dc4a 100644 --- a/satpy/readers/fci_l1c_fdhsi.py +++ b/satpy/readers/fci_l1c_fdhsi.py @@ -128,7 +128,7 @@ def get_dataset(self, key, info=None): logger.debug('Reading {} from {}'.format(key.name, self.filename)) # Get the dataset # Get metadata for given dataset - measured = self.get_channel_group_path(key.name) + measured = self.get_channel_measured_group_path(key.name) data = self[measured + "/effective_radiance"] attrs = data.attrs.copy() @@ -171,16 +171,16 @@ def get_dataset(self, key, info=None): return res - def get_channel_group_path(self, channel): - """Get channel group path.""" - group_path = 'data/{}/measured'.format(channel) + def get_channel_measured_group_path(self, channel): + """Get the channel's measured group path.""" + measured_group_path = 'data/{}/measured'.format(channel) - return group_path + return measured_group_path def calc_area_extent(self, key): """Calculate area extent for a dataset.""" # Get metadata for given dataset - measured = self.get_channel_group_path(key.name) + measured = self.get_channel_measured_group_path(key.name) # Get start/end line and column of loaded swath. nlines, ncols = self[measured + "/effective_radiance/shape"] @@ -307,7 +307,7 @@ def calibrate_counts_to_rad(self, data, key): def calibrate_rad_to_bt(self, radiance, key): """IR channel calibration.""" - measured = self.get_channel_group_path(key.name) + measured = self.get_channel_measured_group_path(key.name) # using the method from RADTOBR and PUG vc = self[measured + "/radiance_to_bt_conversion_coefficient_wavenumber"] @@ -338,7 +338,7 @@ def calibrate_rad_to_bt(self, radiance, key): def calibrate_rad_to_refl(self, radiance, key): """VIS channel calibration.""" - measured = self.get_channel_group_path(key.name) + measured = self.get_channel_measured_group_path(key.name) cesi = self[measured + "/channel_effective_solar_irradiance"]