Skip to content

Commit

Permalink
Merge pull request #2637 from pnuu/fci-32-bit-datasets
Browse files Browse the repository at this point in the history 
Keep FCI data as 32-bit floats
  • Loading branch information
@mraspaud
mraspaud committed Nov 16, 2023
2 parents 35f9f96 + 1930845 commit a0b8a7b
Show file tree
Hide file tree
Showing 2 changed files with 35 additions and 32 deletions.
28 changes: 14 additions & 14 deletions satpy/readers/fci_l1c_nc.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -330,13 +330,13 @@ def _get_dataset_measurand(self, key, info=None):

fv = attrs.pop(
"FillValue",
default_fillvals.get(data.dtype.str[1:], np.nan))
vr = attrs.get("valid_range", [-np.inf, np.inf])
default_fillvals.get(data.dtype.str[1:], np.float32(np.nan)))
vr = attrs.get("valid_range", [np.float32(-np.inf), np.float32(np.inf)])
if key["calibration"] == "counts":
attrs["_FillValue"] = fv
nfv = data.dtype.type(fv)
else:
nfv = np.nan
nfv = np.float32(np.nan)
data = data.where((data >= vr[0]) & (data <= vr[1]), nfv)

res = self.calibrate(data, key)
Expand Down Expand Up @@ -632,16 +632,15 @@ def calibrate_counts_to_rad(self, data, key):
def calibrate_rad_to_bt(self, radiance, key):
"""IR channel calibration."""
# using the method from PUG section Converting from Effective Radiance to Brightness Temperature for IR Channels

measured = self.get_channel_measured_group_path(key["name"])

vc = self.get_and_cache_npxr(measured + "/radiance_to_bt_conversion_coefficient_wavenumber")
vc = self.get_and_cache_npxr(measured + "/radiance_to_bt_conversion_coefficient_wavenumber").astype(np.float32)

a = self.get_and_cache_npxr(measured + "/radiance_to_bt_conversion_coefficient_a")
b = self.get_and_cache_npxr(measured + "/radiance_to_bt_conversion_coefficient_b")
a = self.get_and_cache_npxr(measured + "/radiance_to_bt_conversion_coefficient_a").astype(np.float32)
b = self.get_and_cache_npxr(measured + "/radiance_to_bt_conversion_coefficient_b").astype(np.float32)

c1 = self.get_and_cache_npxr(measured + "/radiance_to_bt_conversion_constant_c1")
c2 = self.get_and_cache_npxr(measured + "/radiance_to_bt_conversion_constant_c2")
c1 = self.get_and_cache_npxr(measured + "/radiance_to_bt_conversion_constant_c1").astype(np.float32)
c2 = self.get_and_cache_npxr(measured + "/radiance_to_bt_conversion_constant_c2").astype(np.float32)

for v in (vc, a, b, c1, c2):
if v == v.attrs.get("FillValue",
Expand All @@ -652,26 +651,27 @@ def calibrate_rad_to_bt(self, radiance, key):
v.attrs.get("long_name",
"at least one necessary coefficient"),
measured))
return radiance * np.nan
return radiance * np.float32(np.nan)

nom = c2 * vc
denom = a * np.log(1 + (c1 * vc ** 3) / radiance)
denom = a * np.log(1 + (c1 * vc ** np.float32(3.)) / radiance)

res = nom / denom - b / a

return res

def calibrate_rad_to_refl(self, radiance, key):
"""VIS channel calibration."""
measured = self.get_channel_measured_group_path(key["name"])

cesi = self.get_and_cache_npxr(measured + "/channel_effective_solar_irradiance")
cesi = self.get_and_cache_npxr(measured + "/channel_effective_solar_irradiance").astype(np.float32)

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 radiance * np.nan
return radiance * np.float32(np.nan)

sun_earth_distance = np.mean(
self.get_and_cache_npxr("state/celestial/earth_sun_distance")) / 149597870.7 # [AU]
Expand All @@ -683,7 +683,7 @@ def calibrate_rad_to_refl(self, radiance, key):
"".format(sun_earth_distance))
sun_earth_distance = 1

res = 100 * radiance * np.pi * sun_earth_distance ** 2 / cesi
res = 100 * radiance * np.float32(np.pi) * np.float32(sun_earth_distance) ** np.float32(2) / cesi
return res


Expand Down
39 changes: 21 additions & 18 deletions satpy/tests/reader_tests/test_fci_l1c_nc.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -103,18 +103,21 @@ def _get_test_calib_for_channel_ir(data, meas_path):
from pyspectral.blackbody import C_SPEED as c
from pyspectral.blackbody import H_PLANCK as h
from pyspectral.blackbody import K_BOLTZMANN as k
data[meas_path + "/radiance_to_bt_conversion_coefficient_wavenumber"] = FakeH5Variable(da.array(955))
data[meas_path + "/radiance_to_bt_conversion_coefficient_a"] = FakeH5Variable(da.array(1))
data[meas_path + "/radiance_to_bt_conversion_coefficient_b"] = FakeH5Variable(da.array(0.4))
data[meas_path + "/radiance_to_bt_conversion_constant_c1"] = FakeH5Variable(da.array(1e11 * 2 * h * c ** 2))
data[meas_path + "/radiance_to_bt_conversion_constant_c2"] = FakeH5Variable(da.array(1e2 * h * c / k))
data[meas_path + "/radiance_to_bt_conversion_coefficient_wavenumber"] = FakeH5Variable(
da.array(955.0, dtype=np.float32))
data[meas_path + "/radiance_to_bt_conversion_coefficient_a"] = FakeH5Variable(da.array(1.0, dtype=np.float32))
data[meas_path + "/radiance_to_bt_conversion_coefficient_b"] = FakeH5Variable(da.array(0.4, dtype=np.float32))
data[meas_path + "/radiance_to_bt_conversion_constant_c1"] = FakeH5Variable(
da.array(1e11 * 2 * h * c ** 2, dtype=np.float32))
data[meas_path + "/radiance_to_bt_conversion_constant_c2"] = FakeH5Variable(
da.array(1e2 * h * c / k, dtype=np.float32))
return data


