From 2ea2cb3b3104cf031d14da66f1ae743f6223a413 Mon Sep 17 00:00:00 2001 From: William Roberts <38170479+wroberts4@users.noreply.github.com> Date: Tue, 30 Apr 2019 16:22:01 -0500 Subject: [PATCH 01/13] Add unit conversion comments and improve logic --- satpy/readers/virr_l1b.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/satpy/readers/virr_l1b.py b/satpy/readers/virr_l1b.py index ddba0545fa..523896b31f 100644 --- a/satpy/readers/virr_l1b.py +++ b/satpy/readers/virr_l1b.py @@ -68,9 +68,11 @@ def get_dataset(self, dataset_id, ds_info): file_key = self.geolocation_prefix + ds_info.get('file_key', dataset_id.name) if self.platform_id == 'FY3B': file_key = file_key.replace('Data/', '') - data = self.get(file_key) - if data is None: - logging.error('File key "{0}" could not be found in file {1}'.format(file_key, self.filename)) + try: + data = self[file_key] + except KeyError: + LOG.error('File key "{0}" could not be found in file {1}'.format(file_key, self.filename)) + raise band_index = ds_info.get('band_index') if band_index is not None: data = data[band_index] @@ -79,6 +81,7 @@ def get_dataset(self, dataset_id, ds_info): if 'E' in dataset_id.name: slope = self[self.l1b_prefix + 'Emissive_Radiance_Scales'].data[:, band_index][:, np.newaxis] intercept = self[self.l1b_prefix + 'Emissive_Radiance_Offsets'].data[:, band_index][:, np.newaxis] + # Converts cm^-1 (wavenumbers) and (mW/m^2)/(str/cm^-1) (radiance data) to SI units m^-1, mW*m^-3*str^-1. radiance_data = rad2temp(self['/attr/' + self.wave_number][band_index] * 100, (data * slope + intercept) * 1e-5) data = xr.DataArray(da.from_array(radiance_data, data.chunks), From bf35d1af2ba3df7772bee50bb6343c0b6c8c6cf1 Mon Sep 17 00:00:00 2001 From: William Roberts <38170479+wroberts4@users.noreply.github.com> Date: Tue, 30 Apr 2019 16:29:03 -0500 Subject: [PATCH 02/13] Add LOG = logging.getLogger(__name__) --- satpy/readers/virr_l1b.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/satpy/readers/virr_l1b.py b/satpy/readers/virr_l1b.py index 523896b31f..b2da02fd0d 100644 --- a/satpy/readers/virr_l1b.py +++ b/satpy/readers/virr_l1b.py @@ -49,6 +49,9 @@ import logging +LOG = logging.getLogger(__name__) + + class VIRR_L1B(HDF5FileHandler): """VIRR_L1B reader.""" From 1d0b3b9039e4d0fa7e3ffa118de3c54202e3d48e Mon Sep 17 00:00:00 2001 From: William Roberts <38170479+wroberts4@users.noreply.github.com> Date: Tue, 30 Apr 2019 16:30:09 -0500 Subject: [PATCH 03/13] Fix styling --- satpy/readers/virr_l1b.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/satpy/readers/virr_l1b.py b/satpy/readers/virr_l1b.py index b2da02fd0d..d00a42d398 100644 --- a/satpy/readers/virr_l1b.py +++ b/satpy/readers/virr_l1b.py @@ -48,7 +48,6 @@ import dask.array as da import logging - LOG = logging.getLogger(__name__) @@ -84,7 +83,8 @@ def get_dataset(self, dataset_id, ds_info): if 'E' in dataset_id.name: slope = self[self.l1b_prefix + 'Emissive_Radiance_Scales'].data[:, band_index][:, np.newaxis] intercept = self[self.l1b_prefix + 'Emissive_Radiance_Offsets'].data[:, band_index][:, np.newaxis] - # Converts cm^-1 (wavenumbers) and (mW/m^2)/(str/cm^-1) (radiance data) to SI units m^-1, mW*m^-3*str^-1. + # Converts cm^-1 (wavenumbers) and (mW/m^2)/(str/cm^-1) (radiance data) + # to SI units m^-1, mW*m^-3*str^-1. radiance_data = rad2temp(self['/attr/' + self.wave_number][band_index] * 100, (data * slope + intercept) * 1e-5) data = xr.DataArray(da.from_array(radiance_data, data.chunks), From 742f92e9e0dc863bc62ace8ff661b359135ffb5d Mon Sep 17 00:00:00 2001 From: William Roberts <38170479+wroberts4@users.noreply.github.com> Date: Wed, 29 May 2019 13:56:03 -0500 Subject: [PATCH 04/13] Add kwargs to innit and improve logging --- satpy/readers/virr_l1b.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/satpy/readers/virr_l1b.py b/satpy/readers/virr_l1b.py index 51c22f391e..67963fb3d4 100644 --- a/satpy/readers/virr_l1b.py +++ b/satpy/readers/virr_l1b.py @@ -53,7 +53,10 @@ class VIRR_L1B(HDF5FileHandler): """VIRR_L1B reader.""" - def __init__(self, filename, filename_info, filetype_info): + def __init__(self, filename, filename_info, filetype_info, **kwargs): + # Prevents an error for kwargs being passed to readers. + if kwargs: + logger.warning("kwargs provided to VIRR_L1B, but kwargs aren't used: '{}'".format(kwargs)) super(VIRR_L1B, self).__init__(filename, filename_info, filetype_info) LOG.debug('day/night flag for {0}: {1}'.format(filename, self['/attr/Day Or Night Flag'])) self.geolocation_prefix = filetype_info['geolocation_prefix'] @@ -69,7 +72,11 @@ def get_dataset(self, dataset_id, ds_info): file_key = self.geolocation_prefix + ds_info.get('file_key', dataset_id.name) if self.platform_id == 'FY3B': file_key = file_key.replace('Data/', '') - data = self[file_key] + try: + data = self[file_key] + except KeyError: + LOG.error('File key "{0}" could not be found in file {1}'.format(file_key, self.filename)) + raise band_index = ds_info.get('band_index') if band_index is not None: data = data[band_index] From 8f8ac475dd784eb8e157288f61b4c2161cced26f Mon Sep 17 00:00:00 2001 From: William Roberts <38170479+wroberts4@users.noreply.github.com> Date: Wed, 29 May 2019 16:43:57 -0500 Subject: [PATCH 05/13] Change logger to LOG --- satpy/readers/virr_l1b.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/satpy/readers/virr_l1b.py b/satpy/readers/virr_l1b.py index 67963fb3d4..2639acd896 100644 --- a/satpy/readers/virr_l1b.py +++ b/satpy/readers/virr_l1b.py @@ -56,7 +56,7 @@ class VIRR_L1B(HDF5FileHandler): def __init__(self, filename, filename_info, filetype_info, **kwargs): # Prevents an error for kwargs being passed to readers. if kwargs: - logger.warning("kwargs provided to VIRR_L1B, but kwargs aren't used: '{}'".format(kwargs)) + LOG.warning("kwargs provided to VIRR_L1B, but kwargs aren't used: '{}'".format(kwargs)) super(VIRR_L1B, self).__init__(filename, filename_info, filetype_info) LOG.debug('day/night flag for {0}: {1}'.format(filename, self['/attr/Day Or Night Flag'])) self.geolocation_prefix = filetype_info['geolocation_prefix'] From e5f50b8d8c7164581391cc073623d5ee5874c039 Mon Sep 17 00:00:00 2001 From: William Roberts <38170479+wroberts4@users.noreply.github.com> Date: Wed, 29 May 2019 16:49:13 -0500 Subject: [PATCH 06/13] Remove redundant warning --- satpy/readers/virr_l1b.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/satpy/readers/virr_l1b.py b/satpy/readers/virr_l1b.py index 2639acd896..9d2ed2313f 100644 --- a/satpy/readers/virr_l1b.py +++ b/satpy/readers/virr_l1b.py @@ -54,9 +54,6 @@ class VIRR_L1B(HDF5FileHandler): """VIRR_L1B reader.""" def __init__(self, filename, filename_info, filetype_info, **kwargs): - # Prevents an error for kwargs being passed to readers. - if kwargs: - LOG.warning("kwargs provided to VIRR_L1B, but kwargs aren't used: '{}'".format(kwargs)) super(VIRR_L1B, self).__init__(filename, filename_info, filetype_info) LOG.debug('day/night flag for {0}: {1}'.format(filename, self['/attr/Day Or Night Flag'])) self.geolocation_prefix = filetype_info['geolocation_prefix'] From 59bf927cbf3a10acfa9b38cdd311b863cf51ef8e Mon Sep 17 00:00:00 2001 From: William Roberts <38170479+wroberts4@users.noreply.github.com> Date: Thu, 30 May 2019 08:33:05 -0500 Subject: [PATCH 07/13] Remove kwargs from init --- satpy/readers/virr_l1b.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/satpy/readers/virr_l1b.py b/satpy/readers/virr_l1b.py index 9d2ed2313f..aefbb0a189 100644 --- a/satpy/readers/virr_l1b.py +++ b/satpy/readers/virr_l1b.py @@ -53,7 +53,7 @@ class VIRR_L1B(HDF5FileHandler): """VIRR_L1B reader.""" - def __init__(self, filename, filename_info, filetype_info, **kwargs): + def __init__(self, filename, filename_info, filetype_info): super(VIRR_L1B, self).__init__(filename, filename_info, filetype_info) LOG.debug('day/night flag for {0}: {1}'.format(filename, self['/attr/Day Or Night Flag'])) self.geolocation_prefix = filetype_info['geolocation_prefix'] From 73e0f64802092786588f94a0a257ea168ceff8e2 Mon Sep 17 00:00:00 2001 From: wroberts Date: Fri, 31 May 2019 15:36:05 -0500 Subject: [PATCH 08/13] Make slope 1 if it's 0 --- satpy/readers/virr_l1b.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/satpy/readers/virr_l1b.py b/satpy/readers/virr_l1b.py index aefbb0a189..f8a4e4cebf 100644 --- a/satpy/readers/virr_l1b.py +++ b/satpy/readers/virr_l1b.py @@ -81,6 +81,9 @@ def get_dataset(self, dataset_id, ds_info): (data <= self[file_key + '/attr/valid_range'][1])) if 'E' in dataset_id.name: slope = self[self.l1b_prefix + 'Emissive_Radiance_Scales'].data[:, band_index][:, np.newaxis] + # 0 slope is invalid. + slope_mask = slope == 0 + slope[slope_mask] = 1 intercept = self[self.l1b_prefix + 'Emissive_Radiance_Offsets'].data[:, band_index][:, np.newaxis] # Converts cm^-1 (wavenumbers) and (mW/m^2)/(str/cm^-1) (radiance data) # to SI units m^-1, mW*m^-3*str^-1. @@ -95,12 +98,20 @@ def get_dataset(self, dataset_id, ds_info): data.data = bt_data elif 'R' in dataset_id.name: slope = self['/attr/RefSB_Cal_Coefficients'][0::2] + # 0 slope is invalid. + slope_mask = slope == 0 + slope[slope_mask] = 1 intercept = self['/attr/RefSB_Cal_Coefficients'][1::2] data = data * slope[band_index] + intercept[band_index] else: + slope = self[file_key + '/attr/Slope'] + # 0 slope is invalid. + slope_mask = slope == 0 + slope[slope_mask] = 1 + intercept = self[file_key + '/attr/Intercept'] data = data.where((data >= self[file_key + '/attr/valid_range'][0]) & (data <= self[file_key + '/attr/valid_range'][1])) - data = self[file_key + '/attr/Intercept'] + self[file_key + '/attr/Slope'] * data + data = slope * data + intercept new_dims = {old: new for old, new in zip(data.dims, ('y', 'x'))} data = data.rename(new_dims) data.attrs.update({'platform_name': self['/attr/Satellite Name'], From fbcb9d5d092dbe53b05bb06f086c42fe919dee8a Mon Sep 17 00:00:00 2001 From: wroberts Date: Fri, 31 May 2019 16:18:33 -0500 Subject: [PATCH 09/13] Fix slope_mask not working on floats --- satpy/readers/virr_l1b.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/satpy/readers/virr_l1b.py b/satpy/readers/virr_l1b.py index f8a4e4cebf..c2a64bc1b4 100644 --- a/satpy/readers/virr_l1b.py +++ b/satpy/readers/virr_l1b.py @@ -83,7 +83,10 @@ def get_dataset(self, dataset_id, ds_info): slope = self[self.l1b_prefix + 'Emissive_Radiance_Scales'].data[:, band_index][:, np.newaxis] # 0 slope is invalid. slope_mask = slope == 0 - slope[slope_mask] = 1 + if isinstance(slope_mask, bool): + slope = 1 if slope_mask else slope + else: + slope[slope_mask] = 1 intercept = self[self.l1b_prefix + 'Emissive_Radiance_Offsets'].data[:, band_index][:, np.newaxis] # Converts cm^-1 (wavenumbers) and (mW/m^2)/(str/cm^-1) (radiance data) # to SI units m^-1, mW*m^-3*str^-1. @@ -100,14 +103,20 @@ def get_dataset(self, dataset_id, ds_info): slope = self['/attr/RefSB_Cal_Coefficients'][0::2] # 0 slope is invalid. slope_mask = slope == 0 - slope[slope_mask] = 1 + if isinstance(slope_mask, bool): + slope = 1 if slope_mask else slope + else: + slope[slope_mask] = 1 intercept = self['/attr/RefSB_Cal_Coefficients'][1::2] data = data * slope[band_index] + intercept[band_index] else: slope = self[file_key + '/attr/Slope'] # 0 slope is invalid. slope_mask = slope == 0 - slope[slope_mask] = 1 + if isinstance(slope_mask, bool): + slope = 1 if slope_mask else slope + else: + slope[slope_mask] = 1 intercept = self[file_key + '/attr/Intercept'] data = data.where((data >= self[file_key + '/attr/valid_range'][0]) & (data <= self[file_key + '/attr/valid_range'][1])) From 58a1eedfaa4bb756183c2464895869d4009631cf Mon Sep 17 00:00:00 2001 From: wroberts Date: Mon, 3 Jun 2019 11:07:15 -0500 Subject: [PATCH 10/13] Make slope correction cleaner --- satpy/readers/virr_l1b.py | 37 +++++++++++-------------------------- 1 file changed, 11 insertions(+), 26 deletions(-) diff --git a/satpy/readers/virr_l1b.py b/satpy/readers/virr_l1b.py index c2a64bc1b4..39e79af61f 100644 --- a/satpy/readers/virr_l1b.py +++ b/satpy/readers/virr_l1b.py @@ -45,6 +45,7 @@ from pyspectral.blackbody import blackbody_wn_rad2temp as rad2temp import numpy as np import dask.array as da +from xarray import DataArray import logging LOG = logging.getLogger(__name__) @@ -69,24 +70,15 @@ def get_dataset(self, dataset_id, ds_info): file_key = self.geolocation_prefix + ds_info.get('file_key', dataset_id.name) if self.platform_id == 'FY3B': file_key = file_key.replace('Data/', '') - try: - data = self[file_key] - except KeyError: - LOG.error('File key "{0}" could not be found in file {1}'.format(file_key, self.filename)) - raise + data = self[file_key] band_index = ds_info.get('band_index') if band_index is not None: data = data[band_index] data = data.where((data >= self[file_key + '/attr/valid_range'][0]) & (data <= self[file_key + '/attr/valid_range'][1])) if 'E' in dataset_id.name: - slope = self[self.l1b_prefix + 'Emissive_Radiance_Scales'].data[:, band_index][:, np.newaxis] - # 0 slope is invalid. - slope_mask = slope == 0 - if isinstance(slope_mask, bool): - slope = 1 if slope_mask else slope - else: - slope[slope_mask] = 1 + slope = self._correct_slope(self[self.l1b_prefix + 'Emissive_Radiance_Scales']. + data[:, band_index][:,np.newaxis]) intercept = self[self.l1b_prefix + 'Emissive_Radiance_Offsets'].data[:, band_index][:, np.newaxis] # Converts cm^-1 (wavenumbers) and (mW/m^2)/(str/cm^-1) (radiance data) # to SI units m^-1, mW*m^-3*str^-1. @@ -100,23 +92,11 @@ def get_dataset(self, dataset_id, ds_info): # new versions of pyspectral can do dask arrays data.data = bt_data elif 'R' in dataset_id.name: - slope = self['/attr/RefSB_Cal_Coefficients'][0::2] - # 0 slope is invalid. - slope_mask = slope == 0 - if isinstance(slope_mask, bool): - slope = 1 if slope_mask else slope - else: - slope[slope_mask] = 1 + slope = self._correct_slope(self['/attr/RefSB_Cal_Coefficients'][0::2]) intercept = self['/attr/RefSB_Cal_Coefficients'][1::2] data = data * slope[band_index] + intercept[band_index] else: - slope = self[file_key + '/attr/Slope'] - # 0 slope is invalid. - slope_mask = slope == 0 - if isinstance(slope_mask, bool): - slope = 1 if slope_mask else slope - else: - slope[slope_mask] = 1 + slope = self._correct_slope(self[file_key + '/attr/Slope']) intercept = self[file_key + '/attr/Intercept'] data = data.where((data >= self[file_key + '/attr/valid_range'][0]) & (data <= self[file_key + '/attr/valid_range'][1])) @@ -135,6 +115,11 @@ def get_dataset(self, dataset_id, ds_info): data.attrs.update({'units': '1'}) return data + def _correct_slope(self, slope): + # 0 slope is invalid. Note: slope can be a scalar or array. + slope = DataArray(slope) + return slope.where(slope != 0, 1) + @property def start_time(self): start_time = self['/attr/Observing Beginning Date'] + 'T' + self['/attr/Observing Beginning Time'] + 'Z' From d093c51033a66b6e4b5eb570a6940a5704990a5b Mon Sep 17 00:00:00 2001 From: wroberts Date: Mon, 3 Jun 2019 11:18:20 -0500 Subject: [PATCH 11/13] Use np.where instead of xarray.DataArray.where --- satpy/readers/virr_l1b.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/satpy/readers/virr_l1b.py b/satpy/readers/virr_l1b.py index 39e79af61f..7e9197bcf7 100644 --- a/satpy/readers/virr_l1b.py +++ b/satpy/readers/virr_l1b.py @@ -45,7 +45,6 @@ from pyspectral.blackbody import blackbody_wn_rad2temp as rad2temp import numpy as np import dask.array as da -from xarray import DataArray import logging LOG = logging.getLogger(__name__) @@ -78,7 +77,7 @@ def get_dataset(self, dataset_id, ds_info): (data <= self[file_key + '/attr/valid_range'][1])) if 'E' in dataset_id.name: slope = self._correct_slope(self[self.l1b_prefix + 'Emissive_Radiance_Scales']. - data[:, band_index][:,np.newaxis]) + data[:, band_index][:, np.newaxis]) intercept = self[self.l1b_prefix + 'Emissive_Radiance_Offsets'].data[:, band_index][:, np.newaxis] # Converts cm^-1 (wavenumbers) and (mW/m^2)/(str/cm^-1) (radiance data) # to SI units m^-1, mW*m^-3*str^-1. @@ -117,8 +116,7 @@ def get_dataset(self, dataset_id, ds_info): def _correct_slope(self, slope): # 0 slope is invalid. Note: slope can be a scalar or array. - slope = DataArray(slope) - return slope.where(slope != 0, 1) + return np.where(slope == 0, 1, slope) @property def start_time(self): From 4292b93513dec7496c3569233d2e64f9dc58fb86 Mon Sep 17 00:00:00 2001 From: wroberts Date: Mon, 3 Jun 2019 13:46:38 -0500 Subject: [PATCH 12/13] Use da.where instead of np.where --- satpy/readers/virr_l1b.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/satpy/readers/virr_l1b.py b/satpy/readers/virr_l1b.py index 7e9197bcf7..7f42caf432 100644 --- a/satpy/readers/virr_l1b.py +++ b/satpy/readers/virr_l1b.py @@ -82,8 +82,7 @@ def get_dataset(self, dataset_id, ds_info): # Converts cm^-1 (wavenumbers) and (mW/m^2)/(str/cm^-1) (radiance data) # to SI units m^-1, mW*m^-3*str^-1. wave_number = self['/attr/' + self.wave_number][band_index] * 100 - bt_data = rad2temp(wave_number, - (data.data * slope + intercept) * 1e-5) + bt_data = rad2temp(wave_number, (data.data * slope + intercept) * 1e-5) if isinstance(bt_data, np.ndarray): # old versions of pyspectral produce numpy arrays data.data = da.from_array(bt_data, chunks=data.data.chunks) @@ -99,7 +98,7 @@ def get_dataset(self, dataset_id, ds_info): intercept = self[file_key + '/attr/Intercept'] data = data.where((data >= self[file_key + '/attr/valid_range'][0]) & (data <= self[file_key + '/attr/valid_range'][1])) - data = slope * data + intercept + data = data * slope + intercept new_dims = {old: new for old, new in zip(data.dims, ('y', 'x'))} data = data.rename(new_dims) data.attrs.update({'platform_name': self['/attr/Satellite Name'], @@ -116,7 +115,7 @@ def get_dataset(self, dataset_id, ds_info): def _correct_slope(self, slope): # 0 slope is invalid. Note: slope can be a scalar or array. - return np.where(slope == 0, 1, slope) + return da.where(slope == 0, 1, slope) @property def start_time(self): From ccdb40166b3ad6879d295c4c608f9c9fd40096e6 Mon Sep 17 00:00:00 2001 From: wroberts Date: Mon, 3 Jun 2019 14:24:14 -0500 Subject: [PATCH 13/13] Add test to check if data returned by get_dataset is a dask array --- satpy/tests/reader_tests/test_virr_l1b.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/satpy/tests/reader_tests/test_virr_l1b.py b/satpy/tests/reader_tests/test_virr_l1b.py index a962f678ee..f4dbd6ce15 100644 --- a/satpy/tests/reader_tests/test_virr_l1b.py +++ b/satpy/tests/reader_tests/test_virr_l1b.py @@ -124,10 +124,12 @@ def _fy3_helper(self, platform_name, reader, Emissive_units): 'E1': 496.542155, 'E2': 297.444511, 'E3': 288.956557, 'solar_zenith_angle': .1, 'satellite_zenith_angle': .1, 'solar_azimuth_angle': .1, 'satellite_azimuth_angle': .1, 'longitude': 10} - datasets = reader.load([band for band, val in band_values.items()]) + datasets = reader.load([band for band in band_values]) for dataset in datasets: + # Object returned by get_dataset. ds = datasets[dataset.name] attributes = ds.attrs + self.assertTrue(isinstance(ds.data, da.Array)) self.assertEqual('VIRR', attributes['sensor']) self.assertEqual(platform_name, attributes['platform_name']) self.assertEqual(datetime.datetime(2018, 12, 25, 21, 41, 47, 90000), attributes['start_time']) @@ -166,7 +168,7 @@ def test_fy3b_file(self): self.assertTrue(FY3B_reader.file_handlers) self._fy3_helper('FY3B', FY3B_reader, 'milliWstts/m^2/cm^(-1)/steradian') - def test_FY3C_file(self): + def test_fy3c_file(self): from satpy.readers import load_reader FY3C_reader = load_reader(self.reader_configs) FY3C_files = FY3C_reader.select_files_from_pathnames(['tf2018359143912.FY3C-L_VIRRX_GEOXX.HDF',