From c642b3afa495a99e925d5805f9af57d86501fab1 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Fri, 11 Mar 2022 10:19:07 -0600 Subject: [PATCH 01/13] Add option to replace saturated MODIS L1b band 2 values with max reflectance --- satpy/readers/hdfeos_base.py | 6 +++--- satpy/readers/modis_l1b.py | 19 ++++++++++------ satpy/tests/reader_tests/_modis_fixtures.py | 7 ++++++ satpy/tests/reader_tests/test_modis_l1b.py | 24 +++++++++++++++++++++ 4 files changed, 47 insertions(+), 9 deletions(-) diff --git a/satpy/readers/hdfeos_base.py b/satpy/readers/hdfeos_base.py index 731611765c..cd47daa71f 100644 --- a/satpy/readers/hdfeos_base.py +++ b/satpy/readers/hdfeos_base.py @@ -91,7 +91,7 @@ def _find_and_run_interpolation(interpolation_functions, src_resolution, dst_res class HDFEOSBaseFileReader(BaseFileHandler): """Base file handler for HDF EOS data for both L1b and L2 products.""" - def __init__(self, filename, filename_info, filetype_info): + def __init__(self, filename, filename_info, filetype_info, **kwargs): """Initialize the base reader.""" BaseFileHandler.__init__(self, filename, filename_info, filetype_info) try: @@ -291,9 +291,9 @@ class HDFEOSGeoReader(HDFEOSBaseFileReader): 'solar_zenith_angle': ('SolarZenith', 'Solar_Zenith'), } - def __init__(self, filename, filename_info, filetype_info): + def __init__(self, filename, filename_info, filetype_info, **kwargs): """Initialize the geographical reader.""" - HDFEOSBaseFileReader.__init__(self, filename, filename_info, filetype_info) + HDFEOSBaseFileReader.__init__(self, filename, filename_info, filetype_info, **kwargs) self.cache = {} @staticmethod diff --git a/satpy/readers/modis_l1b.py b/satpy/readers/modis_l1b.py index f0b5a74e8a..509343431d 100644 --- a/satpy/readers/modis_l1b.py +++ b/satpy/readers/modis_l1b.py @@ -62,9 +62,10 @@ class HDFEOSBandReader(HDFEOSBaseFileReader): "Q": 250, "H": 500} - def __init__(self, filename, filename_info, filetype_info): + def __init__(self, filename, filename_info, filetype_info, mask_saturated=True, **kwargs): """Init the file handler.""" - HDFEOSBaseFileReader.__init__(self, filename, filename_info, filetype_info) + HDFEOSBaseFileReader.__init__(self, filename, filename_info, filetype_info, **kwargs) + self._mask_saturated = mask_saturated ds = self.metadata['INVENTORYMETADATA'][ 'COLLECTIONDESCRIPTIONCLASS']['SHORTNAME']['VALUE'] @@ -121,7 +122,13 @@ def get_dataset(self, key, info): # 65500 NAD closed upper limit array = array.where(array >= np.float32(valid_range[0])) - array = array.where(array <= np.float32(valid_range[1])) + if self._mask_saturated: + array = array.where(array <= np.float32(valid_range[1])) + else: + valid_max = np.float32(valid_range[1]) + array = array.where((array != 65533) & (array != 65528), valid_max) + array = array.where(array <= valid_max) + array = array.where(from_sds(uncertainty, chunks=CHUNK_SIZE)[index, :, :] < 15) if key['calibration'] == 'brightness_temperature': @@ -182,10 +189,10 @@ def get_dataset(self, key, info): class MixedHDFEOSReader(HDFEOSGeoReader, HDFEOSBandReader): """A file handler for the files that have both regular bands and geographical information in them.""" - def __init__(self, filename, filename_info, filetype_info): + def __init__(self, filename, filename_info, filetype_info, **kwargs): """Init the file handler.""" - HDFEOSGeoReader.__init__(self, filename, filename_info, filetype_info) - HDFEOSBandReader.__init__(self, filename, filename_info, filetype_info) + HDFEOSGeoReader.__init__(self, filename, filename_info, filetype_info, **kwargs) + HDFEOSBandReader.__init__(self, filename, filename_info, filetype_info, **kwargs) def get_dataset(self, key, info): """Get the dataset.""" diff --git a/satpy/tests/reader_tests/_modis_fixtures.py b/satpy/tests/reader_tests/_modis_fixtures.py index 29a2c34842..1fdb596046 100644 --- a/satpy/tests/reader_tests/_modis_fixtures.py +++ b/satpy/tests/reader_tests/_modis_fixtures.py @@ -75,6 +75,13 @@ def _generate_angle_data(resolution: int) -> np.ndarray: def _generate_visible_data(resolution: int, num_bands: int, dtype=np.uint16) -> np.ndarray: shape = _shape_for_resolution(resolution) data = np.zeros((num_bands, shape[0], shape[1]), dtype=dtype) + + # add fill value to every band + data[:, -1, -1] = 65535 + + # add band 2 saturation and can't aggregate fill values + data[1, -1, -2] = 65533 + data[1, -1, -3] = 65528 return data diff --git a/satpy/tests/reader_tests/test_modis_l1b.py b/satpy/tests/reader_tests/test_modis_l1b.py index 981fac39bb..245eb91395 100644 --- a/satpy/tests/reader_tests/test_modis_l1b.py +++ b/satpy/tests/reader_tests/test_modis_l1b.py @@ -162,3 +162,27 @@ def test_load_vis(self, modis_l1b_nasa_mod021km_file): assert dataset.shape == _shape_for_resolution(1000) assert dataset.attrs['resolution'] == 1000 _check_shared_metadata(dataset) + + @pytest.mark.parametrize("mask_saturated", [False, True]) + def test_load_vis_saturation(self, mask_saturated, modis_l1b_nasa_mod021km_file): + """Test loading visible band.""" + scene = Scene(reader='modis_l1b', filenames=modis_l1b_nasa_mod021km_file, + reader_kwargs={"mask_saturated": mask_saturated}) + dataset_name = '2' + scene.load([dataset_name]) + dataset = scene[dataset_name] + assert dataset.shape == _shape_for_resolution(1000) + assert dataset.attrs['resolution'] == 1000 + _check_shared_metadata(dataset) + + # check saturation fill values + data = dataset.values + assert np.isnan(data[-1, -1]) # normal fill value + if mask_saturated: + assert np.isnan(data[-1, -2]) # saturation + assert np.isnan(data[-1, -3]) # can't aggregate + else: + # test data factor/offset are 1/0 + # albedos are converted to % + assert data[-1, -2] >= 32767 * 100.0 # saturation + assert data[-1, -3] >= 32767 * 100.0 # can't aggregate From b68fe962229ae1e134e8b16f6a20c649549ef19f Mon Sep 17 00:00:00 2001 From: David Hoese Date: Fri, 11 Mar 2022 10:23:58 -0600 Subject: [PATCH 02/13] Refactor modis_l1b valid range checks --- satpy/readers/modis_l1b.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/satpy/readers/modis_l1b.py b/satpy/readers/modis_l1b.py index 509343431d..10988d39c7 100644 --- a/satpy/readers/modis_l1b.py +++ b/satpy/readers/modis_l1b.py @@ -121,14 +121,7 @@ def get_dataset(self, key, info): # 65501 - 65523 (reserved for future use) # 65500 NAD closed upper limit - array = array.where(array >= np.float32(valid_range[0])) - if self._mask_saturated: - array = array.where(array <= np.float32(valid_range[1])) - else: - valid_max = np.float32(valid_range[1]) - array = array.where((array != 65533) & (array != 65528), valid_max) - array = array.where(array <= valid_max) - + array = self._mask_invalid_or_fill_saturated(array, valid_range) array = array.where(from_sds(uncertainty, chunks=CHUNK_SIZE)[index, :, :] < 15) if key['calibration'] == 'brightness_temperature': @@ -185,6 +178,17 @@ def get_dataset(self, key, info): self._add_satpy_metadata(key, projectable) return projectable + def _mask_invalid_or_fill_saturated(self, array, valid_range): + valid_min = np.float32(valid_range[0]) + valid_max = np.float32(valid_range[1]) + array = array.where(array >= valid_min) + if self._mask_saturated: + array = array.where(array <= valid_max) + else: + array = array.where((array != 65533) & (array != 65528), valid_max) + array = array.where(array <= valid_max) + return array + class MixedHDFEOSReader(HDFEOSGeoReader, HDFEOSBandReader): """A file handler for the files that have both regular bands and geographical information in them.""" From fb5cfb4d5d7b59904e4906de04d812ba19254c97 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Fri, 11 Mar 2022 10:25:28 -0600 Subject: [PATCH 03/13] Move modis l1b fill value block comment to mask method --- satpy/readers/modis_l1b.py | 41 +++++++++++++++++++------------------- 1 file changed, 20 insertions(+), 21 deletions(-) diff --git a/satpy/readers/modis_l1b.py b/satpy/readers/modis_l1b.py index 10988d39c7..eb1c3fdbe7 100644 --- a/satpy/readers/modis_l1b.py +++ b/satpy/readers/modis_l1b.py @@ -100,27 +100,6 @@ def get_dataset(self, key, info): array = xr.DataArray(from_sds(subdata, chunks=CHUNK_SIZE)[index, :, :], dims=['y', 'x']).astype(np.float32) valid_range = var_attrs['valid_range'] - - # Fill values: - # Data Value Meaning - # 65535 Fill Value (includes reflective band data at night mode - # and completely missing L1A scans) - # 65534 L1A DN is missing within a scan - # 65533 Detector is saturated - # 65532 Cannot compute zero point DN, e.g., SV is saturated - # 65531 Detector is dead (see comments below) - # 65530 RSB dn** below the minimum of the scaling range - # 65529 TEB radiance or RSB dn** exceeds the maximum of the - # scaling range - # 65528 Aggregation algorithm failure - # 65527 Rotation of Earth view Sector from nominal science - # collection position - # 65526 Calibration coefficient b1 could not be computed - # 65525 Subframe is dead - # 65524 Both sides of the PCLW electronics on simultaneously - # 65501 - 65523 (reserved for future use) - # 65500 NAD closed upper limit - array = self._mask_invalid_or_fill_saturated(array, valid_range) array = array.where(from_sds(uncertainty, chunks=CHUNK_SIZE)[index, :, :] < 15) @@ -179,6 +158,26 @@ def get_dataset(self, key, info): return projectable def _mask_invalid_or_fill_saturated(self, array, valid_range): + """Replace fill values with NaN or replace saturation with max reflectance.""" + # Fill values: + # Data Value Meaning + # 65535 Fill Value (includes reflective band data at night mode + # and completely missing L1A scans) + # 65534 L1A DN is missing within a scan + # 65533 Detector is saturated + # 65532 Cannot compute zero point DN, e.g., SV is saturated + # 65531 Detector is dead (see comments below) + # 65530 RSB dn** below the minimum of the scaling range + # 65529 TEB radiance or RSB dn** exceeds the maximum of the + # scaling range + # 65528 Aggregation algorithm failure + # 65527 Rotation of Earth view Sector from nominal science + # collection position + # 65526 Calibration coefficient b1 could not be computed + # 65525 Subframe is dead + # 65524 Both sides of the PCLW electronics on simultaneously + # 65501 - 65523 (reserved for future use) + # 65500 NAD closed upper limit valid_min = np.float32(valid_range[0]) valid_max = np.float32(valid_range[1]) array = array.where(array >= valid_min) From e93ba67d7223b52c566bce7788eba1ba1ad2e0f1 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Fri, 11 Mar 2022 10:34:04 -0600 Subject: [PATCH 04/13] More modis l1b refactoring --- satpy/readers/modis_l1b.py | 61 +++++++++++++++++++++----------------- 1 file changed, 34 insertions(+), 27 deletions(-) diff --git a/satpy/readers/modis_l1b.py b/satpy/readers/modis_l1b.py index eb1c3fdbe7..fa86ae21b5 100644 --- a/satpy/readers/modis_l1b.py +++ b/satpy/readers/modis_l1b.py @@ -89,12 +89,10 @@ def get_dataset(self, key, info): for dataset in datasets: subdata = self.sd.select(dataset) var_attrs = subdata.attributes() - band_names = var_attrs["band_names"].split(",") - - # get the relative indices of the desired channel try: - index = band_names.index(key['name']) + index = self._get_band_index(var_attrs, key["name"]) except ValueError: + # no band dimension continue uncertainty = self.sd.select(dataset + "_Uncert_Indexes") array = xr.DataArray(from_sds(subdata, chunks=CHUNK_SIZE)[index, :, :], @@ -102,29 +100,7 @@ def get_dataset(self, key, info): valid_range = var_attrs['valid_range'] array = self._mask_invalid_or_fill_saturated(array, valid_range) array = array.where(from_sds(uncertainty, chunks=CHUNK_SIZE)[index, :, :] < 15) - - if key['calibration'] == 'brightness_temperature': - projectable = calibrate_bt(array, var_attrs, index, key['name']) - info.setdefault('units', 'K') - info.setdefault('standard_name', 'toa_brightness_temperature') - elif key['calibration'] == 'reflectance': - projectable = calibrate_refl(array, var_attrs, index) - info.setdefault('units', '%') - info.setdefault('standard_name', - 'toa_bidirectional_reflectance') - elif key['calibration'] == 'radiance': - projectable = calibrate_radiance(array, var_attrs, index) - info.setdefault('units', var_attrs.get('radiance_units')) - info.setdefault('standard_name', - 'toa_outgoing_radiance_per_unit_wavelength') - elif key['calibration'] == 'counts': - projectable = calibrate_counts(array, var_attrs, index) - info.setdefault('units', 'counts') - info.setdefault('standard_name', 'counts') # made up - else: - raise ValueError("Unknown calibration for " - "key: {}".format(key)) - projectable.attrs = info + projectable = self._calibrate_data(key, info, array, var_attrs, index) # if ((platform_name == 'Aqua' and key['name'] in ["6", "27", "36"]) or # (platform_name == 'Terra' and key['name'] in ["29"])): @@ -157,6 +133,12 @@ def get_dataset(self, key, info): self._add_satpy_metadata(key, projectable) return projectable + def _get_band_index(self, var_attrs, name): + """Get the relative indices of the desired channel.""" + band_names = var_attrs["band_names"].split(",") + index = band_names.index(name) + return index + def _mask_invalid_or_fill_saturated(self, array, valid_range): """Replace fill values with NaN or replace saturation with max reflectance.""" # Fill values: @@ -188,6 +170,31 @@ def _mask_invalid_or_fill_saturated(self, array, valid_range): array = array.where(array <= valid_max) return array + def _calibrate_data(self, key, info, array, var_attrs, index): + if key['calibration'] == 'brightness_temperature': + projectable = calibrate_bt(array, var_attrs, index, key['name']) + info.setdefault('units', 'K') + info.setdefault('standard_name', 'toa_brightness_temperature') + elif key['calibration'] == 'reflectance': + projectable = calibrate_refl(array, var_attrs, index) + info.setdefault('units', '%') + info.setdefault('standard_name', + 'toa_bidirectional_reflectance') + elif key['calibration'] == 'radiance': + projectable = calibrate_radiance(array, var_attrs, index) + info.setdefault('units', var_attrs.get('radiance_units')) + info.setdefault('standard_name', + 'toa_outgoing_radiance_per_unit_wavelength') + elif key['calibration'] == 'counts': + projectable = calibrate_counts(array, var_attrs, index) + info.setdefault('units', 'counts') + info.setdefault('standard_name', 'counts') # made up + else: + raise ValueError("Unknown calibration for " + "key: {}".format(key)) + projectable.attrs = info + return projectable + class MixedHDFEOSReader(HDFEOSGeoReader, HDFEOSBandReader): """A file handler for the files that have both regular bands and geographical information in them.""" From d25145fe2adb23fa5d963121be75e472543ba058 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Fri, 11 Mar 2022 10:35:37 -0600 Subject: [PATCH 05/13] Move possible var names for modis resolutions to class level --- satpy/readers/modis_l1b.py | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/satpy/readers/modis_l1b.py b/satpy/readers/modis_l1b.py index fa86ae21b5..4d6b83f0ed 100644 --- a/satpy/readers/modis_l1b.py +++ b/satpy/readers/modis_l1b.py @@ -62,6 +62,16 @@ class HDFEOSBandReader(HDFEOSBaseFileReader): "Q": 250, "H": 500} + res_to_possible_variable_names = { + 1000: ['EV_250_Aggr1km_RefSB', + 'EV_500_Aggr1km_RefSB', + 'EV_1KM_RefSB', + 'EV_1KM_Emissive'], + 500: ['EV_250_Aggr500_RefSB', + 'EV_500_RefSB'], + 250: ['EV_250_RefSB'], + } + def __init__(self, filename, filename_info, filetype_info, mask_saturated=True, **kwargs): """Init the file handler.""" HDFEOSBaseFileReader.__init__(self, filename, filename_info, filetype_info, **kwargs) @@ -73,19 +83,10 @@ def __init__(self, filename, filename_info, filetype_info, mask_saturated=True, def get_dataset(self, key, info): """Read data from file and return the corresponding projectables.""" - datadict = { - 1000: ['EV_250_Aggr1km_RefSB', - 'EV_500_Aggr1km_RefSB', - 'EV_1KM_RefSB', - 'EV_1KM_Emissive'], - 500: ['EV_250_Aggr500_RefSB', - 'EV_500_RefSB'], - 250: ['EV_250_RefSB']} - if self.resolution != key['resolution']: return - datasets = datadict[self.resolution] + datasets = self.res_to_possible_variable_names[self.resolution] for dataset in datasets: subdata = self.sd.select(dataset) var_attrs = subdata.attributes() From 3742cffe2d4b363afd789124208499e69b69ad15 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Fri, 11 Mar 2022 10:53:18 -0600 Subject: [PATCH 06/13] More refactoring of MODIS L1B code --- satpy/readers/modis_l1b.py | 92 +++++++++++++++++++++----------------- 1 file changed, 50 insertions(+), 42 deletions(-) diff --git a/satpy/readers/modis_l1b.py b/satpy/readers/modis_l1b.py index 4d6b83f0ed..6dc82f1bda 100644 --- a/satpy/readers/modis_l1b.py +++ b/satpy/readers/modis_l1b.py @@ -85,59 +85,62 @@ def get_dataset(self, key, info): """Read data from file and return the corresponding projectables.""" if self.resolution != key['resolution']: return + subdata, uncertainty, var_attrs, band_index = self._get_band_data_uncertainty_attrs_index(key["name"]) + array = xr.DataArray(from_sds(subdata, chunks=CHUNK_SIZE)[band_index, :, :], + dims=['y', 'x']).astype(np.float32) + valid_range = var_attrs['valid_range'] + array = self._mask_invalid_or_fill_saturated(array, valid_range) + array = self._mask_uncertain_pixels(array, uncertainty, band_index) + projectable = self._calibrate_data(key, info, array, var_attrs, band_index) + + # if ((platform_name == 'Aqua' and key['name'] in ["6", "27", "36"]) or + # (platform_name == 'Terra' and key['name'] in ["29"])): + # height, width = projectable.shape + # row_indices = projectable.mask.sum(1) == width + # if row_indices.sum() != height: + # projectable.mask[row_indices, :] = True + + # Get the orbit number + # if not satscene.orbit: + # mda = self.data.attributes()["CoreMetadata.0"] + # orbit_idx = mda.index("ORBITNUMBER") + # satscene.orbit = mda[orbit_idx + 111:orbit_idx + 116] + + # Trimming out dead sensor lines (detectors) on terra: + # (in addition channel 27, 30, 34, 35, and 36 are nosiy) + # if satscene.satname == "terra": + # for band in ["29"]: + # if not satscene[band].is_loaded() or satscene[band].data.mask.all(): + # continue + # width = satscene[band].data.shape[1] + # height = satscene[band].data.shape[0] + # indices = satscene[band].data.mask.sum(1) < width + # if indices.sum() == height: + # continue + # satscene[band] = satscene[band].data[indices, :] + # satscene[band].area = geometry.SwathDefinition( + # lons=satscene[band].area.lons[indices, :], + # lats=satscene[band].area.lats[indices, :]) + self._add_satpy_metadata(key, projectable) + return projectable + def _get_band_data_uncertainty_attrs_index(self, band_name): datasets = self.res_to_possible_variable_names[self.resolution] for dataset in datasets: subdata = self.sd.select(dataset) var_attrs = subdata.attributes() try: - index = self._get_band_index(var_attrs, key["name"]) + band_index = self._get_band_index(var_attrs, band_name) except ValueError: - # no band dimension + # can't find band in list of bands continue uncertainty = self.sd.select(dataset + "_Uncert_Indexes") - array = xr.DataArray(from_sds(subdata, chunks=CHUNK_SIZE)[index, :, :], - dims=['y', 'x']).astype(np.float32) - valid_range = var_attrs['valid_range'] - array = self._mask_invalid_or_fill_saturated(array, valid_range) - array = array.where(from_sds(uncertainty, chunks=CHUNK_SIZE)[index, :, :] < 15) - projectable = self._calibrate_data(key, info, array, var_attrs, index) - - # if ((platform_name == 'Aqua' and key['name'] in ["6", "27", "36"]) or - # (platform_name == 'Terra' and key['name'] in ["29"])): - # height, width = projectable.shape - # row_indices = projectable.mask.sum(1) == width - # if row_indices.sum() != height: - # projectable.mask[row_indices, :] = True - - # Get the orbit number - # if not satscene.orbit: - # mda = self.data.attributes()["CoreMetadata.0"] - # orbit_idx = mda.index("ORBITNUMBER") - # satscene.orbit = mda[orbit_idx + 111:orbit_idx + 116] - - # Trimming out dead sensor lines (detectors) on terra: - # (in addition channel 27, 30, 34, 35, and 36 are nosiy) - # if satscene.satname == "terra": - # for band in ["29"]: - # if not satscene[band].is_loaded() or satscene[band].data.mask.all(): - # continue - # width = satscene[band].data.shape[1] - # height = satscene[band].data.shape[0] - # indices = satscene[band].data.mask.sum(1) < width - # if indices.sum() == height: - # continue - # satscene[band] = satscene[band].data[indices, :] - # satscene[band].area = geometry.SwathDefinition( - # lons=satscene[band].area.lons[indices, :], - # lats=satscene[band].area.lats[indices, :]) - self._add_satpy_metadata(key, projectable) - return projectable - - def _get_band_index(self, var_attrs, name): + return subdata, uncertainty, var_attrs, band_index + + def _get_band_index(self, var_attrs, band_name): """Get the relative indices of the desired channel.""" band_names = var_attrs["band_names"].split(",") - index = band_names.index(name) + index = band_names.index(band_name) return index def _mask_invalid_or_fill_saturated(self, array, valid_range): @@ -171,6 +174,11 @@ def _mask_invalid_or_fill_saturated(self, array, valid_range): array = array.where(array <= valid_max) return array + def _mask_uncertain_pixels(self, array, uncertainty, band_index): + band_uncertainty = from_sds(uncertainty, chunks=CHUNK_SIZE)[band_index, :, :] + array = array.where(band_uncertainty < 15) + return array + def _calibrate_data(self, key, info, array, var_attrs, index): if key['calibration'] == 'brightness_temperature': projectable = calibrate_bt(array, var_attrs, index, key['name']) From 4a0bea638f3f30ca55c172d8079509517d3404ad Mon Sep 17 00:00:00 2001 From: David Hoese Date: Fri, 11 Mar 2022 11:05:48 -0600 Subject: [PATCH 07/13] Add documentation to modis_l1b docstring on mask_saturated keyword argument --- satpy/readers/modis_l1b.py | 24 +++++++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/satpy/readers/modis_l1b.py b/satpy/readers/modis_l1b.py index 6dc82f1bda..354a55f7f9 100644 --- a/satpy/readers/modis_l1b.py +++ b/satpy/readers/modis_l1b.py @@ -25,11 +25,33 @@ a pattern similar to the following one: .. parsed-literal:: + M[O/Y]D02[1/H/Q]KM.A[date].[time].[collection].[processing_time].hdf Other patterns where "collection" and/or "proccessing_time" are missing might also work (see the readers yaml file for details). Geolocation files (MOD03) are also supported. - +The IMAPP direct broadcast naming format is also supported with names like: +``a1.12226.1846.1000m.hdf``. + +Saturation Handling +------------------- + +Band 2 of the MODIS sensor is available in 250m, 500m, and 1km resolutions. +The band data may include a special fill value to indicate when the detector +was saturated in the 250m version of the data. When the data is aggregated to +coarser resolutions this saturation fill value is converted to a +"can't aggregate" fill value. By default, Satpy will replace these fill values +with NaN to indicate they are invalid. This is typically undesired when +generating images for the data as they appear as "holes" in bright clouds. +To control this the keyword argument ``mask_saturated`` can be passed and set +to ``False`` to set these two fill values to the maximum valid value. + +.. code-block:: python + + scene = satpy.Scene(filenames=filenames, + reader='modis_l1b', + reader_kwargs={'mask_saturated': False}) + scene.load(['2']) Geolocation files ----------------- From 62c2a1cd89f474dc137c3449e4cf0dd67dea169b Mon Sep 17 00:00:00 2001 From: David Hoese Date: Fri, 11 Mar 2022 12:08:06 -0600 Subject: [PATCH 08/13] Skip uncertainty checks in modis_l1b if not masking saturated pixels --- satpy/readers/modis_l1b.py | 2 ++ satpy/tests/reader_tests/_modis_fixtures.py | 11 ++++++++++- 2 files changed, 12 insertions(+), 1 deletion(-) diff --git a/satpy/readers/modis_l1b.py b/satpy/readers/modis_l1b.py index 354a55f7f9..76ebc62029 100644 --- a/satpy/readers/modis_l1b.py +++ b/satpy/readers/modis_l1b.py @@ -197,6 +197,8 @@ def _mask_invalid_or_fill_saturated(self, array, valid_range): return array def _mask_uncertain_pixels(self, array, uncertainty, band_index): + if not self._mask_saturated: + return array band_uncertainty = from_sds(uncertainty, chunks=CHUNK_SIZE)[band_index, :, :] array = array.where(band_uncertainty < 15) return array diff --git a/satpy/tests/reader_tests/_modis_fixtures.py b/satpy/tests/reader_tests/_modis_fixtures.py index 1fdb596046..ab09635d45 100644 --- a/satpy/tests/reader_tests/_modis_fixtures.py +++ b/satpy/tests/reader_tests/_modis_fixtures.py @@ -85,6 +85,14 @@ def _generate_visible_data(resolution: int, num_bands: int, dtype=np.uint16) -> return data +def _generate_visible_uncertainty_data(shape: tuple) -> np.ndarray: + uncertainty = np.zeros(shape, dtype=np.uint8) + uncertainty[:, -1, -1] = 15 # fill value + uncertainty[:, -1, -2] = 15 # saturated + uncertainty[:, -1, -3] = 15 # can't aggregate + return uncertainty + + def _get_lonlat_variable_info(resolution: int) -> dict: lon_5km, lat_5km = _generate_lonlat_data(resolution) return { @@ -122,6 +130,7 @@ def _get_angles_variable_info(resolution: int) -> dict: def _get_visible_variable_info(var_name: str, resolution: int, bands: list[str]): num_bands = len(bands) data = _generate_visible_data(resolution, len(bands)) + uncertainty = _generate_visible_uncertainty_data(data.shape) dim_factor = RES_TO_REPEAT_FACTOR[resolution] * 2 band_dim_name = f"Band_{resolution}_{num_bands}_RefSB:MODIS_SWATH_Type_L1B" row_dim_name = f'{dim_factor}*nscans:MODIS_SWATH_Type_L1B' @@ -143,7 +152,7 @@ def _get_visible_variable_info(var_name: str, resolution: int, bands: list[str]) }, }, var_name + '_Uncert_Indexes': { - 'data': np.zeros(data.shape, dtype=np.uint8), + 'data': uncertainty, 'type': SDC.UINT8, 'fill_value': 255, 'attrs': { From 445ce79b579004f8dc0269c2081b1858326f01ee Mon Sep 17 00:00:00 2001 From: David Hoese Date: Thu, 31 Mar 2022 10:37:07 -0500 Subject: [PATCH 09/13] Add more information on fill value handling to MODIS L1b reader docstrings --- satpy/readers/modis_l1b.py | 52 +++++++++++++++++++++++--------------- 1 file changed, 32 insertions(+), 20 deletions(-) diff --git a/satpy/readers/modis_l1b.py b/satpy/readers/modis_l1b.py index 76ebc62029..f26b817879 100644 --- a/satpy/readers/modis_l1b.py +++ b/satpy/readers/modis_l1b.py @@ -53,6 +53,12 @@ reader_kwargs={'mask_saturated': False}) scene.load(['2']) +Note that the saturation fill value can appear in other bands (ex. bands 7-19) +in addition to band 2. Also, the "can't aggregate" fill value is a generic +"catch all" for any problems encountered when aggregating high resolution bands +to lower resolutions. Filling this with the max valid value could replace +non-saturated invalid pixels with valid values. + Geolocation files ----------------- @@ -166,26 +172,32 @@ def _get_band_index(self, var_attrs, band_name): return index def _mask_invalid_or_fill_saturated(self, array, valid_range): - """Replace fill values with NaN or replace saturation with max reflectance.""" - # Fill values: - # Data Value Meaning - # 65535 Fill Value (includes reflective band data at night mode - # and completely missing L1A scans) - # 65534 L1A DN is missing within a scan - # 65533 Detector is saturated - # 65532 Cannot compute zero point DN, e.g., SV is saturated - # 65531 Detector is dead (see comments below) - # 65530 RSB dn** below the minimum of the scaling range - # 65529 TEB radiance or RSB dn** exceeds the maximum of the - # scaling range - # 65528 Aggregation algorithm failure - # 65527 Rotation of Earth view Sector from nominal science - # collection position - # 65526 Calibration coefficient b1 could not be computed - # 65525 Subframe is dead - # 65524 Both sides of the PCLW electronics on simultaneously - # 65501 - 65523 (reserved for future use) - # 65500 NAD closed upper limit + """Replace fill values with NaN or replace saturation with max reflectance. + + If the file handler was created with ``mask_saturated`` set to + ``True`` then all invalid/fill values are set to NaN. If ``False`` + then the fill values 65528 and 65533 are set to the maximum valid + value. These values correspond to "can't aggregate" and "saturation". + + Fill values: + + * 65535 Fill Value (includes reflective band data at night mode + and completely missing L1A scans) + * 65534 L1A DN is missing within a scan + * 65533 Detector is saturated + * 65532 Cannot compute zero point DN, e.g., SV is saturated + * 65531 Detector is dead (see comments below) + * 65530 RSB dn** below the minimum of the scaling range + * 65529 TEB radiance or RSB dn exceeds the maximum of the scaling range + * 65528 Aggregation algorithm failure + * 65527 Rotation of Earth view Sector from nominal science collection position + * 65526 Calibration coefficient b1 could not be computed + * 65525 Subframe is dead + * 65524 Both sides of the PCLW electronics on simultaneously + * 65501 - 65523 (reserved for future use) + * 65500 NAD closed upper limit + + """ valid_min = np.float32(valid_range[0]) valid_max = np.float32(valid_range[1]) array = array.where(array >= valid_min) From 17c26dbdb577d2363c4e08cd7d5b167297d162a8 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Tue, 5 Apr 2022 06:55:39 -0500 Subject: [PATCH 10/13] Use super() instead of class reference in satpy/readers/modis_l1b.py Co-authored-by: Martin Raspaud --- satpy/readers/modis_l1b.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/satpy/readers/modis_l1b.py b/satpy/readers/modis_l1b.py index f26b817879..f389e3bae8 100644 --- a/satpy/readers/modis_l1b.py +++ b/satpy/readers/modis_l1b.py @@ -102,7 +102,7 @@ class HDFEOSBandReader(HDFEOSBaseFileReader): def __init__(self, filename, filename_info, filetype_info, mask_saturated=True, **kwargs): """Init the file handler.""" - HDFEOSBaseFileReader.__init__(self, filename, filename_info, filetype_info, **kwargs) + super().__init__(self, filename, filename_info, filetype_info, **kwargs) self._mask_saturated = mask_saturated ds = self.metadata['INVENTORYMETADATA'][ From 856fc9ecfe9c43e1f88b54308dfe9f2dec692cf6 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Tue, 5 Apr 2022 13:25:26 -0500 Subject: [PATCH 11/13] Refactor modis_l1b reader to have less complex methods --- satpy/readers/modis_l1b.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/satpy/readers/modis_l1b.py b/satpy/readers/modis_l1b.py index f389e3bae8..26bc354c23 100644 --- a/satpy/readers/modis_l1b.py +++ b/satpy/readers/modis_l1b.py @@ -113,7 +113,10 @@ def get_dataset(self, key, info): """Read data from file and return the corresponding projectables.""" if self.resolution != key['resolution']: return - subdata, uncertainty, var_attrs, band_index = self._get_band_data_uncertainty_attrs_index(key["name"]) + var_name, band_index = self._get_band_variable_name_and_index(key["name"]) + subdata = self.sd.select(var_name) + var_attrs = subdata.attributes() + uncertainty = self.sd.select(var_name + "_Uncert_Indexes") array = xr.DataArray(from_sds(subdata, chunks=CHUNK_SIZE)[band_index, :, :], dims=['y', 'x']).astype(np.float32) valid_range = var_attrs['valid_range'] @@ -152,18 +155,17 @@ def get_dataset(self, key, info): self._add_satpy_metadata(key, projectable) return projectable - def _get_band_data_uncertainty_attrs_index(self, band_name): - datasets = self.res_to_possible_variable_names[self.resolution] - for dataset in datasets: - subdata = self.sd.select(dataset) + def _get_band_variable_name_and_index(self, band_name): + variable_names = self.res_to_possible_variable_names[self.resolution] + for variable_name in variable_names: + subdata = self.sd.select(variable_name) var_attrs = subdata.attributes() try: band_index = self._get_band_index(var_attrs, band_name) except ValueError: # can't find band in list of bands continue - uncertainty = self.sd.select(dataset + "_Uncert_Indexes") - return subdata, uncertainty, var_attrs, band_index + return variable_name, band_index def _get_band_index(self, var_attrs, band_name): """Get the relative indices of the desired channel.""" From d5a8f6261e5ce63d09774114b5e640cc9078c6c1 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Tue, 5 Apr 2022 13:41:30 -0500 Subject: [PATCH 12/13] Refactor modis_l1b for cleaner masking logic --- satpy/readers/modis_l1b.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/satpy/readers/modis_l1b.py b/satpy/readers/modis_l1b.py index 26bc354c23..366135146d 100644 --- a/satpy/readers/modis_l1b.py +++ b/satpy/readers/modis_l1b.py @@ -120,7 +120,11 @@ def get_dataset(self, key, info): array = xr.DataArray(from_sds(subdata, chunks=CHUNK_SIZE)[band_index, :, :], dims=['y', 'x']).astype(np.float32) valid_range = var_attrs['valid_range'] - array = self._mask_invalid_or_fill_saturated(array, valid_range) + valid_min = np.float32(valid_range[0]) + valid_max = np.float32(valid_range[1]) + if not self._mask_saturated: + array = self._fill_saturated(array, valid_max) + array = self._mask_invalid(array, valid_min, valid_max) array = self._mask_uncertain_pixels(array, uncertainty, band_index) projectable = self._calibrate_data(key, info, array, var_attrs, band_index) @@ -173,8 +177,8 @@ def _get_band_index(self, var_attrs, band_name): index = band_names.index(band_name) return index - def _mask_invalid_or_fill_saturated(self, array, valid_range): - """Replace fill values with NaN or replace saturation with max reflectance. + def _fill_saturated(self, array, valid_max): + """Replace saturation-related values with max reflectance. If the file handler was created with ``mask_saturated`` set to ``True`` then all invalid/fill values are set to NaN. If ``False`` @@ -200,15 +204,11 @@ def _mask_invalid_or_fill_saturated(self, array, valid_range): * 65500 NAD closed upper limit """ - valid_min = np.float32(valid_range[0]) - valid_max = np.float32(valid_range[1]) - array = array.where(array >= valid_min) - if self._mask_saturated: - array = array.where(array <= valid_max) - else: - array = array.where((array != 65533) & (array != 65528), valid_max) - array = array.where(array <= valid_max) - return array + return array.where((array != 65533) & (array != 65528), valid_max) + + def _mask_invalid(self, array, valid_min, valid_max): + """Replace fill values with NaN.""" + return array.where((array >= valid_min) & (array <= valid_max)) def _mask_uncertain_pixels(self, array, uncertainty, band_index): if not self._mask_saturated: From 5e8a088ea473c70765cc91c1ae364128cd2a971d Mon Sep 17 00:00:00 2001 From: David Hoese Date: Fri, 8 Apr 2022 06:55:15 -0500 Subject: [PATCH 13/13] Fix extra "self" argument in modis file handler --- satpy/readers/modis_l1b.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/satpy/readers/modis_l1b.py b/satpy/readers/modis_l1b.py index 366135146d..51d6108f3e 100644 --- a/satpy/readers/modis_l1b.py +++ b/satpy/readers/modis_l1b.py @@ -102,7 +102,7 @@ class HDFEOSBandReader(HDFEOSBaseFileReader): def __init__(self, filename, filename_info, filetype_info, mask_saturated=True, **kwargs): """Init the file handler.""" - super().__init__(self, filename, filename_info, filetype_info, **kwargs) + super().__init__(filename, filename_info, filetype_info, **kwargs) self._mask_saturated = mask_saturated ds = self.metadata['INVENTORYMETADATA'][