def _get_test_calib_for_channel_vis(data, meas):
data["state/celestial/earth_sun_distance"] = FakeH5Variable(
da.repeat(da.array([149597870.7]), 6000), dims=("x"))
data[meas + "/channel_effective_solar_irradiance"] = FakeH5Variable(da.array(50))
data[meas + "/channel_effective_solar_irradiance"] = FakeH5Variable(da.array((50.0), dtype=np.float32))
return data


Expand All @@ -124,7 +127,7 @@ def _get_test_calib_data_for_channel(data, ch_str):
_get_test_calib_for_channel_ir(data, meas_path)
elif ch_str.startswith("vis") or ch_str.startswith("nir"):
_get_test_calib_for_channel_vis(data, meas_path)
data[meas_path + "/radiance_unit_conversion_coefficient"] = xr.DataArray(da.array(1234.56))
data[meas_path + "/radiance_unit_conversion_coefficient"] = xr.DataArray(da.array(1234.56, dtype=np.float32))


def _get_test_image_data_for_channel(data, ch_str, n_rows_cols):
Expand All @@ -145,8 +148,8 @@ def _get_test_image_data_for_channel(data, ch_str, n_rows_cols):
dims=("y", "x"),
attrs={
"valid_range": [0, 8191],
"warm_scale_factor": 2,
"warm_add_offset": -300,
"warm_scale_factor": np.float32(2.0),
"warm_add_offset": np.float32(-300.0),
**common_attrs
}
)
Expand All @@ -156,8 +159,8 @@ def _get_test_image_data_for_channel(data, ch_str, n_rows_cols):
dims=("y", "x"),
attrs={
"valid_range": [0, 4095],
"warm_scale_factor": 1,
"warm_add_offset": 0,
"warm_scale_factor": np.float32(1.0),
"warm_add_offset": np.float32(0.0),
**common_attrs
}
)
Expand Down Expand Up @@ -521,10 +524,10 @@ def test_load_radiance(self, reader_configs, fh_param,
fh_param["channels"]["terran_grid_type"]):
assert res[ch].shape == (GRID_TYPE_INFO_FOR_TEST_CONTENT[grid_type]["nrows"],
GRID_TYPE_INFO_FOR_TEST_CONTENT[grid_type]["ncols"])
assert res[ch].dtype == np.float64
assert res[ch].dtype == np.float32
assert res[ch].attrs["calibration"] == "radiance"
assert res[ch].attrs["units"] == "mW m-2 sr-1 (cm-1)-1"
assert res[ch].attrs["radiance_unit_conversion_coefficient"] == 1234.56
assert res[ch].attrs["radiance_unit_conversion_coefficient"].values == np.float32(1234.56)
if ch == "ir_38":
numpy.testing.assert_array_equal(res[ch][-1], 15)
numpy.testing.assert_array_equal(res[ch][0], 9700)
Expand All @@ -544,7 +547,7 @@ def test_load_reflectance(self, reader_configs, fh_param,
for ch, grid_type in zip(fh_param["channels"]["solar"], fh_param["channels"]["solar_grid_type"]):
assert res[ch].shape == (GRID_TYPE_INFO_FOR_TEST_CONTENT[grid_type]["nrows"],
GRID_TYPE_INFO_FOR_TEST_CONTENT[grid_type]["ncols"])
assert res[ch].dtype == np.float64
assert res[ch].dtype == np.float32
assert res[ch].attrs["calibration"] == "reflectance"
assert res[ch].attrs["units"] == "%"
numpy.testing.assert_array_almost_equal(res[ch], 100 * 15 * 1 * np.pi / 50)
Expand All @@ -564,15 +567,15 @@ def test_load_bt(self, reader_configs, caplog, fh_param,
for ch, grid_type in zip(fh_param["channels"]["terran"], fh_param["channels"]["terran_grid_type"]):
assert res[ch].shape == (GRID_TYPE_INFO_FOR_TEST_CONTENT[grid_type]["nrows"],
GRID_TYPE_INFO_FOR_TEST_CONTENT[grid_type]["ncols"])
assert res[ch].dtype == np.float64
assert res[ch].dtype == np.float32
assert res[ch].attrs["calibration"] == "brightness_temperature"
assert res[ch].attrs["units"] == "K"

if ch == "ir_38":
numpy.testing.assert_array_almost_equal(res[ch][-1], 209.68274099)
numpy.testing.assert_array_almost_equal(res[ch][0], 1888.851296)
numpy.testing.assert_array_almost_equal(res[ch][-1], np.float32(209.68275))
numpy.testing.assert_array_almost_equal(res[ch][0], np.float32(1888.8513))
else:
numpy.testing.assert_array_almost_equal(res[ch], 209.68274099)
numpy.testing.assert_array_almost_equal(res[ch], np.float32(209.68275))

@pytest.mark.parametrize("fh_param", [(lazy_fixture("FakeFCIFileHandlerFDHSI_fixture")),
(lazy_fixture("FakeFCIFileHandlerHRFI_fixture"))])
Expand Down

0 comments on commit a0b8a7b

Please sign in to comment.