diff --git a/doc/source/readers.rst b/doc/source/readers.rst index d1469d8e13..c1ecd139ec 100644 --- a/doc/source/readers.rst +++ b/doc/source/readers.rst @@ -211,11 +211,10 @@ This is described in the developer guide, see :doc:`dev_guide/custom_reader`. Implemented readers =================== +SEVIRI L1.5 data readers +------------------------ -xRIT-based readers ------------------- - -.. automodule:: satpy.readers.hrit_base +.. automodule:: satpy.readers.seviri_base :noindex: SEVIRI HRIT format reader @@ -224,6 +223,26 @@ SEVIRI HRIT format reader .. automodule:: satpy.readers.seviri_l1b_hrit :noindex: +SEVIRI Native format reader +^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. automodule:: satpy.readers.seviri_l1b_native + :noindex: + +SEVIRI netCDF format reader +^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. automodule:: satpy.readers.seviri_l1b_nc + :noindex: + + +Other xRIT-based readers +------------------------ + +.. automodule:: satpy.readers.hrit_base + :noindex: + + JMA HRIT format reader ^^^^^^^^^^^^^^^^^^^^^^ diff --git a/satpy/readers/seviri_base.py b/satpy/readers/seviri_base.py index 2569050452..eae96f40e1 100644 --- a/satpy/readers/seviri_base.py +++ b/satpy/readers/seviri_base.py @@ -15,11 +15,109 @@ # # You should have received a copy of the GNU General Public License along with # satpy. If not, see . -"""Utilities and helper classes for MSG HRIT/Native data reading. +"""Common functionality for SEVIRI L1.5 data readers. + +Introduction +------------ + +*The Spinning Enhanced Visible and InfraRed Imager (SEVIRI) is the primary +instrument on Meteosat Second Generation (MSG) and has the capacity to observe +the Earth in 12 spectral channels.* + +*Level 1.5 corresponds to image data that has been corrected for all unwanted +radiometric and geometric effects, has been geolocated using a standardised +projection, and has been calibrated and radiance-linearised.* +(From the EUMETSAT documentation) + +Satpy provides the following readers for SEVIRI L1.5 data in different formats: + +- Native: :mod:`satpy.readers.seviri_l1b_native` +- HRIT: :mod:`satpy.readers.seviri_l1b_hrit` +- netCDF: :mod:`satpy.readers.seviri_l1b_nc` + + +Calibration +----------- + +This section describes how to control the calibration of SEVIRI L1.5 data. + + +Calibration to radiance +^^^^^^^^^^^^^^^^^^^^^^^ + +The SEVIRI L1.5 data readers allow for choosing between two file-internal +calibration coefficients to convert counts to radiances: + + - Nominal for all channels (default) + - GSICS where available (IR currently) and nominal for the remaining + channels (VIS & HRV currently) + +In order to change the default behaviour, use the ``reader_kwargs`` keyword +argument upon Scene creation:: + + import satpy + scene = satpy.Scene(filenames, + reader='seviri_l1b_...', + reader_kwargs={'calib_mode': 'GSICS'}) + scene.load(['VIS006', 'IR_108']) + +Furthermore, it is possible to specify external calibration coefficients +for the conversion from counts to radiances. External coefficients take +precedence over internal coefficients, but you can also mix internal and +external coefficients: If external calibration coefficients are specified +for only a subset of channels, the remaining channels will be calibrated +using the chosen file-internal coefficients (nominal or GSICS). + +Calibration coefficients must be specified in [mW m-2 sr-1 (cm-1)-1]. + +In the following example we use external calibration coefficients for the +``VIS006`` & ``IR_108`` channels, and nominal coefficients for the +remaining channels:: + + coefs = {'VIS006': {'gain': 0.0236, 'offset': -1.20}, + 'IR_108': {'gain': 0.2156, 'offset': -10.4}} + scene = satpy.Scene(filenames, + reader='seviri_l1b_...', + reader_kwargs={'ext_calib_coefs': coefs}) + scene.load(['VIS006', 'VIS008', 'IR_108', 'IR_120']) + +In the next example we use external calibration coefficients for the +``VIS006`` & ``IR_108`` channels, GSICS coefficients where available +(other IR channels) and nominal coefficients for the rest:: + + coefs = {'VIS006': {'gain': 0.0236, 'offset': -1.20}, + 'IR_108': {'gain': 0.2156, 'offset': -10.4}} + scene = satpy.Scene(filenames, + reader='seviri_l1b_...', + reader_kwargs={'calib_mode': 'GSICS', + 'ext_calib_coefs': coefs}) + scene.load(['VIS006', 'VIS008', 'IR_108', 'IR_120']) + + +Calibration to reflectance +^^^^^^^^^^^^^^^^^^^^^^^^^^ + +When loading solar channels, the SEVIRI L1.5 data readers apply a correction for +the Sun-Earth distance variation throughout the year - as recommended by +the EUMETSAT document +`Conversion from radiances to reflectances for SEVIRI warm channels`_. +In the unlikely situation that this correction is not required, it can be +removed on a per-channel basis using +:func:`satpy.readers.utils.remove_earthsun_distance_correction`. + References: - MSG Level 1.5 Image Data Format Description - https://www.eumetsat.int/website/wcm/idc/idcplg?IdcService=GET_FILE&dDocName=PDF_TEN_05105_MSG_IMG_DATA&RevisionSelectionMethod=LatestReleased&Rendition=Web + - `MSG Level 1.5 Image Data Format Description`_ + - `Radiometric Calibration of MSG SEVIRI Level 1.5 Image Data in Equivalent Spectral Blackbody Radiance`_ + +.. _Conversion from radiances to reflectances for SEVIRI warm channels: + https://www-cdn.eumetsat.int/files/2020-04/pdf_msg_seviri_rad2refl.pdf + +.. _MSG Level 1.5 Image Data Format Description: + https://www-cdn.eumetsat.int/files/2020-05/pdf_ten_05105_msg_img_data.pdf + +.. _Radiometric Calibration of MSG SEVIRI Level 1.5 Image Data in Equivalent Spectral Blackbody Radiance: + https://www-cdn.eumetsat.int/files/2020-04/pdf_ten_msg_seviri_rad_calib.pdf """ @@ -334,23 +432,29 @@ def images_used(self): mpef_product_header = MpefProductHeader().get() -class SEVIRICalibrationHandler(object): - """Calibration handler for SEVIRI HRIT- and native-formats.""" +class SEVIRICalibrationAlgorithm: + """SEVIRI calibration algorithms.""" + + def __init__(self, platform_id, scan_time): + """Initialize the calibration algorithm.""" + self._platform_id = platform_id + self._scan_time = scan_time - def _convert_to_radiance(self, data, gain, offset): + def convert_to_radiance(self, data, gain, offset): """Calibrate to radiance.""" + data = data.where(data > 0) return (data * gain + offset).clip(0.0, None) def _erads2bt(self, data, channel_name): """Convert effective radiance to brightness temperature.""" - cal_info = CALIB[self.platform_id][channel_name] + cal_info = CALIB[self._platform_id][channel_name] alpha = cal_info["ALPHA"] beta = cal_info["BETA"] - wavenumber = CALIB[self.platform_id][channel_name]["VC"] + wavenumber = CALIB[self._platform_id][channel_name]["VC"] return (self._tl15(data, wavenumber) - beta) / alpha - def _ir_calibrate(self, data, channel_name, cal_type): + def ir_calibrate(self, data, channel_name, cal_type): """Calibrate to brightness temperature.""" if cal_type == 1: # spectral radiances @@ -364,7 +468,7 @@ def _ir_calibrate(self, data, channel_name, cal_type): def _srads2bt(self, data, channel_name): """Convert spectral radiance to brightness temperature.""" a__, b__, c__ = BTFIT[channel_name] - wavenumber = CALIB[self.platform_id][channel_name]["VC"] + wavenumber = CALIB[self._platform_id][channel_name]["VC"] temp = self._tl15(data, wavenumber) return a__ * temp * temp + b__ * temp + c__ @@ -374,14 +478,95 @@ def _tl15(self, data, wavenumber): return ((C2 * wavenumber) / np.log((1.0 / data) * C1 * wavenumber ** 3 + 1.0)) - def _vis_calibrate(self, data, solar_irradiance): + def vis_calibrate(self, data, solar_irradiance): """Calibrate to reflectance. This uses the method described in Conversion from radiances to reflectances for SEVIRI warm channels: https://tinyurl.com/y67zhphm """ reflectance = np.pi * data * 100.0 / solar_irradiance - return apply_earthsun_distance_correction(reflectance, self.start_time) + return apply_earthsun_distance_correction(reflectance, self._scan_time) + + +class SEVIRICalibrationHandler: + """Calibration handler for SEVIRI HRIT-, native- and netCDF-formats. + + Handles selection of calibration coefficients and calls the appropriate + calibration algorithm. + """ + + def __init__(self, platform_id, channel_name, coefs, calib_mode, scan_time): + """Initialize the calibration handler.""" + self._platform_id = platform_id + self._channel_name = channel_name + self._coefs = coefs + self._calib_mode = calib_mode.upper() + self._scan_time = scan_time + self._algo = SEVIRICalibrationAlgorithm( + platform_id=self._platform_id, + scan_time=self._scan_time + ) + + valid_modes = ('NOMINAL', 'GSICS') + if self._calib_mode not in valid_modes: + raise ValueError( + 'Invalid calibration mode: {}. Choose one of {}'.format( + self._calib_mode, valid_modes) + ) + + def calibrate(self, data, calibration): + """Calibrate the given data.""" + if calibration == 'counts': + res = data + elif calibration in ['radiance', 'reflectance', + 'brightness_temperature']: + gain, offset = self.get_gain_offset() + res = self._algo.convert_to_radiance( + data.astype(np.float32), gain, offset + ) + else: + raise ValueError( + 'Invalid calibration {} for channel {}'.format( + calibration, self._channel_name + ) + ) + + if calibration == 'reflectance': + solar_irradiance = CALIB[self._platform_id][self._channel_name]["F"] + res = self._algo.vis_calibrate(res, solar_irradiance) + elif calibration == 'brightness_temperature': + res = self._algo.ir_calibrate( + res, self._channel_name, self._coefs['radiance_type'] + ) + + return res + + def get_gain_offset(self): + """Get gain & offset for calibration from counts to radiance. + + Choices for internal coefficients are nominal or GSICS. If no + GSICS coefficients are available for a certain channel, fall back to + nominal coefficients. External coefficients take precedence over + internal coefficients. + """ + coefs = self._coefs['coefs'] + + # Select internal coefficients for the given calibration mode + internal_gain = coefs['NOMINAL']['gain'] + internal_offset = coefs['NOMINAL']['offset'] + if self._calib_mode == 'GSICS': + gsics_gain = coefs['GSICS']['gain'] + gsics_offset = coefs['GSICS']['offset'] * gsics_gain + if gsics_gain != 0 and gsics_offset != 0: + # If no GSICS coefficients are available for a certain channel, + # they are set to zero in the file. + internal_gain = gsics_gain + internal_offset = gsics_offset + + # Override with external coefficients, if any. + gain = coefs['EXTERNAL'].get('gain', internal_gain) + offset = coefs['EXTERNAL'].get('offset', internal_offset) + return gain, offset def chebyshev(coefs, time, domain): @@ -432,6 +617,24 @@ def calculate_area_extent(area_dict): return tuple(aex) +def create_coef_dict(coefs_nominal, coefs_gsics, radiance_type, ext_coefs): + """Create coefficient dictionary expected by calibration class.""" + return { + 'coefs': { + 'NOMINAL': { + 'gain': coefs_nominal[0], + 'offset': coefs_nominal[1], + }, + 'GSICS': { + 'gain': coefs_gsics[0], + 'offset': coefs_gsics[1] + }, + 'EXTERNAL': ext_coefs + }, + 'radiance_type': radiance_type + } + + def get_service_mode(ssp_lon): """Get information about whether we're dealing with data from the FES, RSS or IODC service mode.""" seviri_service_modes = {'0.0': {'name': 'FES', 'desc': 'Full Earth scanning service'}, diff --git a/satpy/readers/seviri_l1b_hrit.py b/satpy/readers/seviri_l1b_hrit.py index 6d6c23fc39..b35438252a 100644 --- a/satpy/readers/seviri_l1b_hrit.py +++ b/satpy/readers/seviri_l1b_hrit.py @@ -21,7 +21,7 @@ ------------ The ``seviri_l1b_hrit`` reader reads and calibrates MSG-SEVIRI L1.5 image data in HRIT format. The format is explained -in the `MSG Level 1.5 Image Format Description`_. The files are usually named as +in the `MSG Level 1.5 Image Data Format Description`_. The files are usually named as follows: .. code-block:: none @@ -104,7 +104,7 @@ * The ``orbital_parameters`` attribute provides the nominal and actual satellite position, as well as the projection centre. * You can choose between nominal and GSICS calibration coefficients or even specify your own coefficients, see - :class:`HRITMSGFileHandler`. + :class:`satpy.readers.seviri_base`. * The ``raw_metadata`` attribute provides raw metadata from the prologue, epilogue and segment header. By default, arrays with more than 100 elements are excluded in order to limit memory usage. This threshold can be adjusted, see :class:`HRITMSGFileHandler`. @@ -119,34 +119,6 @@ scn['IR_108']['y'] = mi scn['IR_108'].sel(time=np.datetime64('2019-03-01T12:06:13.052000000')) -Notes: - When loading solar channels, this reader applies a correction for the - Sun-Earth distance variation throughout the year - as recommended by - the EUMETSAT document: - 'Conversion from radiances to reflectances for SEVIRI warm channels' - In the unlikely situation that this correction is not required, it can be - removed on a per-channel basis using the - satpy.readers.utils.remove_earthsun_distance_correction(channel, utc_time) - function. - - -References: - - `MSG Level 1.5 Image Format Description`_ - - `Radiometric Calibration of MSG SEVIRI Level 1.5 Image Data in Equivalent Spectral Blackbody Radiance`_ - - `Conversion from radiances to reflectances for SEVIRI warm channels`_ - - -.. _MSG Level 1.5 Image Format Description: http://www.eumetsat.int/website/wcm/idc/idcplg?IdcService=GET_FILE&dDocName= - PDF_TEN_05105_MSG_IMG_DATA&RevisionSelectionMethod=LatestReleased&Rendition=Web - -.. _Radiometric Calibration of MSG SEVIRI Level 1.5 Image Data in Equivalent Spectral Blackbody Radiance: - https://www.eumetsat.int/website/wcm/idc/idcplg?IdcService=GET_FILE&dDocName=PDF_TEN_MSG_SEVIRI_RAD_CALIB& - RevisionSelectionMethod=LatestReleased&Rendition=Web - -.. _Conversion from radiances to reflectances for SEVIRI warm channels: - https://www.eumetsat.int/website/wcm/idc/idcplg?IdcService=GET_FILE&dDocName=PDF_MSG_SEVIRI_RAD2REFL& - RevisionSelectionMethod=LatestReleased&Rendition=Web - """ from __future__ import division @@ -160,16 +132,18 @@ import pyproj import xarray as xr -import satpy.readers.utils as utils from pyresample import geometry +from satpy import CHUNK_SIZE +import satpy.readers.utils as utils from satpy.readers.eum_base import recarray2dict, time_cds_short from satpy.readers.hrit_base import (HRITFileHandler, ancillary_text, annotation_header, base_hdr_map, image_data_function) -from satpy.readers.seviri_base import (CALIB, CHANNEL_NAMES, SATNUM, - VIS_CHANNELS, SEVIRICalibrationHandler, + +from satpy.readers.seviri_base import (CHANNEL_NAMES, SATNUM, + SEVIRICalibrationHandler, chebyshev, get_cds_time, HRV_NUM_COLUMNS, - pad_data_horizontally) + pad_data_horizontally, create_coef_dict) from satpy.readers.seviri_l1b_native_hdr import (hrit_epilogue, hrit_prologue, impf_configuration) from satpy.readers._geos_area import get_area_extent, get_area_definition @@ -433,55 +407,12 @@ def reduce(self, max_size): return self._reduce(self.epilogue, max_size=max_size) -class HRITMSGFileHandler(HRITFileHandler, SEVIRICalibrationHandler): +class HRITMSGFileHandler(HRITFileHandler): """SEVIRI HRIT format reader. **Calibration** - It is possible to choose between two file-internal calibration coefficients for the conversion - from counts to radiances: - - - Nominal for all channels (default) - - GSICS for IR channels and nominal for VIS channels - - In order to change the default behaviour, use the ``reader_kwargs`` upon Scene creation:: - - import satpy - import glob - - filenames = glob.glob('H-000-MSG3*') - scene = satpy.Scene(filenames, - reader='seviri_l1b_hrit', - reader_kwargs={'calib_mode': 'GSICS'}) - scene.load(['VIS006', 'IR_108']) - - Furthermore, it is possible to specify external calibration coefficients for the conversion from - counts to radiances. They must be specified in [mW m-2 sr-1 (cm-1)-1]. External coefficients - take precedence over internal coefficients. If external calibration coefficients are specified - for only a subset of channels, the remaining channels will be calibrated using the chosen - file-internal coefficients (nominal or GSICS). - - In the following example we use external calibration coefficients for the ``VIS006`` & - ``IR_108`` channels, and nominal coefficients for the remaining channels:: - - coefs = {'VIS006': {'gain': 0.0236, 'offset': -1.20}, - 'IR_108': {'gain': 0.2156, 'offset': -10.4}} - scene = satpy.Scene(filenames, - reader='seviri_l1b_hrit', - reader_kwargs={'ext_calib_coefs': coefs}) - scene.load(['VIS006', 'VIS008', 'IR_108', 'IR_120']) - - In the next example we use we use external calibration coefficients for the ``VIS006`` & - ``IR_108`` channels, nominal coefficients for the remaining VIS channels and GSICS coefficients - for the remaining IR channels:: - - coefs = {'VIS006': {'gain': 0.0236, 'offset': -1.20}, - 'IR_108': {'gain': 0.2156, 'offset': -10.4}} - scene = satpy.Scene(filenames, - reader='seviri_l1b_hrit', - reader_kwargs={'calib_mode': 'GSICS', - 'ext_calib_coefs': coefs}) - scene.load(['VIS006', 'VIS008', 'IR_108', 'IR_120']) + See :mod:`satpy.readers.seviri_base`. **Raw Metadata** @@ -520,14 +451,10 @@ def __init__(self, filename, filename_info, filetype_info, self.prologue = prologue.prologue self.epilogue = epilogue.epilogue self._filename_info = filename_info - self.ext_calib_coefs = ext_calib_coefs if ext_calib_coefs is not None else {} self.mda_max_array_size = mda_max_array_size self.fill_hrv = fill_hrv - calib_mode_choices = ('NOMINAL', 'GSICS') - if calib_mode.upper() not in calib_mode_choices: - raise ValueError('Invalid calibration mode: {}. Choose one of {}'.format( - calib_mode, calib_mode_choices)) - self.calib_mode = calib_mode.upper() + self.calib_mode = calib_mode + self.ext_calib_coefs = ext_calib_coefs or {} self._get_header() @@ -751,49 +678,30 @@ def pad_hrv_data(self, res): def calibrate(self, data, calibration): """Calibrate the data.""" tic = datetime.now() - channel_name = self.channel_name - - if calibration == 'counts': - res = data - elif calibration in ['radiance', 'reflectance', 'brightness_temperature']: - # Choose calibration coefficients - # a) Internal: Nominal or GSICS? - band_idx = self.mda['spectral_channel_id'] - 1 - if self.calib_mode != 'GSICS' or self.channel_name in VIS_CHANNELS: - # you cant apply GSICS values to the VIS channels - coefs = self.prologue["RadiometricProcessing"]["Level15ImageCalibration"] - int_gain = coefs['CalSlope'][band_idx] - int_offset = coefs['CalOffset'][band_idx] - else: - coefs = self.prologue["RadiometricProcessing"]['MPEFCalFeedback'] - int_gain = coefs['GSICSCalCoeff'][band_idx] - int_offset = coefs['GSICSOffsetCount'][band_idx] * int_gain - - # b) Internal or external? External takes precedence. - gain = self.ext_calib_coefs.get(self.channel_name, {}).get('gain', int_gain) - offset = self.ext_calib_coefs.get(self.channel_name, {}).get('offset', int_offset) - - # Convert to radiance - data = data.where(data > 0) - res = self._convert_to_radiance(data.astype(np.float32), gain, offset) - line_mask = self.mda['image_segment_line_quality']['line_validity'] >= 2 - line_mask &= self.mda['image_segment_line_quality']['line_validity'] <= 3 - line_mask &= self.mda['image_segment_line_quality']['line_radiometric_quality'] == 4 - line_mask &= self.mda['image_segment_line_quality']['line_geometric_quality'] == 4 - res *= np.choose(line_mask, [1, np.nan])[:, np.newaxis].astype(np.float32) - - if calibration == 'reflectance': - solar_irradiance = CALIB[self.platform_id][channel_name]["F"] - res = self._vis_calibrate(res, solar_irradiance) - - elif calibration == 'brightness_temperature': - cal_type = self.prologue['ImageDescription'][ - 'Level15ImageProduction']['PlannedChanProcessing'][self.mda['spectral_channel_id']] - res = self._ir_calibrate(res, channel_name, cal_type) - + calib = SEVIRICalibrationHandler( + platform_id=self.platform_id, + channel_name=self.channel_name, + coefs=self._get_calib_coefs(self.channel_name), + calib_mode=self.calib_mode, + scan_time=self.start_time + ) + res = calib.calibrate(data, calibration) + if calibration in ['radiance', 'reflectance', 'brightness_temperature']: + res = self._mask_bad_quality(res) logger.debug("Calibration time " + str(datetime.now() - tic)) return res + def _mask_bad_quality(self, data): + """Mask scanlines with bad quality.""" + # Based on missing (2) or corrupted (3) data + line_mask = self.mda['image_segment_line_quality']['line_validity'] >= 2 + line_mask &= self.mda['image_segment_line_quality']['line_validity'] <= 3 + # Do not use (4) + line_mask &= self.mda['image_segment_line_quality']['line_radiometric_quality'] == 4 + line_mask &= self.mda['image_segment_line_quality']['line_geometric_quality'] == 4 + data *= np.choose(line_mask, [1, np.nan])[:, np.newaxis].astype(np.float32) + return data + def _get_raw_mda(self): """Compile raw metadata to be included in the dataset attributes.""" # Metadata from segment header (excluding items which vary among the different segments) @@ -811,3 +719,39 @@ def _get_timestamps(self): """Read scanline timestamps from the segment header.""" tline = self.mda['image_segment_line_quality']['line_mean_acquisition'] return get_cds_time(days=tline['days'], msecs=tline['milliseconds']) + + def _get_calib_coefs(self, channel_name): + """Get coefficients for calibration from counts to radiance.""" + band_idx = self.mda['spectral_channel_id'] - 1 + coefs_nominal = self.prologue["RadiometricProcessing"][ + "Level15ImageCalibration"] + coefs_gsics = self.prologue["RadiometricProcessing"]['MPEFCalFeedback'] + radiance_types = self.prologue['ImageDescription'][ + 'Level15ImageProduction']['PlannedChanProcessing'] + return create_coef_dict( + coefs_nominal=( + coefs_nominal['CalSlope'][band_idx], + coefs_nominal['CalOffset'][band_idx] + ), + coefs_gsics=( + coefs_gsics['GSICSCalCoeff'][band_idx], + coefs_gsics['GSICSOffsetCount'][band_idx] + ), + ext_coefs=self.ext_calib_coefs.get(channel_name, {}), + radiance_type=radiance_types[band_idx] + ) + + +def pad_data(data, final_size, east_bound, west_bound): + """Pad the data given east and west bounds and the desired size.""" + nlines = final_size[0] + if west_bound - east_bound != data.shape[1] - 1: + raise IndexError('East and west bounds do not match data shape') + padding_east = da.zeros((nlines, east_bound - 1), + dtype=data.dtype, chunks=CHUNK_SIZE) + padding_west = da.zeros((nlines, (final_size[1] - west_bound)), + dtype=data.dtype, chunks=CHUNK_SIZE) + if np.issubdtype(data.dtype, np.floating): + padding_east = padding_east * np.nan + padding_west = padding_west * np.nan + return np.hstack((padding_east, data, padding_west)) diff --git a/satpy/readers/seviri_l1b_native.py b/satpy/readers/seviri_l1b_native.py index 85a964c8be..ad5ee34b89 100644 --- a/satpy/readers/seviri_l1b_native.py +++ b/satpy/readers/seviri_l1b_native.py @@ -17,27 +17,11 @@ # satpy. If not, see . """SEVIRI native format reader. -Notes: - When loading solar channels, this reader applies a correction for the - Sun-Earth distance variation throughout the year - as recommended by - the EUMETSAT document: - 'Conversion from radiances to reflectances for SEVIRI warm channels' - In the unlikely situation that this correction is not required, it can be - removed on a per-channel basis using the - satpy.readers.utils.remove_earthsun_distance_correction(channel, utc_time) - function. - References: - `MSG Level 1.5 Native Format File Definition`_ - - `MSG Level 1.5 Image Data Format Description`_ - - `Conversion from radiances to reflectances for SEVIRI warm channels`_ -.. _MSG Level 1.5 Native Format File Definition +.. _MSG Level 1.5 Native Format File Definition: https://www-cdn.eumetsat.int/files/2020-04/pdf_fg15_msg-native-format-15.pdf -.. _MSG Level 1.5 Image Data Format Description - https://www-cdn.eumetsat.int/files/2020-05/pdf_ten_05105_msg_img_data.pdf -.. _Conversion from radiances to reflectances for SEVIRI warm channels: - https://www-cdn.eumetsat.int/files/2020-04/pdf_msg_seviri_rad2refl.pdf """ @@ -54,11 +38,12 @@ from satpy.readers.file_handlers import BaseFileHandler from satpy.readers.eum_base import recarray2dict -from satpy.readers.seviri_base import (SEVIRICalibrationHandler, - CHANNEL_NAMES, CALIB, SATNUM, - dec10216, VISIR_NUM_COLUMNS, - VISIR_NUM_LINES, HRV_NUM_COLUMNS, HRV_NUM_LINES, - VIS_CHANNELS, get_service_mode, pad_data_horizontally, pad_data_vertically) +from satpy.readers.seviri_base import ( + SEVIRICalibrationHandler, CHANNEL_NAMES, SATNUM, dec10216, + VISIR_NUM_COLUMNS, VISIR_NUM_LINES, HRV_NUM_COLUMNS, HRV_NUM_LINES, + create_coef_dict, get_service_mode, pad_data_horizontally, + pad_data_vertically +) from satpy.readers.seviri_l1b_native_hdr import (GSDTRecords, native_header, native_trailer) from satpy.readers._geos_area import get_area_definition @@ -66,12 +51,12 @@ logger = logging.getLogger('native_msg') -class NativeMSGFileHandler(BaseFileHandler, SEVIRICalibrationHandler): +class NativeMSGFileHandler(BaseFileHandler): """SEVIRI native format reader. - The Level1.5 Image data calibration method can be changed by adding the - required mode to the Scene object instantiation kwargs eg - kwargs = {"calib_mode": "gsics",} + **Calibration** + + See :mod:`satpy.readers.seviri_base`. **Padding channel data to full disk** @@ -85,13 +70,15 @@ class NativeMSGFileHandler(BaseFileHandler, SEVIRICalibrationHandler): reader_kwargs={'fill_disk': False}) """ - def __init__(self, filename, filename_info, filetype_info, calib_mode='nominal', fill_disk=False): + def __init__(self, filename, filename_info, filetype_info, + calib_mode='nominal', fill_disk=False, ext_calib_coefs=None): """Initialize the reader.""" super(NativeMSGFileHandler, self).__init__(filename, filename_info, filetype_info) self.platform_name = None self.calib_mode = calib_mode + self.ext_calib_coefs = ext_calib_coefs or {} self.fill_disk = fill_disk # Declare required variables. @@ -471,51 +458,43 @@ def _get_hrv_channel(self): def calibrate(self, data, dataset_id): """Calibrate the data.""" tic = datetime.now() + channel_name = dataset_id['name'] + calib = SEVIRICalibrationHandler( + platform_id=self.platform_id, + channel_name=channel_name, + coefs=self._get_calib_coefs(channel_name), + calib_mode=self.calib_mode, + scan_time=self.start_time + ) + res = calib.calibrate(data, dataset_id['calibration']) + logger.debug("Calibration time " + str(datetime.now() - tic)) + return res - data15hdr = self.header['15_DATA_HEADER'] - calibration = dataset_id['calibration'] - channel = dataset_id['name'] - + def _get_calib_coefs(self, channel_name): + """Get coefficients for calibration from counts to radiance.""" # even though all the channels may not be present in the file, # the header does have calibration coefficients for all the channels # hence, this channel index needs to refer to full channel list - i = list(CHANNEL_NAMES.values()).index(channel) - - if calibration == 'counts': - return data - - if calibration in ['radiance', 'reflectance', 'brightness_temperature']: - # determine the required calibration coefficients to use - # for the Level 1.5 Header - if (self.calib_mode.upper() != 'GSICS' and self.calib_mode.upper() != 'NOMINAL'): - raise NotImplementedError( - 'Unknown Calibration mode : Please check') - - # NB GSICS doesn't have calibration coeffs for VIS channels - if (self.calib_mode.upper() != 'GSICS' or channel in VIS_CHANNELS): - coeffs = data15hdr[ - 'RadiometricProcessing']['Level15ImageCalibration'] - gain = coeffs['CalSlope'][i] - offset = coeffs['CalOffset'][i] - else: - coeffs = data15hdr[ - 'RadiometricProcessing']['MPEFCalFeedback'] - gain = coeffs['GSICSCalCoeff'][i] - offset = coeffs['GSICSOffsetCount'][i] - offset = offset * gain - res = self._convert_to_radiance(data, gain, offset) - - if calibration == 'reflectance': - solar_irradiance = CALIB[self.platform_id][channel]["F"] - res = self._vis_calibrate(res, solar_irradiance) - - elif calibration == 'brightness_temperature': - cal_type = data15hdr['ImageDescription'][ - 'Level15ImageProduction']['PlannedChanProcessing'][i] - res = self._ir_calibrate(res, channel, cal_type) - - logger.debug("Calibration time " + str(datetime.now() - tic)) - return res + band_idx = list(CHANNEL_NAMES.values()).index(channel_name) + + coefs_nominal = self.header['15_DATA_HEADER'][ + 'RadiometricProcessing']['Level15ImageCalibration'] + coefs_gsics = self.header['15_DATA_HEADER'][ + 'RadiometricProcessing']['MPEFCalFeedback'] + radiance_types = self.header['15_DATA_HEADER']['ImageDescription'][ + 'Level15ImageProduction']['PlannedChanProcessing'] + return create_coef_dict( + coefs_nominal=( + coefs_nominal['CalSlope'][band_idx], + coefs_nominal['CalOffset'][band_idx] + ), + coefs_gsics=( + coefs_gsics['GSICSCalCoeff'][band_idx], + coefs_gsics['GSICSOffsetCount'][band_idx] + ), + ext_coefs=self.ext_calib_coefs.get(channel_name, {}), + radiance_type=radiance_types[band_idx] + ) class ImageBoundaries: diff --git a/satpy/readers/seviri_l1b_nc.py b/satpy/readers/seviri_l1b_nc.py index ce76e082c2..de60610739 100644 --- a/satpy/readers/seviri_l1b_nc.py +++ b/satpy/readers/seviri_l1b_nc.py @@ -15,35 +15,11 @@ # # You should have received a copy of the GNU General Public License along with # satpy. If not, see . -"""SEVIRI netcdf format reader. - -Notes: - When loading solar channels, this reader applies a correction for the - Sun-Earth distance variation throughout the year - as recommended by - the EUMETSAT document: - 'Conversion from radiances to reflectances for SEVIRI warm channels' - In the unlikely situation that this correction is not required, it can be - removed on a per-channel basis using the - satpy.readers.utils.remove_earthsun_distance_correction(channel, utc_time) - function. - -References: - - - `MSG Level 1.5 Image Data Format Description`_ - - `Conversion from radiances to reflectances for SEVIRI warm channels`_ - -.. _MSG Level 1.5 Image Data Format Description: - https://www.eumetsat.int/website/wcm/idc/idcplg?IdcService=GET_FILE&dDocName=PDF_TEN_05105_MSG_IMG_DATA& - RevisionSelectionMethod=LatestReleased&Rendition=Web - -.. _Conversion from radiances to reflectances for SEVIRI warm channels: - https://www.eumetsat.int/website/wcm/idc/idcplg?IdcService=GET_FILE&dDocName=PDF_MSG_SEVIRI_RAD2REFL& - RevisionSelectionMethod=LatestReleased&Rendition=Web -""" +"""SEVIRI netcdf format reader.""" from satpy.readers.file_handlers import BaseFileHandler from satpy.readers.seviri_base import (SEVIRICalibrationHandler, - CHANNEL_NAMES, CALIB, SATNUM) + CHANNEL_NAMES, SATNUM) import xarray as xr from satpy.readers._geos_area import get_area_definition @@ -52,12 +28,22 @@ import datetime -class NCSEVIRIFileHandler(BaseFileHandler, SEVIRICalibrationHandler): - """File handler for NC seviri files.""" +class NCSEVIRIFileHandler(BaseFileHandler): + """File handler for NC seviri files. + + **Calibration** + + See :mod:`satpy.readers.seviri_base`. Note that there is only one set of + calibration coefficients available in the netCDF files and therefore there + is no `calib_mode` argument. + + """ - def __init__(self, filename, filename_info, filetype_info): + def __init__(self, filename, filename_info, filetype_info, + ext_calib_coefs=None): """Init the file handler.""" super(NCSEVIRIFileHandler, self).__init__(filename, filename_info, filetype_info) + self.ext_calib_coefs = ext_calib_coefs or {} self.nc = None self.mda = {} self.reference = datetime.datetime(1958, 1, 1) @@ -113,7 +99,6 @@ def _read_file(self): def get_dataset(self, dataset_id, dataset_info): """Get the dataset.""" channel = dataset_id['name'] - i = list(CHANNEL_NAMES.values()).index(channel) if (channel == 'HRV'): self.nc = self.nc.rename({'num_columns_hrv': 'x', 'num_rows_hrv': 'y'}) @@ -127,33 +112,16 @@ def get_dataset(self, dataset_id, dataset_info): pass dataset = self.nc[dataset_info['nc_key']] - dataset.attrs.update(dataset_info) - - # Calibrate the data as needed - # MPEF MSG calibration coeffiencts (gain and count) - offset = dataset.attrs['add_offset'].astype('float32') - gain = dataset.attrs['scale_factor'].astype('float32') self.platform_id = int(self.nc.attrs['satellite_id']) - cal_type = self.nc['planned_chan_processing'].values[i] # Correct for the scan line order dataset = dataset.sel(y=slice(None, None, -1)) - if dataset_id['calibration'] == 'counts': - dataset.attrs['_FillValue'] = 0 - - if dataset_id['calibration'] in ['radiance', 'reflectance', 'brightness_temperature']: - dataset = dataset.where(dataset != 0).astype('float32') - dataset = self._convert_to_radiance(dataset, gain, offset) - - if dataset_id['calibration'] == 'reflectance': - solar_irradiance = CALIB[int(self.platform_id)][channel]["F"] - dataset = self._vis_calibrate(dataset, solar_irradiance) - - elif dataset_id['calibration'] == 'brightness_temperature': - dataset = self._ir_calibrate(dataset, channel, cal_type) + # Calibrate the data as needed + dataset = self.calibrate(dataset, dataset_id) + # Update dataset attributes dataset.attrs.update(self.nc[dataset_info['nc_key']].attrs) dataset.attrs.update(dataset_info) dataset.attrs['platform_name'] = "Meteosat-" + SATNUM[self.platform_id] @@ -170,6 +138,41 @@ def get_dataset(self, dataset_id, dataset_info): return dataset + def calibrate(self, dataset, dataset_id): + """Calibrate the data.""" + channel = dataset_id['name'] + calibration = dataset_id['calibration'] + + if dataset_id['calibration'] == 'counts': + dataset.attrs['_FillValue'] = 0 + + calib = SEVIRICalibrationHandler( + platform_id=int(self.platform_id), + channel_name=channel, + coefs=self._get_calib_coefs(dataset, channel), + calib_mode='NOMINAL', + scan_time=self.start_time + ) + + return calib.calibrate(dataset, calibration) + + def _get_calib_coefs(self, dataset, channel): + """Get coefficients for calibration from counts to radiance.""" + band_idx = list(CHANNEL_NAMES.values()).index(channel) + offset = dataset.attrs['add_offset'].astype('float32') + gain = dataset.attrs['scale_factor'].astype('float32') + # Only one calibration available here + return { + 'coefs': { + 'NOMINAL': { + 'gain': gain, + 'offset': offset + }, + 'EXTERNAL': self.ext_calib_coefs.get(channel, {}) + }, + 'radiance_type': self.nc['planned_chan_processing'].values[band_idx] + } + def get_area_def(self, dataset_id): """Get the area def.""" pdict = {} diff --git a/satpy/tests/reader_tests/test_seviri_base.py b/satpy/tests/reader_tests/test_seviri_base.py index 90602944fb..2d49fb1568 100644 --- a/satpy/tests/reader_tests/test_seviri_base.py +++ b/satpy/tests/reader_tests/test_seviri_base.py @@ -18,6 +18,7 @@ """Test the MSG common (native and hrit format) functionionalities.""" import unittest + import numpy as np import xarray as xr import dask.array as da diff --git a/satpy/tests/reader_tests/test_seviri_l1b_calibration.py b/satpy/tests/reader_tests/test_seviri_l1b_calibration.py index c1a0dc8bed..1f9a01f6d0 100644 --- a/satpy/tests/reader_tests/test_seviri_l1b_calibration.py +++ b/satpy/tests/reader_tests/test_seviri_l1b_calibration.py @@ -17,28 +17,38 @@ # satpy. If not, see . """Unittesting the native msg reader.""" +from datetime import datetime import unittest + import numpy as np +import pytest import xarray as xr -from datetime import datetime -from satpy.readers.seviri_base import SEVIRICalibrationHandler - -COUNTS_INPUT = np.array([[377., 377., 377., 376., 375.], - [376., 375., 376., 374., 374.], - [374., 373., 373., 374., 374.], - [347., 345., 345., 348., 347.], - [306., 306., 307., 307., 308.]], dtype=np.float32) - -RADIANCES_OUTPUT = np.array([[66.84162903, 66.84162903, 66.84162903, 66.63659668, - 66.4315567], - [66.63659668, 66.4315567, 66.63659668, 66.22652435, - 66.22652435], - [66.22652435, 66.02148438, 66.02148438, 66.22652435, - 66.22652435], - [60.69055939, 60.28048706, 60.28048706, 60.89559937, - 60.69055939], - [52.28409576, 52.28409576, 52.48912811, 52.48912811, - 52.69416809]], dtype=np.float32) + +from satpy.readers.seviri_base import ( + SEVIRICalibrationHandler, SEVIRICalibrationAlgorithm +) + + +COUNTS_INPUT = xr.DataArray( + np.array([[377., 377., 377., 376., 375.], + [376., 375., 376., 374., 374.], + [374., 373., 373., 374., 374.], + [347., 345., 345., 348., 347.], + [306., 306., 307., 307., 308.]], dtype=np.float32) +) + +RADIANCES_OUTPUT = xr.DataArray( + np.array([[66.84162903, 66.84162903, 66.84162903, 66.63659668, + 66.4315567], + [66.63659668, 66.4315567, 66.63659668, 66.22652435, + 66.22652435], + [66.22652435, 66.02148438, 66.02148438, 66.22652435, + 66.22652435], + [60.69055939, 60.28048706, 60.28048706, 60.89559937, + 60.69055939], + [52.28409576, 52.28409576, 52.48912811, 52.48912811, + 52.69416809]], dtype=np.float32) +) GAIN = 0.20503567620766011 OFFSET = -10.456819486590666 @@ -49,152 +59,303 @@ CHANNEL_NAME = 'IR_108' PLATFORM_ID = 323 # Met-10 -TBS_OUTPUT1 = np.array([[269.29684448, 269.29684448, 269.29684448, 269.13296509, - 268.96871948], - [269.13296509, 268.96871948, 269.13296509, 268.80422974, - 268.80422974], - [268.80422974, 268.63937378, 268.63937378, 268.80422974, - 268.80422974], - [264.23751831, 263.88912964, 263.88912964, 264.41116333, - 264.23751831], - [256.77682495, 256.77682495, 256.96743774, 256.96743774, - 257.15756226]], dtype=np.float32) - - -TBS_OUTPUT2 = np.array([[268.94519043, 268.94519043, 268.94519043, 268.77984619, - 268.61422729], - [268.77984619, 268.61422729, 268.77984619, 268.44830322, - 268.44830322], - [268.44830322, 268.28204346, 268.28204346, 268.44830322, - 268.44830322], - [263.84396362, 263.49285889, 263.49285889, 264.01898193, - 263.84396362], - [256.32858276, 256.32858276, 256.52044678, 256.52044678, - 256.71188354]], dtype=np.float32) - +TBS_OUTPUT1 = xr.DataArray( + np.array([[269.29684448, 269.29684448, 269.29684448, 269.13296509, + 268.96871948], + [269.13296509, 268.96871948, 269.13296509, 268.80422974, + 268.80422974], + [268.80422974, 268.63937378, 268.63937378, 268.80422974, + 268.80422974], + [264.23751831, 263.88912964, 263.88912964, 264.41116333, + 264.23751831], + [256.77682495, 256.77682495, 256.96743774, 256.96743774, + 257.15756226]], dtype=np.float32) +) + +TBS_OUTPUT2 = xr.DataArray( + np.array([[268.94519043, 268.94519043, 268.94519043, 268.77984619, + 268.61422729], + [268.77984619, 268.61422729, 268.77984619, 268.44830322, + 268.44830322], + [268.44830322, 268.28204346, 268.28204346, 268.44830322, + 268.44830322], + [263.84396362, 263.49285889, 263.49285889, 264.01898193, + 263.84396362], + [256.32858276, 256.32858276, 256.52044678, 256.52044678, + 256.71188354]], dtype=np.float32) +) VIS008_SOLAR_IRRADIANCE = 73.1807 -VIS008_RADIANCE = np.array([[0.62234485, 0.59405649, 0.59405649, 0.59405649, 0.59405649], - [0.59405649, 0.62234485, 0.62234485, 0.59405649, 0.62234485], - [0.76378691, 0.79207528, 0.79207528, 0.76378691, 0.79207528], - [3.30974245, 3.33803129, 3.33803129, 3.25316572, 3.47947311], - [7.52471399, 7.83588648, 8.2602129, 8.57138538, 8.99571133]], dtype=np.float32) -VIS008_RADIANCE = xr.DataArray(VIS008_RADIANCE) - -VIS008_REFLECTANCE = np.array([[2.8066392, 2.6790648, 2.6790648, 2.6790648, - 2.6790648], - [2.6790648, 2.8066392, 2.8066392, 2.6790648, - 2.8066392], - [3.444512, 3.572086, 3.572086, 3.444512, - 3.572086], - [14.926213, 15.053792, 15.053792, 14.671064, - 15.691662], - [33.934814, 35.33813, 37.251755, 38.655075, - 40.56869]], dtype=np.float32) - -# -- - -CAL_DTYPE = np.array([[(0.0208876, -1.06526761), (0.0278805, -1.42190546), - (0.0235881, -1.20299312), (0.00365867, -0.18659201), - (0.00831811, -0.42422367), (0.03862197, -1.96972038), - (0.12674432, -6.46396025), (0.10396091, -5.30200645), - (0.20503568, -10.45681949), (0.22231115, -11.33786848), - (0.1576069, -8.03795174), (0.0373969, -1.90724192)]], - dtype=[('CalSlope', '>f8'), ('CalOffset', '>f8')]) - -IR_108_RADIANCES = np.ma.array([[133.06815651, 133.68326355, 134.29837059, 134.91347763, - 135.52858467], - [136.14369171, 136.75879875, 137.37390579, 137.98901283, - 138.60411987], - [139.21922691, 139.83433395, 140.44944099, 141.06454803, - 141.67965507]], - mask=False, dtype=np.float64) - -VIS006_RADIANCES = np.ma.array([[13.55605239, 13.61871519, 13.68137799, 13.74404079, - 13.80670359], - [13.86936639, 13.93202919, 13.99469199, 14.05735479, - 14.12001759], - [14.18268039, 14.24534319, 14.30800599, 14.37066879, - 14.43333159]], mask=False, dtype=np.float64) - -VIS006_REFLECTANCES = np.array([[65.00454035, 65.30502359, 65.60550682, 65.90599006, - 66.2064733], - [66.50695654, 66.80743977, 67.10792301, 67.40840625, - 67.70888949], - [68.00937272, 68.30985596, 68.6103392, 68.91082244, - 69.21130567]], dtype=np.float64) - -IR_108_TBS = np.array([[311.77913132, 312.11070275, 312.44143083, 312.77132215, - 313.10038322], - [313.42862046, 313.75604023, 314.0826488, 314.40845236, - 314.73345704], - [315.05766888, 315.38109386, 315.70373788, 316.02560677, - 316.34670629]], dtype=np.float64) - - -# Calibration type = Effective radiances -CALIBRATION_TYPE = np.array( - [[2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]], dtype=np.uint8) - - -# This should preferably be put in a helper-module -# Fixme! -def assertNumpyArraysEqual(self, other): - if self.shape != other.shape: - raise AssertionError("Shapes don't match") - if not np.allclose(self, other): - raise AssertionError("Elements don't match!") - - -class TestSEVIRICalibrationHandler(unittest.TestCase): - - """Test the SEVIRICalibrationHandler class in the msg_base module""" +VIS008_RADIANCE = xr.DataArray( + np.array([[0.62234485, 0.59405649, 0.59405649, 0.59405649, 0.59405649], + [0.59405649, 0.62234485, 0.62234485, 0.59405649, 0.62234485], + [0.76378691, 0.79207528, 0.79207528, 0.76378691, 0.79207528], + [3.30974245, 3.33803129, 3.33803129, 3.25316572, 3.47947311], + [7.52471399, 7.83588648, 8.2602129, 8.57138538, 8.99571133]], + dtype=np.float32) +) - def setUp(self): - """Setup the SEVIRI Calibration handler for testing.""" +VIS008_REFLECTANCE = xr.DataArray( + np.array([[2.8066392, 2.6790648, 2.6790648, 2.6790648, 2.6790648], + [2.6790648, 2.8066392, 2.8066392, 2.6790648, 2.8066392], + [3.444512, 3.572086, 3.572086, 3.444512, 3.572086], + [14.926213, 15.053792, 15.053792, 14.671064, 15.691662], + [33.934814, 35.33813, 37.251755, 38.655075, 40.56869]], + dtype=np.float32) +) - hdr = {} - hdr['15_DATA_HEADER'] = {} - hdr['15_DATA_HEADER']['RadiometricProcessing'] = { - 'Level15ImageCalibration': CAL_DTYPE} - hdr['15_DATA_HEADER']['ImageDescription'] = {} - hdr['15_DATA_HEADER']['ImageDescription']['Level15ImageProduction'] = { - 'PlannedChanProcessing': CALIBRATION_TYPE} +class TestSEVIRICalibrationAlgorithm(unittest.TestCase): + """Unit Tests for SEVIRI calibration algorithm.""" - self.handler = SEVIRICalibrationHandler() - self.handler.platform_id = PLATFORM_ID - self.handler.start_time = datetime(2020, 8, 15, 13, 0, 40) + def setUp(self): + """Set up the SEVIRI Calibration algorithm for testing.""" + self.algo = SEVIRICalibrationAlgorithm( + platform_id=PLATFORM_ID, + scan_time=datetime(2020, 8, 15, 13, 0, 40) + ) def test_convert_to_radiance(self): - """Test the conversion from counts to radiance method""" - - data = COUNTS_INPUT - gain = GAIN - offset = OFFSET - result = self.handler._convert_to_radiance(data, gain, offset) - assertNumpyArraysEqual(result, RADIANCES_OUTPUT) + """Test the conversion from counts to radiances.""" + result = self.algo.convert_to_radiance(COUNTS_INPUT, GAIN, OFFSET) + xr.testing.assert_allclose(result, RADIANCES_OUTPUT) + self.assertEqual(result.dtype, np.float32) def test_ir_calibrate(self): + """Test conversion from radiance to brightness temperature.""" + result = self.algo.ir_calibrate(RADIANCES_OUTPUT, + CHANNEL_NAME, CAL_TYPE1) + xr.testing.assert_allclose(result, TBS_OUTPUT1, rtol=1E-5) + self.assertEqual(result.dtype, np.float32) - result = self.handler._ir_calibrate(RADIANCES_OUTPUT, - CHANNEL_NAME, CAL_TYPE1) - assertNumpyArraysEqual(result, TBS_OUTPUT1) - - result = self.handler._ir_calibrate(RADIANCES_OUTPUT, - CHANNEL_NAME, CAL_TYPE2) - assertNumpyArraysEqual(result, TBS_OUTPUT2) + result = self.algo.ir_calibrate(RADIANCES_OUTPUT, + CHANNEL_NAME, CAL_TYPE2) + xr.testing.assert_allclose(result, TBS_OUTPUT2, rtol=1E-5) with self.assertRaises(NotImplementedError): - self.handler._ir_calibrate(RADIANCES_OUTPUT, - CHANNEL_NAME, CAL_TYPEBAD) + self.algo.ir_calibrate(RADIANCES_OUTPUT, CHANNEL_NAME, CAL_TYPEBAD) def test_vis_calibrate(self): - result = self.handler._vis_calibrate(VIS008_RADIANCE, - VIS008_SOLAR_IRRADIANCE) - assertNumpyArraysEqual(result, VIS008_REFLECTANCE) + """Test conversion from radiance to reflectance.""" + result = self.algo.vis_calibrate(VIS008_RADIANCE, + VIS008_SOLAR_IRRADIANCE) + xr.testing.assert_allclose(result, VIS008_REFLECTANCE) self.assertTrue(result.sun_earth_distance_correction_applied) - - def tearDown(self): - pass + self.assertEqual(result.dtype, np.float32) + + +class TestSeviriCalibrationHandler: + """Unit tests for SEVIRI calibration handler.""" + + def test_init(self): + """Test initialization of the calibration handler.""" + with pytest.raises(ValueError): + SEVIRICalibrationHandler( + platform_id=None, + channel_name=None, + coefs=None, + calib_mode='invalid', + scan_time=None + ) + + def _get_calibration_handler(self, calib_mode='NOMINAL', ext_coefs=None): + """Provide a calibration handler.""" + return SEVIRICalibrationHandler( + platform_id=324, + channel_name='IR_108', + coefs={ + 'coefs': { + 'NOMINAL': { + 'gain': 10, + 'offset': -1 + }, + 'GSICS': { + 'gain': 20, + 'offset': -2 + }, + 'EXTERNAL': ext_coefs or {} + }, + 'radiance_type': 1 + }, + calib_mode=calib_mode, + scan_time=None + ) + + def test_calibrate_exceptions(self): + """Test exceptions raised by the calibration handler.""" + calib = self._get_calibration_handler() + with pytest.raises(ValueError): + calib.calibrate(None, 'invalid') + + @pytest.mark.parametrize( + ('calib_mode', 'ext_coefs', 'expected'), + [ + ('NOMINAL', {}, (10, -1)), + ('GSICS', {}, (20, -40)), + ('GSICS', {'gain': 30, 'offset': -3}, (30, -3)), + ('NOMINAL', {'gain': 30, 'offset': -3}, (30, -3)) + ] + ) + def test_get_gain_offset(self, calib_mode, ext_coefs, expected): + """Test selection of gain and offset.""" + calib = self._get_calibration_handler(calib_mode=calib_mode, + ext_coefs=ext_coefs) + coefs = calib.get_gain_offset() + assert coefs == expected + + +class TestFileHandlerCalibrationBase: + """Base class for file handler calibration tests.""" + + platform_id = 324 + gains_nominal = np.arange(1, 13) + offsets_nominal = np.arange(-1, -13, -1) + # No GSICS coefficients for VIS channels -> set to zero + gains_gsics = [0, 0, 0, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 0] + offsets_gsics = [0, 0, 0, -0.4, -0.5, -0.6, -0.7, -0.8, -0.9, -1.0, -1.1, 0] + radiance_types = 2 * np.ones(12) + scan_time = datetime(2020, 1, 1) + external_coefs = { + 'VIS006': {'gain': 10, 'offset': -10}, + 'IR_108': {'gain': 20, 'offset': -20}, + 'HRV': {'gain': 5, 'offset': -5} + } + spectral_channel_ids = {'VIS006': 1, 'IR_108': 9, 'HRV': 12} + expected = { + 'VIS006': { + 'counts': { + 'NOMINAL': xr.DataArray( + [[0, 10], + [100, 255]], + dims=('y', 'x') + ) + }, + 'radiance': { + 'NOMINAL': xr.DataArray( + [[np.nan, 9], + [99, 254]], + dims=('y', 'x') + ), + 'GSICS': xr.DataArray( + [[np.nan, 9], + [99, 254]], + dims=('y', 'x') + ), + 'EXTERNAL': xr.DataArray( + [[np.nan, 90], + [990, 2540]], + dims=('y', 'x') + ) + }, + 'reflectance': { + 'NOMINAL': xr.DataArray( + [[np.nan, 40.47923], + [445.27155, 1142.414]], + dims=('y', 'x') + ), + 'EXTERNAL': xr.DataArray( + [[np.nan, 404.7923], + [4452.7153, 11424.14]], + dims=('y', 'x') + ) + } + }, + 'IR_108': { + 'counts': { + 'NOMINAL': xr.DataArray( + [[0, 10], + [100, 255]], + dims=('y', 'x') + ) + }, + 'radiance': { + 'NOMINAL': xr.DataArray( + [[np.nan, 81], + [891, 2286]], + dims=('y', 'x') + ), + 'GSICS': xr.DataArray( + [[np.nan, 8.19], + [89.19, 228.69]], + dims=('y', 'x') + ), + 'EXTERNAL': xr.DataArray( + [[np.nan, 180], + [1980, 5080]], + dims=('y', 'x') + ) + }, + 'brightness_temperature': { + 'NOMINAL': xr.DataArray( + [[np.nan, 279.82318], + [543.2585, 812.77167]], + dims=('y', 'x') + ), + 'GSICS': xr.DataArray( + [[np.nan, 189.20985], + [285.53293, 356.06668]], + dims=('y', 'x') + ), + 'EXTERNAL': xr.DataArray( + [[np.nan, 335.14236], + [758.6249, 1262.7567]], + dims=('y', 'x') + ), + } + }, + 'HRV': { + 'counts': { + 'NOMINAL': xr.DataArray( + [[0, 10], + [100, 255]], + dims=('y', 'x') + ) + }, + 'radiance': { + 'NOMINAL': xr.DataArray( + [[np.nan, 108], + [1188, 3048]], + dims=('y', 'x') + ), + 'GSICS': xr.DataArray( + [[np.nan, 108], + [1188, 3048]], + dims=('y', 'x') + ), + 'EXTERNAL': xr.DataArray( + [[np.nan, 45], + [495, 1270]], + dims=('y', 'x') + ) + }, + 'reflectance': { + 'NOMINAL': xr.DataArray( + [[np.nan, 401.28372], + [4414.121, 11325.118]], + dims=('y', 'x') + ), + 'EXTERNAL': xr.DataArray( + [[np.nan, 167.20154], + [1839.217, 4718.799]], + dims=('y', 'x') + ) + } + } + } + + @pytest.fixture(name='counts') + def counts(self): + """Provide fake image counts.""" + return xr.DataArray( + [[0, 10], + [100, 255]], + dims=('y', 'x') + ) + + def _get_expected( + self, channel, calibration, calib_mode, use_ext_coefs + ): + if use_ext_coefs: + return self.expected[channel][calibration]['EXTERNAL'] + return self.expected[channel][calibration][calib_mode] diff --git a/satpy/tests/reader_tests/test_seviri_l1b_hrit.py b/satpy/tests/reader_tests/test_seviri_l1b_hrit.py index b469ad612a..529dc8a1ac 100644 --- a/satpy/tests/reader_tests/test_seviri_l1b_hrit.py +++ b/satpy/tests/reader_tests/test_seviri_l1b_hrit.py @@ -22,12 +22,17 @@ from datetime import datetime import numpy as np +import pytest import xarray as xr -from satpy.readers.seviri_l1b_hrit import (HRITMSGFileHandler, HRITMSGPrologueFileHandler, HRITMSGEpilogueFileHandler, - NoValidOrbitParams) -from satpy.readers.seviri_base import CHANNEL_NAMES, VIS_CHANNELS +from satpy.readers.seviri_l1b_hrit import ( + HRITMSGFileHandler, HRITMSGPrologueFileHandler, HRITMSGEpilogueFileHandler, + NoValidOrbitParams +) from satpy.tests.utils import make_dataid +from satpy.tests.reader_tests.test_seviri_l1b_calibration import ( + TestFileHandlerCalibrationBase +) from numpy import testing as npt @@ -343,89 +348,6 @@ def test_read_band(self, memmap): res = self.reader.read_band('VIS006', None) self.assertEqual(res.shape, (464, 3712)) - @mock.patch('satpy.readers.hrit_base.HRITFileHandler.__init__', return_value=None) - @mock.patch('satpy.readers.seviri_l1b_hrit.HRITMSGFileHandler._get_header', autospec=True) - @mock.patch('satpy.readers.seviri_base.SEVIRICalibrationHandler._convert_to_radiance') - def test_calibrate(self, _convert_to_radiance, get_header, *mocks): - """Test selection of calibration coefficients.""" - shp = (10, 10) - counts = xr.DataArray(np.zeros(shp)) - nominal_gain = np.arange(1, 13) - nominal_offset = np.arange(-1, -13, -1) - gsics_gain = np.arange(0.1, 1.3, 0.1) - gsics_offset = np.arange(-0.1, -1.3, -0.1) - - # Mock prologue & epilogue - pro = mock.MagicMock(prologue={'RadiometricProcessing': { - 'Level15ImageCalibration': {'CalSlope': nominal_gain, - 'CalOffset': nominal_offset}, - 'MPEFCalFeedback': {'GSICSCalCoeff': gsics_gain, - 'GSICSOffsetCount': gsics_offset} - }}) - epi = mock.MagicMock(epilogue=None) - - # Mock header readout - mda = {'image_segment_line_quality': {'line_validity': np.zeros(shp[0]), - 'line_radiometric_quality': np.zeros(shp[0]), - 'line_geometric_quality': np.zeros(shp[0])}} - - def get_header_patched(self): - self.mda = mda - - get_header.side_effect = get_header_patched - - # Test selection of calibration coefficients - # - # a) Default: Nominal calibration - reader = HRITMSGFileHandler(filename=None, filename_info=None, filetype_info=None, - prologue=pro, epilogue=epi) - for ch_id, ch_name in CHANNEL_NAMES.items(): - reader.channel_name = ch_name - reader.mda['spectral_channel_id'] = ch_id - reader.calibrate(data=counts, calibration='radiance') - _convert_to_radiance.assert_called_with(mock.ANY, nominal_gain[ch_id - 1], - nominal_offset[ch_id - 1]) - - # b) GSICS calibration for IR channels, nominal calibration for VIS channels - reader = HRITMSGFileHandler(filename=None, filename_info=None, filetype_info=None, - prologue=pro, epilogue=epi, calib_mode='GSICS') - for ch_id, ch_name in CHANNEL_NAMES.items(): - if ch_name in VIS_CHANNELS: - gain, offset = nominal_gain[ch_id - 1], nominal_offset[ch_id - 1] - else: - gain, offset = gsics_gain[ch_id - 1], gsics_offset[ch_id - 1] - offset = offset * gain - - reader.channel_name = ch_name - reader.mda['spectral_channel_id'] = ch_id - reader.calibrate(data=counts, calibration='radiance') - _convert_to_radiance.assert_called_with(mock.ANY, gain, offset) - - # c) External calibration coefficients for selected channels, GSICS coefs for remaining - # IR channels, nominal coefs for remaining VIS channels - coefs = {'VIS006': {'gain': 1.234, 'offset': -0.1}, - 'IR_108': {'gain': 2.345, 'offset': -0.2}} - reader = HRITMSGFileHandler(filename=None, filename_info=None, filetype_info=None, - prologue=pro, epilogue=epi, ext_calib_coefs=coefs, - calib_mode='GSICS') - for ch_id, ch_name in CHANNEL_NAMES.items(): - if ch_name in coefs: - gain, offset = coefs[ch_name]['gain'], coefs[ch_name]['offset'] - elif ch_name not in VIS_CHANNELS: - gain, offset = gsics_gain[ch_id - 1], gsics_offset[ch_id - 1] - offset = offset * gain - else: - gain, offset = nominal_gain[ch_id - 1], nominal_offset[ch_id - 1] - - reader.channel_name = ch_name - reader.mda['spectral_channel_id'] = ch_id - reader.calibrate(data=counts, calibration='radiance') - _convert_to_radiance.assert_called_with(mock.ANY, gain, offset) - - # d) Invalid mode - self.assertRaises(ValueError, HRITMSGFileHandler, filename=None, filename_info=None, - filetype_info=None, prologue=pro, epilogue=epi, calib_mode='invalid') - @mock.patch('satpy.readers.seviri_l1b_hrit.HRITMSGFileHandler._get_timestamps') @mock.patch('satpy.readers.seviri_l1b_hrit.HRITFileHandler.get_dataset') @mock.patch('satpy.readers.seviri_l1b_hrit.HRITMSGFileHandler.calibrate') @@ -738,3 +660,105 @@ def test_reduce(self, reduce_mda): self.reader._reduced = 'red' self.assertEqual(self.reader.reduce(123), 'red') reduce_mda.assert_not_called() + + +class TestHRITMSGCalibration(TestFileHandlerCalibrationBase): + """Unit tests for calibration.""" + + @pytest.fixture(name='file_handler') + def file_handler(self): + """Create a mocked file handler.""" + prolog = { + 'RadiometricProcessing': { + 'Level15ImageCalibration': { + 'CalSlope': self.gains_nominal, + 'CalOffset': self.offsets_nominal, + }, + 'MPEFCalFeedback': { + 'GSICSCalCoeff': self.gains_gsics, + 'GSICSOffsetCount': self.offsets_gsics, + } + }, + 'ImageDescription': { + 'Level15ImageProduction': { + 'PlannedChanProcessing': self.radiance_types + } + } + } + epilog = { + 'ImageProductionStats': { + 'ActualScanningSummary': { + 'ForwardScanStart': self.scan_time + } + } + } + mda = { + 'image_segment_line_quality': { + 'line_validity': np.zeros(2), + 'line_radiometric_quality': np.zeros(2), + 'line_geometric_quality': np.zeros(2) + }, + } + + with mock.patch( + 'satpy.readers.seviri_l1b_hrit.HRITMSGFileHandler.__init__', + return_value=None + ): + fh = HRITMSGFileHandler() + fh.platform_id = self.platform_id + fh.mda = mda + fh.prologue = prolog + fh.epilogue = epilog + return fh + + @pytest.mark.parametrize( + ('channel', 'calibration', 'calib_mode', 'use_ext_coefs'), + [ + # VIS channel, internal coefficients + ('VIS006', 'counts', 'NOMINAL', False), + ('VIS006', 'radiance', 'NOMINAL', False), + ('VIS006', 'radiance', 'GSICS', False), + ('VIS006', 'reflectance', 'NOMINAL', False), + # VIS channel, external coefficients (mode should have no effect) + ('VIS006', 'radiance', 'GSICS', True), + ('VIS006', 'reflectance', 'NOMINAL', True), + # IR channel, internal coefficients + ('IR_108', 'counts', 'NOMINAL', False), + ('IR_108', 'radiance', 'NOMINAL', False), + ('IR_108', 'radiance', 'GSICS', False), + ('IR_108', 'brightness_temperature', 'NOMINAL', False), + ('IR_108', 'brightness_temperature', 'GSICS', False), + # IR channel, external coefficients (mode should have no effect) + ('IR_108', 'radiance', 'NOMINAL', True), + ('IR_108', 'brightness_temperature', 'GSICS', True), + # HRV channel, internal coefficiens + ('HRV', 'counts', 'NOMINAL', False), + ('HRV', 'radiance', 'NOMINAL', False), + ('HRV', 'radiance', 'GSICS', False), + ('HRV', 'reflectance', 'NOMINAL', False), + # HRV channel, external coefficients (mode should have no effect) + ('HRV', 'radiance', 'GSICS', True), + ('HRV', 'reflectance', 'NOMINAL', True), + ] + ) + def test_calibrate( + self, file_handler, counts, channel, calibration, calib_mode, + use_ext_coefs + ): + """Test the calibration.""" + external_coefs = self.external_coefs if use_ext_coefs else {} + expected = self._get_expected( + channel=channel, + calibration=calibration, + calib_mode=calib_mode, + use_ext_coefs=use_ext_coefs + ) + + fh = file_handler + fh.mda['spectral_channel_id'] = self.spectral_channel_ids[channel] + fh.channel_name = channel + fh.calib_mode = calib_mode + fh.ext_calib_coefs = external_coefs + + res = fh.calibrate(counts, calibration) + xr.testing.assert_allclose(res, expected) diff --git a/satpy/tests/reader_tests/test_seviri_l1b_native.py b/satpy/tests/reader_tests/test_seviri_l1b_native.py index 8d06480a6c..1e0ee43f44 100644 --- a/satpy/tests/reader_tests/test_seviri_l1b_native.py +++ b/satpy/tests/reader_tests/test_seviri_l1b_native.py @@ -19,14 +19,18 @@ import unittest from unittest import mock + import numpy as np +import pytest import xarray as xr from satpy.readers.seviri_l1b_native import ( NativeMSGFileHandler, ImageBoundaries, Padder, get_available_channels, ) - +from satpy.tests.reader_tests.test_seviri_l1b_calibration import ( + TestFileHandlerCalibrationBase +) from satpy.tests.utils import make_dataid CHANNEL_INDEX_LIST = ['VIS006', 'VIS008', 'IR_016', 'IR_039', @@ -1020,180 +1024,94 @@ def test_is_roi_roi(self): self.assertEqual(calculated, expected) -class TestNativeMSGCalibrationMode(unittest.TestCase): - """Test NativeMSGFileHandler.get_area_extent. - - The expected results have been verified by manually - inspecting the output of geoferenced imagery. - """ - - @staticmethod - def create_test_header(earth_model, dataset_id, is_full_disk, is_rapid_scan): - """Create Test Header. - - Mocked NativeMSGFileHandler with sufficient attributes for - NativeMSGFileHandler._convert_to_radiance and NativeMSGFileHandler.calibrate to be able to execute. - """ - if dataset_id['name'] == 'HRV': - # reference_grid = 'ReferenceGridHRV' - column_dir_grid_step = 1.0001343488693237 - line_dir_grid_step = 1.0001343488693237 - else: - # reference_grid = 'ReferenceGridVIS_IR' - column_dir_grid_step = 3.0004031658172607 - line_dir_grid_step = 3.0004031658172607 - - if is_full_disk: - north = 3712 - east = 1 - west = 3712 - south = 1 - n_visir_cols = 3712 - n_visir_lines = 3712 - n_hrv_lines = 11136 - elif is_rapid_scan: - north = 3712 - east = 1 - west = 3712 - south = 2321 - n_visir_cols = 3712 - n_visir_lines = 1392 - n_hrv_lines = 4176 - else: - north = 3574 - east = 78 - west = 2591 - south = 1746 - n_visir_cols = 2516 - n_visir_lines = north - south + 1 - n_hrv_lines = 11136 +class TestNativeMSGCalibration(TestFileHandlerCalibrationBase): + """Unit tests for calibration.""" + @pytest.fixture(name='file_handler') + def file_handler(self): + """Create a mocked file handler.""" header = { '15_DATA_HEADER': { - 'ImageDescription': { - 'reference_grid': { - 'ColumnDirGridStep': column_dir_grid_step, - 'LineDirGridStep': line_dir_grid_step, - 'GridOrigin': 2, # south-east corner - }, - 'ProjectionDescription': { - 'LongitudeOfSSP': 0.0 - } - }, - 'GeometricProcessing': { - 'EarthModel': { - 'TypeOfEarthModel': earth_model, - 'EquatorialRadius': 6378169.0, - 'NorthPolarRadius': 6356583.800000001, - 'SouthPolarRadius': 6356583.800000001, - } - }, 'RadiometricProcessing': { 'Level15ImageCalibration': { - 'CalSlope': [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.96, 0.97], - 'CalOffset': [-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0], + 'CalSlope': self.gains_nominal, + 'CalOffset': self.offsets_nominal, }, 'MPEFCalFeedback': { - 'GSICSCalCoeff': [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, - 0.7, 0.8, 0.9, 0.95, 0.96, 0.97], - 'GSICSOffsetCount': [-51.0, -51.0, -51.0, -51.0, -51.0, -51.0, - -51.0, -51.0, -51.0, -51.0, -51.0, -51.0] - }, + 'GSICSCalCoeff': self.gains_gsics, + 'GSICSOffsetCount': self.offsets_gsics + } }, - 'SatelliteStatus': { - 'SatelliteDefinition': { - 'SatelliteId': 324 + 'ImageDescription': { + 'Level15ImageProduction': { + 'PlannedChanProcessing': self.radiance_types + } + }, + 'ImageAcquisition': { + 'PlannedAcquisitionTime': { + 'TrueRepeatCycleStart': self.scan_time } } - }, - '15_SECONDARY_PRODUCT_HEADER': { - 'NorthLineSelectedRectangle': {'Value': north}, - 'EastColumnSelectedRectangle': {'Value': east}, - 'WestColumnSelectedRectangle': {'Value': west}, - 'SouthLineSelectedRectangle': {'Value': south}, - 'SelectedBandIDs': {'Value': 'xxxxxxxxxxxx'}, - 'NumberColumnsVISIR': {'Value': n_visir_cols}, - 'NumberLinesVISIR': {'Value': n_visir_lines}, - 'NumberColumnsHRV': {'Value': 11136}, - 'NumberLinesHRV': {'Value': n_hrv_lines}, } - } + with mock.patch('satpy.readers.seviri_l1b_native.NativeMSGFileHandler.__init__', + return_value=None): + fh = NativeMSGFileHandler() + fh.header = header + fh.platform_id = self.platform_id + return fh + + @pytest.mark.parametrize( + ('channel', 'calibration', 'calib_mode', 'use_ext_coefs'), + [ + # VIS channel, internal coefficients + ('VIS006', 'counts', 'NOMINAL', False), + ('VIS006', 'radiance', 'NOMINAL', False), + ('VIS006', 'radiance', 'GSICS', False), + ('VIS006', 'reflectance', 'NOMINAL', False), + # VIS channel, external coefficients (mode should have no effect) + ('VIS006', 'radiance', 'GSICS', True), + ('VIS006', 'reflectance', 'NOMINAL', True), + # IR channel, internal coefficients + ('IR_108', 'counts', 'NOMINAL', False), + ('IR_108', 'radiance', 'NOMINAL', False), + ('IR_108', 'radiance', 'GSICS', False), + ('IR_108', 'brightness_temperature', 'NOMINAL', False), + ('IR_108', 'brightness_temperature', 'GSICS', False), + # IR channel, external coefficients (mode should have no effect) + ('IR_108', 'radiance', 'NOMINAL', True), + ('IR_108', 'brightness_temperature', 'GSICS', True), + # HRV channel, internal coefficiens + ('HRV', 'counts', 'NOMINAL', False), + ('HRV', 'radiance', 'NOMINAL', False), + ('HRV', 'radiance', 'GSICS', False), + ('HRV', 'reflectance', 'NOMINAL', False), + # HRV channel, external coefficients (mode should have no effect) + ('HRV', 'radiance', 'GSICS', True), + ('HRV', 'reflectance', 'NOMINAL', True), + ] + ) + def test_calibrate( + self, file_handler, counts, channel, calibration, calib_mode, + use_ext_coefs + ): + """Test the calibration.""" + external_coefs = self.external_coefs if use_ext_coefs else {} + expected = self._get_expected( + channel=channel, + calibration=calibration, + calib_mode=calib_mode, + use_ext_coefs=use_ext_coefs + ) - return header - - def calibration_mode_test(self, test_dict, cal_mode): - """Test the Calibration Mode.""" - # dummy data array - data = xr.DataArray([255., 200., 300.]) - - earth_model = test_dict['earth_model'] - dataset_id = test_dict['dataset_id'] - index = CHANNEL_INDEX_LIST.index(dataset_id['name']) - - # determine the cal coeffs needed for the expected data calculation - if cal_mode == 'nominal': - cal_slope = test_dict['CalSlope'][index] - cal_offset = test_dict['CalOffset'][index] - else: - cal_slope_arr = test_dict['GSICSCalCoeff'] - cal_offset_arr = test_dict['GSICSOffsetCount'] - cal_offset = cal_offset_arr[index] * cal_slope_arr[index] - cal_slope = cal_slope_arr[index] - - is_full_disk = test_dict['is_full_disk'] - is_rapid_scan = test_dict['is_rapid_scan'] - header = self.create_test_header(earth_model, dataset_id, is_full_disk, is_rapid_scan) - - with mock.patch('satpy.readers.seviri_l1b_native.np.fromfile') as fromfile: - fromfile.return_value = header - with mock.patch('satpy.readers.seviri_l1b_native.recarray2dict') as recarray2dict: - recarray2dict.side_effect = (lambda x: x) - with mock.patch('satpy.readers.seviri_l1b_native.NativeMSGFileHandler._get_memmap') as _get_memmap: - _get_memmap.return_value = np.arange(3) - with mock.patch('satpy.readers.seviri_l1b_native.NativeMSGFileHandler._read_trailer'): - # Create an instance of the native msg reader - # with the calibration mode to test - fh = NativeMSGFileHandler(None, {}, None, calib_mode=cal_mode) - - # Calculate the expected calibration values using the coeffs - # from the test data set - expected = fh._convert_to_radiance(data, cal_slope, cal_offset) - - # Calculate the calibrated values using the cal coeffs from the - # test header and using the correct calibration mode values - fh.header = header - calculated = fh.calibrate(data, dataset_id) - - return (expected.data, calculated.data) + fh = file_handler + fh.calib_mode = calib_mode + fh.ext_calib_coefs = external_coefs - def test_calibration_mode_nominal(self): - """Test the nominal calibration mode.""" - # Test using the Nominal calibration mode - expected, calculated = self.calibration_mode_test( - TEST_CALIBRATION_MODE, - 'nominal', - ) - assertNumpyArraysEqual(calculated, expected) - - def test_calibration_mode_gsics(self): - """Test the GSICS calibration mode.""" - # Test using the GSICS calibration mode - expected, calculated = self.calibration_mode_test( - TEST_CALIBRATION_MODE, - 'gsics', - ) - assertNumpyArraysEqual(calculated, expected) - - def test_calibration_mode_dummy(self): - """Test a dummy calibration mode.""" - # pass in a calibration mode that is not recognised by the reader - # and an exception will be raised - self.assertRaises(NotImplementedError, self.calibration_mode_test, - TEST_CALIBRATION_MODE, - 'dummy', - ) + dataset_id = make_dataid(name=channel, calibration=calibration) + res = fh.calibrate(counts, dataset_id) + xr.testing.assert_allclose(res, expected) class TestNativeMSGPadder(unittest.TestCase): diff --git a/satpy/tests/reader_tests/test_seviri_l1b_nc.py b/satpy/tests/reader_tests/test_seviri_l1b_nc.py index 93c3741983..58217ee513 100644 --- a/satpy/tests/reader_tests/test_seviri_l1b_nc.py +++ b/satpy/tests/reader_tests/test_seviri_l1b_nc.py @@ -22,9 +22,13 @@ from datetime import datetime import numpy as np +import pytest import xarray as xr from satpy.readers.seviri_l1b_nc import NCSEVIRIFileHandler +from satpy.tests.reader_tests.test_seviri_l1b_calibration import ( + TestFileHandlerCalibrationBase +) from satpy.tests.utils import make_dataid @@ -54,13 +58,15 @@ class TestNCSEVIRIFileHandler(unittest.TestCase): def setUp(self): """Set up the test case.""" + start_time = datetime(2016, 3, 3, 0, 0) with mock.patch.object(NCSEVIRIFileHandler, '_read_file', new=new_read_file): self.reader = NCSEVIRIFileHandler( 'filename', {'platform_shortname': 'MSG3', - 'start_time': datetime(2016, 3, 3, 0, 0), + 'start_time': start_time, 'service': 'MSG'}, {'filetype': 'info'}) + self.reader.deltaSt = start_time def test_get_dataset_remove_attrs(self): """Test getting the hrv dataset.""" @@ -71,3 +77,110 @@ def test_get_dataset_remove_attrs(self): strip_attrs = ["comment", "long_name", "nc_key", "scale_factor", "add_offset", "valid_min", "valid_max"] self.assertFalse(any([k in res.attrs.keys() for k in strip_attrs])) + + +class TestCalibration(TestFileHandlerCalibrationBase): + """Unit tests for calibration.""" + + @pytest.fixture(name='file_handler') + def file_handler(self, counts): + """Create a mocked file handler.""" + with mock.patch( + 'satpy.readers.seviri_l1b_nc.NCSEVIRIFileHandler.__init__', + return_value=None + ): + # Create dataset and set calibration coefficients + ds = xr.Dataset( + { + 'VIS006': counts.copy(), + 'IR_108': counts.copy(), + 'planned_chan_processing': self.radiance_types, + }, + attrs={ + 'satellite_id': self.platform_id + } + ) + ds['VIS006'].attrs.update({ + 'scale_factor': self.gains_nominal[0], + 'add_offset': self.offsets_nominal[0] + }) + ds['IR_108'].attrs.update({ + 'scale_factor': self.gains_nominal[8], + 'add_offset': self.offsets_nominal[8], + }) + + # Add some attributes so that the reader can strip them + strip_attrs = { + 'comment': None, + 'long_name': None, + 'valid_min': None, + 'valid_max': None + } + for name in ['VIS006', 'IR_108']: + ds[name].attrs.update(strip_attrs) + + # Create file handler + fh = NCSEVIRIFileHandler() + fh.nc = ds + fh.deltaSt = self.scan_time + fh.mda = { + 'projection_parameters': { + 'ssp_longitude': 0, + 'h': 12345 + } + } + return fh + + @pytest.mark.parametrize( + ('channel', 'calibration', 'use_ext_coefs'), + [ + # VIS channel, internal coefficients + ('VIS006', 'counts', False), + ('VIS006', 'radiance', False), + ('VIS006', 'reflectance', False), + # VIS channel, external coefficients + ('VIS006', 'radiance', True), + ('VIS006', 'reflectance', True), + # IR channel, internal coefficients + ('IR_108', 'counts', False), + ('IR_108', 'radiance', False), + ('IR_108', 'brightness_temperature', False), + # IR channel, external coefficients + ('IR_108', 'radiance', True), + ('IR_108', 'brightness_temperature', True), + # FUTURE: Enable once HRV reading has been fixed. + # # HRV channel, internal coefficiens + # ('HRV', 'counts', False), + # ('HRV', 'radiance', False), + # ('HRV', 'reflectance', False), + # # HRV channel, external coefficients (mode should have no effect) + # ('HRV', 'radiance', True), + # ('HRV', 'reflectance', True), + ] + ) + def test_calibrate( + self, file_handler, channel, calibration, use_ext_coefs + ): + """Test the calibration.""" + external_coefs = self.external_coefs if use_ext_coefs else {} + expected = self._get_expected( + channel=channel, + calibration=calibration, + calib_mode='NOMINAL', + use_ext_coefs=use_ext_coefs + ) + fh = file_handler + fh.ext_calib_coefs = external_coefs + + dataset_info = { + 'nc_key': channel + } + dataset_id = make_dataid(name=channel, calibration=calibration) + + res = fh.get_dataset(dataset_id, dataset_info) + + # Flip dataset to achieve compatibility with other SEVIRI readers. + # FUTURE: Remove if flipping has been disabled. + res = res.isel(y=slice(None, None, -1)) + + xr.testing.assert_allclose(res, expected)