diff --git a/satpy/etc/readers/mirs.yaml b/satpy/etc/readers/mirs.yaml index b968c849c2..558f8a2254 100644 --- a/satpy/etc/readers/mirs.yaml +++ b/satpy/etc/readers/mirs.yaml @@ -25,4 +25,90 @@ file_types: file_patterns: - 'IMG_SX.{platform_shortname}.D{start_time:%y%j.S%H%M}.E{end_time:%H%M}.B{num}.WE.HR.ORB.nc' -datasets: {} +datasets: + longitude: + name: longitude + file_type: metop_amsu + file_key: Longitude + units: degrees + valid_range: [ -180., 180. ] + standard_name: longitude + latitude: + name: latitude + file_type: metop_amsu + file_key: Latitude + valid_range: [-90., 90.] + units: degrees + standard_name: latitude + rain_rate: + name: RR + description: Rain Rate + file_key: RR + file_type: metop_amsu + units: mm/hr + coordinates: [longitude, latitude] + mask: + name: Sfc_type + file_key: Sfc_type + file_type: metop_amsu + description: Surface Type:0-ocean,1-sea ice,2-land,3-snow + units: "1" + coordinates: [longitude, latitude] + sea_ice: + name: SIce + description: Sea Ice + file_key: SIce + file_type: metop_amsu + units: "%" + coordinates: [longitude, latitude] + snow_cover: + name: Snow + description: Snow Cover + long_name: snow_cover + file_key: Snow + file_type: metop_amsu + units: '1' + coordinates: [longitude, latitude] + total_precipitable_water: + name: TPW + description: Total Precipitable Water + file_key: TPW + file_type: metop_amsu + units: mm + coordinates: [longitude, latitude] + swe: + name: SWE + description: Snow Water Equivalence + file_key: SWE + file_type: metop_amsu + units: cm + coordinates: [longitude, latitude] + cloud_liquid_water: + name: CLW + description: Cloud Liquid Water + file_key: CLW + file_type: metop_amsu + units: mm + coordinates: [longitude, latitude] + skin_temperature: + name: TSkin + description: skin temperature + file_key: TSkin + file_type: metop_amsu + units: K + coordinates: [longitude, latitude] + snow_fall_rate: + name: SFR + description: snow fall rate + file_key: SFR + file_type: metop_amsu + units: mm/hr + coordinates: [longitude, latitude] + bt: + name: BT + file_type: metop_amsu + description: Channel Brightness Temperature for every channel + long_name: Channel Temperature (K) + units: K + valid_range: [0, 50000] + standard_name: brightness_temperature diff --git a/satpy/readers/mirs.py b/satpy/readers/mirs.py index 789ff4f07d..eb8bee2448 100644 --- a/satpy/readers/mirs.py +++ b/satpy/readers/mirs.py @@ -93,16 +93,9 @@ def read_atms_coeff_to_string(fn): def read_atms_limb_correction_coefficients(fn): """Read the limb correction files.""" coeff_str = read_atms_coeff_to_string(fn) - # there should be 22 channels and 96 fov in the coefficient file. (read last line and get values) - n_chn = (coeff_str[-1].replace(" ", " ").split(" ")[1]) - n_fov = (coeff_str[-1].replace(" ", " ").split(" ")[2]) - n_chn = int(n_chn) - n_fov = int(n_fov) - if n_chn < 22: - LOG.warning('Coefficient file has less than 22 channels: %s' % n_chn) - if n_fov < 96: - LOG.warning('Coefficient file has less than 96 fov: %s' % n_fov) - # make it a generator + n_chn = 22 + n_fov = 96 + # make the string a generator coeff_str = (line.strip() for line in coeff_str) all_coeffs = np.zeros((n_chn, n_fov, n_chn), dtype=np.float32) @@ -309,7 +302,7 @@ def _get_coeff_filenames(self): return coeff_fn - def get_metadata(self, ds_info): + def update_metadata(self, ds_info): """Get metadata.""" metadata = {} metadata.update(ds_info) @@ -334,38 +327,71 @@ def _nan_for_dtype(data_arr_dtype): return np.nan @staticmethod - def _scale_data(data_arr, attrs): - # handle scaling - # take special care for integer/category fields - scale_factor = attrs.pop('scale_factor', 1.) - add_offset = attrs.pop('add_offset', 0.) + def _scale_data(data_arr, scale_factor, add_offset): + """Scale data, if needed.""" scaling_needed = not (scale_factor == 1 and add_offset == 0) if scaling_needed: data_arr = data_arr * scale_factor + add_offset - return data_arr, attrs + return data_arr - def _fill_data(self, data_arr, attrs): - try: - global_attr_fill = self.nc.missing_value - except AttributeError: - global_attr_fill = None - fill_value = attrs.pop('_FillValue', global_attr_fill) - - fill_out = self._nan_for_dtype(data_arr.dtype) + def _fill_data(self, data_arr, fill_value, scale_factor, add_offset): + """Fill missing data with NaN.""" if fill_value is not None: + fill_value = self._scale_data(fill_value, scale_factor, add_offset) + fill_out = self._nan_for_dtype(data_arr.dtype) data_arr = data_arr.where(data_arr != fill_value, fill_out) - return data_arr, attrs + return data_arr - def _apply_valid_range(self, data_arr, attrs): - # handle valid_range - valid_range = attrs.pop('valid_range', None) + def _apply_valid_range(self, data_arr, valid_range, scale_factor, add_offset): + """Get and apply valid_range.""" if valid_range is not None: valid_min, valid_max = valid_range + valid_min = self._scale_data(valid_min, scale_factor, add_offset) + valid_max = self._scale_data(valid_max, scale_factor, add_offset) if valid_min is not None and valid_max is not None: data_arr = data_arr.where((data_arr >= valid_min) & (data_arr <= valid_max)) - return data_arr, attrs + return data_arr + + def apply_attributes(self, data, ds_info): + """Combine attributes from file and yaml and apply. + + File attributes should take precedence over yaml if both are present + + """ + try: + global_attr_fill = self.nc.missing_value + except AttributeError: + global_attr_fill = 1.0 + + # let file metadata take precedence over ds_info from yaml, + # but if yaml has more to offer, include it here, but fix + # units. + ds_info.update(data.attrs) + + # special cases + if ds_info['name'] in ["latitude", "longitude"]: + ds_info["standard_name"] = ds_info.get("standard_name", + ds_info['name']) + + # try to assign appropriate units (if "Kelvin" covert to K) + units_convert = {"Kelvin": "K"} + data_unit = ds_info.get('units', None) + ds_info['units'] = units_convert.get(data_unit, data_unit) + + scale = ds_info.pop('scale_factor', 1.0) + offset = ds_info.pop('add_offset', 0.) + fill_value = ds_info.pop("_FillValue", global_attr_fill) + valid_range = ds_info.pop('valid_range', None) + + data = self._scale_data(data, scale, offset) + data = self._fill_data(data, fill_value, scale, offset) + data = self._apply_valid_range(data, valid_range, scale, offset) + + data.attrs = ds_info + + return data, ds_info def get_dataset(self, ds_id, ds_info): """Get datasets.""" @@ -373,23 +399,34 @@ def get_dataset(self, ds_id, ds_info): idx = ds_info['channel_index'] data = self['BT'] data = data.rename(new_name_or_name_dict=ds_info["name"]) + data, ds_info = self.apply_attributes(data, ds_info) if self.sensor.lower() == "atms" and self.limb_correction: sfc_type_mask = self['Sfc_type'] data = limb_correct_atms_bt(data, sfc_type_mask, self._get_coeff_filenames, ds_info) + self.nc = self.nc.merge(data) else: LOG.info("No Limb Correction applied.") data = data[:, :, idx] else: data = self[ds_id['name']] + data, ds_info = self.apply_attributes(data, ds_info) + + data.attrs = self.update_metadata(ds_info) - data.attrs = self.get_metadata(ds_info) return data - def _available_if_this_file_type(self, configured_datasets): + def available_datasets(self, configured_datasets=None): + """Dynamically discover what variables can be loaded from this file. + + See :meth:`satpy.readers.file_handlers.BaseHandler.available_datasets` + for more information. + + """ + handled_vars = set() for is_avail, ds_info in (configured_datasets or []): if is_avail is not None: # some other file handler said it has this dataset @@ -397,7 +434,15 @@ def _available_if_this_file_type(self, configured_datasets): # file handler so let's yield early yield is_avail, ds_info continue - yield self.file_type_matches(ds_info['file_type']), ds_info + + yaml_info = {} + if self.file_type_matches(ds_info['file_type']): + handled_vars.add(ds_info['name']) + yaml_info = ds_info + if ds_info['name'] == 'BT': + yield from self._available_btemp_datasets(yaml_info) + yield True, ds_info + yield from self._available_new_datasets(handled_vars) def _count_channel_repeat_number(self): """Count channel/polarization pair repetition.""" @@ -414,7 +459,7 @@ def _count_channel_repeat_number(self): return chn_total, normals - def _available_btemp_datasets(self): + def _available_btemp_datasets(self, yaml_info): """Create metadata for channel BTs.""" chn_total, normals = self._count_channel_repeat_number() # keep track of current channel count for string description @@ -428,17 +473,17 @@ def _available_btemp_datasets(self): desc_bt = "Channel {} Brightness Temperature at {}GHz {}{}" desc_bt = desc_bt.format(idx, normal_f, normal_p, p_count) - ds_info = { + ds_info = yaml_info.copy() + ds_info.update({ 'file_type': self.filetype_info['file_type'], 'name': new_name, 'description': desc_bt, - 'units': 'K', 'channel_index': idx, 'frequency': "{}GHz".format(normal_f), 'polarization': normal_p, 'dependencies': ('BT', 'Sfc_type'), - 'coordinates': ["longitude", "latitude"] - } + 'coordinates': ['longitude', 'latitude'] + }) yield True, ds_info def _get_ds_info_for_data_arr(self, var_name): @@ -447,9 +492,6 @@ def _get_ds_info_for_data_arr(self, var_name): 'name': var_name, 'coordinates': ["longitude", "latitude"] } - - if var_name in ["longitude", "latitude"]: - ds_info['standard_name'] = var_name return ds_info def _is_2d_yx_data_array(self, data_arr): @@ -457,10 +499,12 @@ def _is_2d_yx_data_array(self, data_arr): has_x_dim = data_arr.dims[1] == "x" return has_y_dim and has_x_dim - def _available_new_datasets(self): + def _available_new_datasets(self, handled_vars): """Metadata for available variables other than BT.""" possible_vars = list(self.nc.items()) + list(self.nc.coords.items()) for var_name, data_arr in possible_vars: + if var_name in handled_vars: + continue if data_arr.ndim != 2: # we don't currently handle non-2D variables continue @@ -471,30 +515,9 @@ def _available_new_datasets(self): ds_info = self._get_ds_info_for_data_arr(var_name) yield True, ds_info - def available_datasets(self, configured_datasets=None): - """Dynamically discover what variables can be loaded from this file. - - See :meth:`satpy.readers.file_handlers.BaseHandler.available_datasets` - for more information. - - """ - yield from self._available_if_this_file_type(configured_datasets) - yield from self._available_new_datasets() - yield from self._available_btemp_datasets() - def __getitem__(self, item): - """Wrap around `self.nc[item]`. - - Some datasets use a 32-bit float scaling factor like the 'x' and 'y' - variables which causes inaccurate unscaled data values. This method - forces the scale factor to a 64-bit float first. - - """ + """Wrap around `self.nc[item]`.""" data = self.nc[item] - attrs = data.attrs.copy() - data, attrs = self._scale_data(data, attrs) - data, attrs = self._fill_data(data, attrs) - data, attrs = self._apply_valid_range(data, attrs) # 'Freq' dimension causes issues in other processing if 'Freq' in data.coords: diff --git a/satpy/tests/reader_tests/test_mirs.py b/satpy/tests/reader_tests/test_mirs.py index 63d48b37fa..df445ab4e7 100644 --- a/satpy/tests/reader_tests/test_mirs.py +++ b/satpy/tests/reader_tests/test_mirs.py @@ -25,12 +25,12 @@ import numpy as np import xarray as xr -AWIPS_FILE = "IMG_SX.M2.D17037.S1601.E1607.B0000001.WE.HR.ORB.nc" +METOP_FILE = "IMG_SX.M2.D17037.S1601.E1607.B0000001.WE.HR.ORB.nc" NPP_MIRS_L2_SWATH = "NPR-MIRS-IMG_v11r6_npp_s201702061601000_e201702061607000_c202012201658410.nc" N20_MIRS_L2_SWATH = "NPR-MIRS-IMG_v11r4_n20_s201702061601000_e201702061607000_c202012201658410.nc" OTHER_MIRS_L2_SWATH = "NPR-MIRS-IMG_v11r4_gpm_s201702061601000_e201702061607000_c202010080001310.nc" -EXAMPLE_FILES = [AWIPS_FILE, NPP_MIRS_L2_SWATH, OTHER_MIRS_L2_SWATH] +EXAMPLE_FILES = [METOP_FILE, NPP_MIRS_L2_SWATH, OTHER_MIRS_L2_SWATH] N_CHANNEL = 3 N_FOV = 96 @@ -50,6 +50,8 @@ DS_IDS = ['RR', 'longitude', 'latitude'] TEST_VARS = ['btemp_88v1', 'btemp_88v2', 'btemp_22h', 'RR', 'Sfc_type'] +DEFAULT_UNITS = {'btemp_88v1': 'K', 'btemp_88v2': 'K', + 'btemp_22h': 'K', 'RR': 'mm/hr', 'Sfc_type': "1"} PLATFORM = {"M2": "metop-a", "NPP": "npp", "GPM": "gpm"} SENSOR = {"m2": "amsu-mhs", "npp": "atms", "gpm": "GPI"} @@ -132,7 +134,7 @@ def _get_datasets_with_attributes(**kwargs): attrs = {'missing_value': -999.} ds = xr.Dataset(ds_vars, attrs=attrs) - + ds = ds.assign_coords({"Freq": FREQ, "Latitude": latitude, "Longitude": longitude}) return ds @@ -173,13 +175,13 @@ def _get_datasets_with_less_attributes(): attrs = {'missing_value': -999.} ds = xr.Dataset(ds_vars, attrs=attrs) - + ds = ds.assign_coords({"Freq": FREQ, "Latitude": latitude, "Longitude": longitude}) return ds def fake_open_dataset(filename, **kwargs): """Create a Dataset similar to reading an actual file with xarray.open_dataset.""" - if filename == AWIPS_FILE: + if filename == METOP_FILE: return _get_datasets_with_less_attributes() return _get_datasets_with_attributes() @@ -197,7 +199,7 @@ def setup_method(self): @pytest.mark.parametrize( ("filenames", "expected_loadables"), [ - ([AWIPS_FILE], 1), + ([METOP_FILE], 1), ([NPP_MIRS_L2_SWATH], 1), ([OTHER_MIRS_L2_SWATH], 1), ] @@ -217,7 +219,7 @@ def test_reader_creation(self, filenames, expected_loadables): @pytest.mark.parametrize( ("filenames", "expected_datasets"), [ - ([AWIPS_FILE], DS_IDS), + ([METOP_FILE], DS_IDS), ([NPP_MIRS_L2_SWATH], DS_IDS), ([OTHER_MIRS_L2_SWATH], DS_IDS), ] @@ -254,6 +256,11 @@ def _check_valid_range(data_arr, test_valid_range): assert data_arr.data.min() >= test_valid_range[0] assert data_arr.data.max() <= test_valid_range[1] + @staticmethod + def _check_fill_value(data_arr, test_fill_value): + assert '_FillValue' not in data_arr.attrs + assert test_fill_value not in data_arr.data + @staticmethod def _check_attrs(data_arr, platform_name): attrs = data_arr.attrs @@ -266,12 +273,7 @@ def _check_attrs(data_arr, platform_name): @pytest.mark.parametrize( ("filenames", "loadable_ids", "platform_name"), [ - ([AWIPS_FILE], TEST_VARS, "metop-a"), - ([NPP_MIRS_L2_SWATH], TEST_VARS, "npp"), - ([N20_MIRS_L2_SWATH], TEST_VARS, "noaa-20"), - ([OTHER_MIRS_L2_SWATH], TEST_VARS, "gpm"), - - ([AWIPS_FILE], TEST_VARS, "metop-a"), + ([METOP_FILE], TEST_VARS, "metop-a"), ([NPP_MIRS_L2_SWATH], TEST_VARS, "npp"), ([N20_MIRS_L2_SWATH], TEST_VARS, "noaa-20"), ([OTHER_MIRS_L2_SWATH], TEST_VARS, "gpm"), @@ -306,9 +308,13 @@ def test_basic_load(self, filenames, loadable_ids, if "valid_range" in input_fake_data.attrs: valid_range = input_fake_data.attrs["valid_range"] self._check_valid_range(data_arr, valid_range) + if "_FillValue" in input_fake_data.attrs: + fill_value = input_fake_data.attrs["_FillValue"] + self._check_fill_value(data_arr, fill_value) sensor = data_arr.attrs['sensor'] if reader_kw.get('limb_correction', True) and sensor == 'atms': fd.assert_called() else: fd.assert_not_called() + assert data_arr.attrs['units'] == DEFAULT_UNITS[var_name]