diff --git a/satpy/etc/writers/awips_tiled.yaml b/satpy/etc/writers/awips_tiled.yaml index 9e41e1c454..58a79fa9ae 100644 --- a/satpy/etc/writers/awips_tiled.yaml +++ b/satpy/etc/writers/awips_tiled.yaml @@ -66,9 +66,12 @@ templates: # value: "${ORGANIZATION}" awips_id: {} # value: "{awips_id}" # special variable created by awips_tiled.py + physical_element: {} +# value: "{physical_element}" #special variable created by awips_tiled.py satellite_id: value: "{platform_name!u}-{sensor!u}" sector_id: {} # special handler in awips_tiled.py + coordinates: x: attributes: @@ -117,6 +120,7 @@ templates: reader: clavrx var_name: data attributes: + units: {} physical_element: value: 'CLAVR-x {name}' clavrx_cloud_type: @@ -126,11 +130,13 @@ templates: attributes: physical_element: raw_value: CLAVR-x Cloud Type + units: {} clavrx_cld_temp_acha: reader: clavrx name: cld_temp_acha var_name: data attributes: + units: {} physical_element: raw_value: CLAVR-x Cloud Top Temperature (ACHA) clavrx_cld_height_acha: @@ -138,6 +144,7 @@ templates: name: cld_height_acha var_name: data attributes: + units: {} physical_element: raw_value: CLAVR-x Cloud Top Height (ACHA) clavrx_cloud_phase: @@ -145,6 +152,7 @@ templates: name: cloud_phase var_name: data attributes: + units: {} physical_element: raw_value: CLAVR-x Cloud Phase clavrx_cld_opd_dcomp: @@ -152,6 +160,7 @@ templates: name: cld_opd_dcomp var_name: data attributes: + units: {} physical_element: raw_value: CLAVR-x Cloud Optical Depth (dcomp) clavrx_clld_opd_nlcomp: @@ -159,6 +168,7 @@ templates: name: cloud_opd_nlcomp var_name: data attributes: + units: {} physical_element: raw_value: CLAVR-x Cloud Optical Depth (nlcomp) clavrx_cld_reff_dcomp: @@ -166,6 +176,7 @@ templates: name: cld_reff_dcomp var_name: data attributes: + units: {} physical_element: raw_value: CLAVR-x Cloud Effective Radius (dcomp) clavrx_cld_reff_nlcomp: @@ -173,6 +184,7 @@ templates: name: cld_reff_nlcomp var_name: data attributes: + units: {} physical_element: raw_value: CLAVR-x Cloud Effective Radius (nlcomp) clavrx_cld_emiss_acha: @@ -180,6 +192,7 @@ templates: name: cld_emiss_acha var_name: data attributes: + units: {} physical_element: raw_value: CLAVR-x Cloud Emissivity (ACHA) clavrx_refl_lunar_dnb_nom: @@ -187,6 +200,7 @@ templates: name: refl_lunar_dnb_nom var_name: data attributes: + units: {} physical_element: raw_value: CLAVR-x Cloud Lunar Reflectance clavrx_rain_rate: @@ -194,6 +208,7 @@ templates: name: rain_rate var_name: data attributes: + units: {} physical_element: raw_value: CLAVR-x Rain Rate diff --git a/satpy/readers/clavrx.py b/satpy/readers/clavrx.py index 5f26d617aa..b44fdd25a5 100644 --- a/satpy/readers/clavrx.py +++ b/satpy/readers/clavrx.py @@ -96,7 +96,7 @@ class _CLAVRxHelper: def _remove_attributes(attrs: dict) -> dict: """Remove attributes that described data before scaling.""" old_attrs = ['unscaled_missing', 'SCALED_MIN', 'SCALED_MAX', - 'SCALED_MISSING', 'actual_missing'] + 'SCALED_MISSING'] for attr_key in old_attrs: attrs.pop(attr_key, None) @@ -105,31 +105,39 @@ def _remove_attributes(attrs: dict) -> dict: @staticmethod def _scale_data(data_arr: xr.DataArray, scale_factor: float, add_offset: float) -> xr.DataArray: """Scale data, if needed.""" - scaling_needed = not (scale_factor == 1 and add_offset == 0) + scaling_needed = not (scale_factor == 1.0 and add_offset == 0.0) if scaling_needed: data_arr = data_arr * scale_factor + add_offset return data_arr @staticmethod - def _get_data(data: xr.DataArray, dataset_id: dict, ds_info: dict) -> xr.DataArray: + def _get_data(data: xr.DataArray, dataset_id: dict) -> xr.DataArray: """Get a dataset.""" if dataset_id.get('resolution'): data.attrs['resolution'] = dataset_id['resolution'] attrs = data.attrs.copy() - fill = attrs.pop('_FillValue', None) - factor = attrs.pop('scale_factor', 1.0) - offset = attrs.pop('add_offset', 0.0) - valid_range = attrs.pop('valid_range', None) - data = data.where(data != fill) - data = _CLAVRxHelper._scale_data(data, factor, offset) - - if valid_range is not None: + fill = attrs.get('_FillValue') + factor = attrs.pop('scale_factor', (np.ones(1, dtype=data.dtype))[0]) + offset = attrs.pop('add_offset', (np.zeros(1, dtype=data.dtype))[0]) + valid_range = attrs.get('valid_range', [None]) + if isinstance(valid_range, np.ndarray): + attrs["valid_range"] = valid_range.tolist() + + flags = not data.attrs.get("SCALED", 1) and any(data.attrs.get("flag_values", [None])) + if not flags: + data = data.where(data != fill) + data = _CLAVRxHelper._scale_data(data, factor, offset) + # don't need _FillValue if it has been applied. + attrs.pop('_FillValue', None) + fill = _CLAVRxHelper._scale_data(fill, factor, offset) + + if all(valid_range): valid_min = _CLAVRxHelper._scale_data(valid_range[0], factor, offset) valid_max = _CLAVRxHelper._scale_data(valid_range[1], factor, offset) data = data.where((data >= valid_min) & (data <= valid_max)) - data.attrs['valid_min'], data.attrs['valid_max'] = valid_min, valid_max + attrs['valid_range'] = [valid_min, valid_max] data.attrs = _CLAVRxHelper._remove_attributes(attrs) @@ -237,6 +245,32 @@ def _read_axi_fixed_grid(filename: str, l1b_attr) -> geometry.AreaDefinition: return area + @staticmethod + def get_metadata(sensor, platform, attrs: dict, ds_info: dict) -> dict: + """Get metadata.""" + i = {} + i.update(attrs) + i.update(ds_info) + + flag_meanings = i.get('flag_meanings', None) + if not i.get('SCALED', 1) and not flag_meanings: + i['flag_meanings'] = '' + i.setdefault('flag_values', [None]) + u = i.get('units') + if u in CF_UNITS: + # CF compliance + i['units'] = CF_UNITS[u] + if u.lower() == "none": + i['units'] = "1" + i['sensor'] = sensor + i['platform_name'] = platform + rps = _get_rows_per_scan(sensor) + if rps: + i['rows_per_scan'] = rps + i['reader'] = 'clavrx' + + return i + class CLAVRXHDF4FileHandler(HDF4FileHandler, _CLAVRxHelper): """A file handler for CLAVRx files.""" @@ -257,36 +291,13 @@ def end_time(self): """Get the end time.""" return self.filename_info.get('end_time', self.start_time) - def get_metadata(self, attrs: dict, ds_info: dict) -> dict: - """Get metadata.""" - i = {} - i.update(attrs) - i.update(ds_info) - - flag_meanings = i.get('flag_meanings', None) - if not i.get('SCALED', 1) and not flag_meanings: - i['flag_meanings'] = '' - i.setdefault('flag_values', [None]) - u = i.get('units') - if u in CF_UNITS: - # CF compliance - i['units'] = CF_UNITS[u] - - i['sensor'] = self.sensor - i['platform'] = i['platform_name'] = self.platform - rps = _get_rows_per_scan(self.sensor) - if rps: - i['rows_per_scan'] = rps - i['reader'] = 'clavrx' - - return i - def get_dataset(self, dataset_id, ds_info): """Get a dataset.""" var_name = ds_info.get('file_key', dataset_id['name']) data = self[var_name] - data = _CLAVRxHelper._get_data(data, dataset_id, ds_info) - data.attrs = self.get_metadata(data.attrs, ds_info) + data = _CLAVRxHelper._get_data(data, dataset_id) + data.attrs = _CLAVRxHelper.get_metadata(self.sensor, self.platform, + data.attrs, ds_info) return data def get_nadir_resolution(self, sensor): @@ -380,8 +391,7 @@ def __init__(self, filename, filename_info, filetype_info): decode_cf=True, mask_and_scale=False, decode_coords=True, - chunks={'pixel_elements_along_scan_direction': CHUNK_SIZE, - 'scan_lines_along_track_direction': CHUNK_SIZE}) + chunks=CHUNK_SIZE) # y,x is used in satpy, bands rather than channel using in xrimage self.nc = self.nc.rename_dims({'scan_lines_along_track_direction': "y", 'pixel_elements_along_scan_direction': "x"}) @@ -389,12 +399,16 @@ def __init__(self, filename, filename_info, filetype_info): self.platform = _get_platform( self.filename_info.get('platform_shortname', None)) self.sensor = self.nc.attrs.get('sensor', None) + # coordinates need scaling and valid_range (mask_and_scale won't work on valid_range) + self.nc.coords["latitude"] = _CLAVRxHelper._get_data(self.nc.coords["latitude"], + {"name": "latitude"}) + self.nc.coords["longitude"] = _CLAVRxHelper._get_data(self.nc.coords["longitude"], + {"name": "longitude"}) def _get_ds_info_for_data_arr(self, var_name): ds_info = { 'file_type': self.filetype_info['file_type'], 'name': var_name, - 'coordinates': ["longitude", "latitude"] } return ds_info @@ -453,38 +467,13 @@ def get_area_def(self, key): l1b_att = str(self.nc.attrs.get('L1B', None)) return _CLAVRxHelper._read_axi_fixed_grid(self.filename, l1b_att) - def get_metadata(self, attrs, ds_info): - """Get metadata.""" - i = {} - i.update(attrs) - i.update(ds_info) - - flag_meanings = i.get('flag_meanings', None) - if not i.get('SCALED', 1) and not flag_meanings: - i['flag_meanings'] = '' - i.setdefault('flag_values', [None]) - - u = i.get('units') - if u in CF_UNITS: - # CF compliance - i['units'] = CF_UNITS[u] - - i['sensor'] = self.sensor - i['platform'] = i['platform_name'] = self.platform - rps = _get_rows_per_scan(self.sensor) - if rps: - i['rows_per_scan'] = rps - i['reader'] = 'clavrx' - - return i - def get_dataset(self, dataset_id, ds_info): """Get a dataset.""" var_name = ds_info.get('name', dataset_id['name']) data = self[var_name] - data = _CLAVRxHelper._get_data(data, dataset_id, ds_info) - data.attrs = self.get_metadata(data.attrs, ds_info) - + data = _CLAVRxHelper._get_data(data, dataset_id) + data.attrs = _CLAVRxHelper.get_metadata(self.sensor, self.platform, + data.attrs, ds_info) return data def __getitem__(self, item): diff --git a/satpy/tests/reader_tests/test_clavrx.py b/satpy/tests/reader_tests/test_clavrx.py index 87bea6a0f9..6eb8c5593e 100644 --- a/satpy/tests/reader_tests/test_clavrx.py +++ b/satpy/tests/reader_tests/test_clavrx.py @@ -98,7 +98,7 @@ def get_test_content(self, filename, filename_info, filetype_info): '_FillValue': -128, 'flag_meanings': 'clear water supercooled mixed ice unknown', 'flag_values': [0, 1, 2, 3, 4, 5], - 'units': '1', + 'units': 'none', }) file_content['variable3/shape'] = DEFAULT_FILE_SHAPE @@ -217,9 +217,8 @@ def test_load_all(self): 'variable3']) self.assertEqual(len(datasets), 3) for v in datasets.values(): - assert 'calibration' not in v.attrs self.assertEqual(v.attrs['units'], '1') - self.assertEqual(v.attrs['platform'], 'npp') + self.assertEqual(v.attrs['platform_name'], 'npp') self.assertEqual(v.attrs['sensor'], 'viirs') self.assertIsInstance(v.attrs['area'], SwathDefinition) self.assertEqual(v.attrs['area'].lons.attrs['rows_per_scan'], 16) @@ -373,9 +372,13 @@ def test_load_all_old_donor(self): datasets = r.load(['variable1', 'variable2', 'variable3']) self.assertEqual(len(datasets), 3) for v in datasets.values(): - assert 'calibration' not in v.attrs + self.assertNotIn('calibration', v.attrs) self.assertEqual(v.attrs['units'], '1') self.assertIsInstance(v.attrs['area'], AreaDefinition) + if v.attrs["name"] == 'variable1': + self.assertIsInstance(v.attrs["valid_range"], list) + else: + self.assertNotIn('valid_range', v.attrs) self.assertIsNotNone(datasets['variable3'].attrs.get('flag_meanings')) def test_load_all_new_donor(self): @@ -406,10 +409,10 @@ def test_load_all_new_donor(self): datasets = r.load(['variable1', 'variable2', 'variable3']) self.assertEqual(len(datasets), 3) for v in datasets.values(): - assert 'calibration' not in v.attrs + self.assertNotIn('calibration', v.attrs) self.assertEqual(v.attrs['units'], '1') self.assertIsInstance(v.attrs['area'], AreaDefinition) self.assertTrue(v.attrs['area'].is_geostationary) - self.assertEqual(v.attrs['platform'], 'himawari8') + self.assertEqual(v.attrs['platform_name'], 'himawari8') self.assertEqual(v.attrs['sensor'], 'ahi') self.assertIsNotNone(datasets['variable3'].attrs.get('flag_meanings')) diff --git a/satpy/tests/reader_tests/test_clavrx_nc.py b/satpy/tests/reader_tests/test_clavrx_nc.py index c705561d54..3dd5d8befe 100644 --- a/satpy/tests/reader_tests/test_clavrx_nc.py +++ b/satpy/tests/reader_tests/test_clavrx_nc.py @@ -73,7 +73,7 @@ def fake_test_content(filename, **kwargs): 'scale_factor': 1., 'add_offset': 0., 'units': '1', - 'valid_range': (-32767, 32767), + 'valid_range': [-32767, 32767], }) # data with fill values @@ -84,7 +84,7 @@ def fake_test_content(filename, **kwargs): 'scale_factor': 1., 'add_offset': 0., 'units': '1', - 'valid_range': (-32767, 32767), + 'valid_range': [-32767, 32767], }) variable2 = variable2.where(variable2 % 2 != 0) @@ -188,7 +188,9 @@ def test_load_all_new_donor(self, filenames, loadable_ids): assert 'calibration' not in v.attrs assert v.attrs['units'] == '1' assert isinstance(v.attrs['area'], AreaDefinition) - assert v.attrs['platform'] == 'himawari8' + assert v.attrs['platform_name'] == 'himawari8' assert v.attrs['sensor'] == 'AHI' assert 'rows_per_scan' not in v.coords.get('longitude').attrs + if v.attrs["name"] in ["variable1", "variable2"]: + assert isinstance(v.attrs["valid_range"], list) assert (datasets['variable3'].attrs.get('flag_meanings')) is not None diff --git a/satpy/tests/writer_tests/test_awips_tiled.py b/satpy/tests/writer_tests/test_awips_tiled.py index 83b8e04be9..950740c8b0 100644 --- a/satpy/tests/writer_tests/test_awips_tiled.py +++ b/satpy/tests/writer_tests/test_awips_tiled.py @@ -22,6 +22,7 @@ import numpy as np import dask.array as da +from pyproj import CRS import pytest @@ -37,8 +38,33 @@ def _check_production_location(ds): assert len(ds.attrs[prod_loc_name]) == 31 -def check_required_common_attributes(ds): +def check_required_properties(unmasked_ds, masked_ds): + """Check various aspects of coordinates and attributes for correctness.""" + _check_scaled_x_coordinate_variable(unmasked_ds, masked_ds) + _check_scaled_y_coordinate_variable(unmasked_ds, masked_ds) + _check_required_common_attributes(unmasked_ds) + + +def _check_required_common_attributes(ds): """Check common properties of the created AWIPS tiles for validity.""" + for attr_name in ('tile_row_offset', 'tile_column_offset', + 'product_tile_height', 'product_tile_width', + 'number_product_tiles', + 'product_rows', 'product_columns'): + assert attr_name in ds.attrs + _check_production_location(ds) + + for data_arr in ds.data_vars.values(): + if data_arr.ndim == 0: + # grid mapping variable + assert 'grid_mapping_name' in data_arr.attrs + continue + assert data_arr.encoding.get('zlib', False) + assert 'grid_mapping' in data_arr.attrs + assert data_arr.attrs['grid_mapping'] in ds + + +def _check_scaled_x_coordinate_variable(ds, masked_ds): assert 'x' in ds.coords x_coord = ds.coords['x'] np.testing.assert_equal(np.diff(x_coord), 1) @@ -46,8 +72,14 @@ def check_required_common_attributes(ds): assert x_attrs.get('standard_name') == 'projection_x_coordinate' assert x_attrs.get('units') == 'meters' assert 'scale_factor' in x_attrs + assert x_attrs['scale_factor'] > 0 assert 'add_offset' in x_attrs + unscaled_x = masked_ds.coords['x'].values + assert (np.diff(unscaled_x) > 0).all() + + +def _check_scaled_y_coordinate_variable(ds, masked_ds): assert 'y' in ds.coords y_coord = ds.coords['y'] np.testing.assert_equal(np.diff(y_coord), 1) @@ -55,23 +87,11 @@ def check_required_common_attributes(ds): assert y_attrs.get('standard_name') == 'projection_y_coordinate' assert y_attrs.get('units') == 'meters' assert 'scale_factor' in y_attrs + assert y_attrs['scale_factor'] < 0 assert 'add_offset' in y_attrs - for attr_name in ('tile_row_offset', 'tile_column_offset', - 'product_tile_height', 'product_tile_width', - 'number_product_tiles', - 'product_rows', 'product_columns'): - assert attr_name in ds.attrs - _check_production_location(ds) - - for data_arr in ds.data_vars.values(): - if data_arr.ndim == 0: - # grid mapping variable - assert 'grid_mapping_name' in data_arr.attrs - continue - assert data_arr.encoding.get('zlib', False) - assert 'grid_mapping' in data_arr.attrs - assert data_arr.attrs['grid_mapping'] in ds + unscaled_y = masked_ds.coords['y'].values + assert (np.diff(unscaled_y) < 0).all() class TestAWIPSTiledWriter: @@ -81,6 +101,8 @@ def setup_method(self): """Create temporary directory to save files to.""" import tempfile self.base_dir = tempfile.mkdtemp() + self.start_time = datetime(2018, 1, 1, 12, 0, 0) + self.end_time = self.start_time + timedelta(minutes=20) def teardown_method(self): """Remove the temporary directory created for a test.""" @@ -95,32 +117,39 @@ def test_init(self): from satpy.writers.awips_tiled import AWIPSTiledWriter AWIPSTiledWriter(base_dir=self.base_dir) - def _get_test_lcc_data(self): - from xarray import DataArray + def _get_test_area(self, shape=(200, 100), crs=None, extents=None): from pyresample.geometry import AreaDefinition - from pyresample.utils import proj4_str_to_dict + if crs is None: + crs = CRS('+proj=lcc +datum=WGS84 +ellps=WGS84 +lon_0=-95. +lat_0=25 +lat_1=25 +units=m +no_defs') + if extents is None: + extents = (-1000., -1500., 1000., 1500.) area_def = AreaDefinition( 'test', 'test', 'test', - proj4_str_to_dict('+proj=lcc +datum=WGS84 +ellps=WGS84 +lon_0=-95. ' - '+lat_0=25 +lat_1=25 +units=m +no_defs'), - 100, - 200, - (-1000., -1500., 1000., 1500.), + crs, + shape[1], + shape[0], + extents, ) - now = datetime(2018, 1, 1, 12, 0, 0) - data = np.linspace(0., 1., 20000, dtype=np.float32).reshape((200, 100)) + return area_def + + def _get_test_data(self, shape=(200, 100), chunks=50): + data = np.linspace(0., 1., shape[0] * shape[1], dtype=np.float32).reshape(shape) + return da.from_array(data, chunks=chunks) + + def _get_test_lcc_data(self, dask_arr, area_def): + from xarray import DataArray ds = DataArray( - da.from_array(data, chunks=50), + dask_arr, attrs=dict( name='test_ds', platform_name='PLAT', sensor='SENSOR', units='1', area=area_def, - start_time=now, - end_time=now + timedelta(minutes=20)) + start_time=self.start_time, + end_time=self.end_time) ) return ds @@ -130,7 +159,9 @@ def test_basic_numbered_1_tile(self, use_save_dataset): """Test creating a single numbered tile.""" import xarray as xr from satpy.writers.awips_tiled import AWIPSTiledWriter - input_data_arr = self._get_test_lcc_data() + data = self._get_test_data() + area_def = self._get_test_area() + input_data_arr = self._get_test_lcc_data(data, area_def) w = AWIPSTiledWriter(base_dir=self.base_dir, compress=True) if use_save_dataset: w.save_dataset(input_data_arr, sector_id='TEST', source_name='TESTS') @@ -141,9 +172,9 @@ def test_basic_numbered_1_tile(self, use_save_dataset): assert len(all_files) == 1 assert os.path.basename(all_files[0]) == 'TESTS_AII_PLAT_SENSOR_test_ds_TEST_T001_20180101_1200.nc' for fn in all_files: - output_ds = xr.open_dataset(fn, mask_and_scale=False) - check_required_common_attributes(output_ds) + unmasked_ds = xr.open_dataset(fn, mask_and_scale=False) output_ds = xr.open_dataset(fn, mask_and_scale=True) + check_required_properties(unmasked_ds, output_ds) scale_factor = output_ds['data'].encoding['scale_factor'] np.testing.assert_allclose(input_data_arr.values, output_ds['data'].data, atol=scale_factor / 2) @@ -162,7 +193,9 @@ def test_basic_numbered_tiles(self, tile_count, tile_size): import dask from satpy.writers.awips_tiled import AWIPSTiledWriter from satpy.tests.utils import CustomScheduler - input_data_arr = self._get_test_lcc_data() + data = self._get_test_data() + area_def = self._get_test_area() + input_data_arr = self._get_test_lcc_data(data, area_def) w = AWIPSTiledWriter(base_dir=self.base_dir, compress=True) save_kwargs = dict( sector_id='TEST', @@ -184,88 +217,71 @@ def test_basic_numbered_tiles(self, tile_count, tile_size): expected_num_files = 0 if should_error else 9 assert len(all_files) == expected_num_files for fn in all_files: - ds = xr.open_dataset(fn, mask_and_scale=False) - check_required_common_attributes(ds) - assert ds.attrs['my_global'] == 'TEST' - assert ds.attrs['sector_id'] == 'TEST' + unmasked_ds = xr.open_dataset(fn, mask_and_scale=False) + masked_ds = xr.open_dataset(fn, mask_and_scale=True) + check_required_properties(unmasked_ds, masked_ds) + assert unmasked_ds.attrs['my_global'] == 'TEST' + assert unmasked_ds.attrs['sector_id'] == 'TEST' + assert 'physical_element' in unmasked_ds.attrs stime = input_data_arr.attrs['start_time'] - assert ds.attrs['start_date_time'] == stime.strftime('%Y-%m-%dT%H:%M:%S') + assert unmasked_ds.attrs['start_date_time'] == stime.strftime('%Y-%m-%dT%H:%M:%S') def test_basic_lettered_tiles(self): """Test creating a lettered grid.""" import xarray as xr from satpy.writers.awips_tiled import AWIPSTiledWriter - from xarray import DataArray - from pyresample.geometry import AreaDefinition - from pyresample.utils import proj4_str_to_dict w = AWIPSTiledWriter(base_dir=self.base_dir, compress=True) - area_def = AreaDefinition( - 'test', - 'test', - 'test', - proj4_str_to_dict('+proj=lcc +datum=WGS84 +ellps=WGS84 +lon_0=-95. ' - '+lat_0=25 +lat_1=25 +units=m +no_defs'), - 1000, - 2000, - (-1000000., -1500000., 1000000., 1500000.), - ) - now = datetime(2018, 1, 1, 12, 0, 0) - ds = DataArray( - da.from_array(np.linspace(0., 1., 2000000, dtype=np.float32).reshape((2000, 1000)), chunks=500), - attrs=dict( - name='test_ds', - platform_name='PLAT', - sensor='SENSOR', - units='1', - area=area_def, - start_time=now, - end_time=now + timedelta(minutes=20)) - ) + data = self._get_test_data(shape=(2000, 1000), chunks=500) + area_def = self._get_test_area(shape=(2000, 1000), + extents=(-1000000., -1500000., 1000000., 1500000.)) + ds = self._get_test_lcc_data(data, area_def) # tile_count should be ignored since we specified lettered_grid w.save_datasets([ds], sector_id='LCC', source_name="TESTS", tile_count=(3, 3), lettered_grid=True) all_files = glob(os.path.join(self.base_dir, 'TESTS_AII*.nc')) assert len(all_files) == 16 for fn in all_files: - ds = xr.open_dataset(fn, mask_and_scale=False) - check_required_common_attributes(ds) - assert ds.attrs['start_date_time'] == now.strftime('%Y-%m-%dT%H:%M:%S') + unmasked_ds = xr.open_dataset(fn, mask_and_scale=False) + masked_ds = xr.open_dataset(fn, mask_and_scale=True) + check_required_properties(unmasked_ds, masked_ds) + assert masked_ds.attrs['start_date_time'] == self.start_time.strftime('%Y-%m-%dT%H:%M:%S') + + def test_basic_lettered_tiles_diff_projection(self): + """Test creating a lettered grid from data with differing projection..""" + import xarray as xr + from satpy.writers.awips_tiled import AWIPSTiledWriter + w = AWIPSTiledWriter(base_dir=self.base_dir, compress=True) + crs = CRS("+proj=lcc +datum=WGS84 +ellps=WGS84 +lon_0=-95. +lat_0=45 +lat_1=45 +units=m +no_defs") + data = self._get_test_data(shape=(2000, 1000), chunks=500) + area_def = self._get_test_area(shape=(2000, 1000), crs=crs, + extents=(-1000000., -1500000., 1000000., 1500000.)) + ds = self._get_test_lcc_data(data, area_def) + # tile_count should be ignored since we specified lettered_grid + w.save_datasets([ds], sector_id='LCC', source_name="TESTS", tile_count=(3, 3), lettered_grid=True) + all_files = sorted(glob(os.path.join(self.base_dir, 'TESTS_AII*.nc'))) + assert len(all_files) == 24 + assert "TC02" in all_files[0] # the first tile should be TC02 + for fn in all_files: + unmasked_ds = xr.open_dataset(fn, mask_and_scale=False) + masked_ds = xr.open_dataset(fn, mask_and_scale=True) + check_required_properties(unmasked_ds, masked_ds) + assert masked_ds.attrs['start_date_time'] == self.start_time.strftime('%Y-%m-%dT%H:%M:%S') def test_lettered_tiles_update_existing(self): """Test updating lettered tiles with additional data.""" import shutil import xarray as xr from satpy.writers.awips_tiled import AWIPSTiledWriter - from xarray import DataArray - from pyresample.geometry import AreaDefinition - from pyresample.utils import proj4_str_to_dict import dask first_base_dir = os.path.join(self.base_dir, 'first') w = AWIPSTiledWriter(base_dir=first_base_dir, compress=True) - area_def = AreaDefinition( - 'test', - 'test', - 'test', - proj4_str_to_dict('+proj=lcc +datum=WGS84 +ellps=WGS84 +lon_0=-95. ' - '+lat_0=25 +lat_1=25 +units=m +no_defs'), - 1000, - 2000, - (-1000000., -1500000., 1000000., 1500000.), - ) - now = datetime(2018, 1, 1, 12, 0, 0) - data = np.linspace(0., 1., 2000000, dtype=np.float32).reshape((2000, 1000)) + shape = (2000, 1000) + data = np.linspace(0., 1., shape[0] * shape[1], dtype=np.float32).reshape(shape) # pixels to be filled in later data[:, -200:] = np.nan - ds = DataArray( - da.from_array(data, chunks=500), - attrs=dict( - name='test_ds', - platform_name='PLAT', - sensor='SENSOR', - units='1', - area=area_def, - start_time=now, - end_time=now + timedelta(minutes=20)) - ) + data = da.from_array(data, chunks=500) + area_def = self._get_test_area(shape=(2000, 1000), + extents=(-1000000., -1500000., 1000000., 1500000.)) + ds = self._get_test_lcc_data(data, area_def) # tile_count should be ignored since we specified lettered_grid w.save_datasets([ds], sector_id='LCC', source_name="TESTS", tile_count=(3, 3), lettered_grid=True) all_files = sorted(glob(os.path.join(first_base_dir, 'TESTS_AII*.nc'))) @@ -280,32 +296,15 @@ def test_lettered_tiles_update_existing(self): # Second writing/updating # Area is about 100 pixels to the right - area_def2 = AreaDefinition( - 'test', - 'test', - 'test', - proj4_str_to_dict('+proj=lcc +datum=WGS84 +ellps=WGS84 +lon_0=-95. ' - '+lat_0=25 +lat_1=25 +units=m +no_defs'), - 1000, - 2000, - (-800000., -1500000., 1200000., 1500000.), - ) + area_def2 = self._get_test_area(shape=(2000, 1000), + extents=(-800000., -1500000., 1200000., 1500000.)) data2 = np.linspace(0., 1., 2000000, dtype=np.float32).reshape((2000, 1000)) # a gap at the beginning where old values remain data2[:, :200] = np.nan # a gap at the end where old values remain data2[:, -400:-300] = np.nan - ds2 = DataArray( - da.from_array(data2, chunks=500), - attrs=dict( - name='test_ds', - platform_name='PLAT', - sensor='SENSOR', - units='1', - area=area_def2, - start_time=now, - end_time=now + timedelta(minutes=20)) - ) + data2 = da.from_array(data2, chunks=500) + ds2 = self._get_test_lcc_data(data2, area_def2) w = AWIPSTiledWriter(base_dir=second_base_dir, compress=True) # HACK: The _copy_to_existing function hangs when opening the output # file multiple times...sometimes. If we limit dask to one worker @@ -341,71 +340,31 @@ def test_lettered_tiles_sector_ref(self): """Test creating a lettered grid using the sector as reference.""" import xarray as xr from satpy.writers.awips_tiled import AWIPSTiledWriter - from xarray import DataArray - from pyresample.geometry import AreaDefinition - from pyresample.utils import proj4_str_to_dict w = AWIPSTiledWriter(base_dir=self.base_dir, compress=True) - area_def = AreaDefinition( - 'test', - 'test', - 'test', - proj4_str_to_dict('+proj=lcc +datum=WGS84 +ellps=WGS84 +lon_0=-95. ' - '+lat_0=25 +lat_1=25 +units=m +no_defs'), - 1000, - 2000, - (-1000000., -1500000., 1000000., 1500000.), - ) - now = datetime(2018, 1, 1, 12, 0, 0) - ds = DataArray( - da.from_array(np.linspace(0., 1., 2000000, dtype=np.float32).reshape((2000, 1000)), chunks=500), - attrs=dict( - name='test_ds', - platform_name='PLAT', - sensor='SENSOR', - units='1', - area=area_def, - start_time=now, - end_time=now + timedelta(minutes=20)) - ) + data = self._get_test_data(shape=(2000, 1000), chunks=500) + area_def = self._get_test_area(shape=(2000, 1000), + extents=(-1000000., -1500000., 1000000., 1500000.)) + ds = self._get_test_lcc_data(data, area_def) w.save_datasets([ds], sector_id='LCC', source_name="TESTS", lettered_grid=True, use_sector_reference=True, use_end_time=True) all_files = glob(os.path.join(self.base_dir, 'TESTS_AII*.nc')) assert len(all_files) == 16 for fn in all_files: - ds = xr.open_dataset(fn, mask_and_scale=False) - check_required_common_attributes(ds) - assert ds.attrs['start_date_time'] == (now + timedelta(minutes=20)).strftime('%Y-%m-%dT%H:%M:%S') + unmasked_ds = xr.open_dataset(fn, mask_and_scale=False) + masked_ds = xr.open_dataset(fn, mask_and_scale=True) + check_required_properties(unmasked_ds, masked_ds) + expected_start = (self.start_time + timedelta(minutes=20)).strftime('%Y-%m-%dT%H:%M:%S') + assert masked_ds.attrs['start_date_time'] == expected_start def test_lettered_tiles_no_fit(self): """Test creating a lettered grid with no data overlapping the grid.""" from satpy.writers.awips_tiled import AWIPSTiledWriter - from xarray import DataArray - from pyresample.geometry import AreaDefinition - from pyresample.utils import proj4_str_to_dict w = AWIPSTiledWriter(base_dir=self.base_dir, compress=True) - area_def = AreaDefinition( - 'test', - 'test', - 'test', - proj4_str_to_dict('+proj=lcc +datum=WGS84 +ellps=WGS84 +lon_0=-95. ' - '+lat_0=25 +lat_1=25 +units=m +no_defs'), - 1000, - 2000, - (4000000., 5000000., 5000000., 6000000.), - ) - now = datetime(2018, 1, 1, 12, 0, 0) - ds = DataArray( - da.from_array(np.linspace(0., 1., 2000000, dtype=np.float32).reshape((2000, 1000)), chunks=500), - attrs=dict( - name='test_ds', - platform_name='PLAT', - sensor='SENSOR', - units='1', - area=area_def, - start_time=now, - end_time=now + timedelta(minutes=20)) - ) + data = self._get_test_data(shape=(2000, 1000), chunks=500) + area_def = self._get_test_area(shape=(2000, 1000), + extents=(4000000., 5000000., 5000000., 6000000.)) + ds = self._get_test_lcc_data(data, area_def) w.save_datasets([ds], sector_id='LCC', source_name="TESTS", tile_count=(3, 3), lettered_grid=True) # No files created all_files = glob(os.path.join(self.base_dir, 'TESTS_AII*.nc')) @@ -414,32 +373,11 @@ def test_lettered_tiles_no_fit(self): def test_lettered_tiles_no_valid_data(self): """Test creating a lettered grid with no valid data.""" from satpy.writers.awips_tiled import AWIPSTiledWriter - from xarray import DataArray - from pyresample.geometry import AreaDefinition - from pyresample.utils import proj4_str_to_dict w = AWIPSTiledWriter(base_dir=self.base_dir, compress=True) - area_def = AreaDefinition( - 'test', - 'test', - 'test', - proj4_str_to_dict('+proj=lcc +datum=WGS84 +ellps=WGS84 +lon_0=-95. ' - '+lat_0=25 +lat_1=25 +units=m +no_defs'), - 1000, - 2000, - (-1000000., -1500000., 1000000., 1500000.), - ) - now = datetime(2018, 1, 1, 12, 0, 0) - ds = DataArray( - da.full((2000, 1000), np.nan, chunks=500, dtype=np.float32), - attrs=dict( - name='test_ds', - platform_name='PLAT', - sensor='SENSOR', - units='1', - area=area_def, - start_time=now, - end_time=now + timedelta(minutes=20)) - ) + data = da.full((2000, 1000), np.nan, chunks=500, dtype=np.float32) + area_def = self._get_test_area(shape=(2000, 1000), + extents=(-1000000., -1500000., 1000000., 1500000.)) + ds = self._get_test_lcc_data(data, area_def) w.save_datasets([ds], sector_id='LCC', source_name="TESTS", tile_count=(3, 3), lettered_grid=True) # No files created - all NaNs should result in no tiles being created all_files = glob(os.path.join(self.base_dir, 'TESTS_AII*.nc')) @@ -448,31 +386,11 @@ def test_lettered_tiles_no_valid_data(self): def test_lettered_tiles_bad_filename(self): """Test creating a lettered grid with a bad filename.""" from satpy.writers.awips_tiled import AWIPSTiledWriter - from xarray import DataArray - from pyresample.geometry import AreaDefinition w = AWIPSTiledWriter(base_dir=self.base_dir, compress=True, filename="{Bad Key}.nc") - area_def = AreaDefinition( - 'test', - 'test', - 'test', - ('+proj=lcc +datum=WGS84 +ellps=WGS84 +lon_0=-95. ' - '+lat_0=25 +lat_1=25 +units=m +no_defs'), - 1000, - 2000, - (-1000000., -1500000., 1000000., 1500000.), - ) - now = datetime(2018, 1, 1, 12, 0, 0) - ds = DataArray( - da.from_array(np.linspace(0., 1., 2000000, dtype=np.float32).reshape((2000, 1000)), chunks=500), - attrs=dict( - name='test_ds', - platform_name='PLAT', - sensor='SENSOR', - units='1', - area=area_def, - start_time=now, - end_time=now + timedelta(minutes=20)) - ) + data = self._get_test_data(shape=(2000, 1000), chunks=500) + area_def = self._get_test_area(shape=(2000, 1000), + extents=(-1000000., -1500000., 1000000., 1500000.)) + ds = self._get_test_lcc_data(data, area_def) with pytest.raises(KeyError): w.save_datasets([ds], sector_id='LCC', @@ -484,33 +402,13 @@ def test_basic_numbered_tiles_rgb(self): """Test creating a multiple numbered tiles with RGB.""" from satpy.writers.awips_tiled import AWIPSTiledWriter import xarray as xr - from xarray import DataArray - from pyresample.geometry import AreaDefinition w = AWIPSTiledWriter(base_dir=self.base_dir, compress=True) - area_def = AreaDefinition( - 'test', - 'test', - 'test', - ('+proj=lcc +datum=WGS84 +ellps=WGS84 +lon_0=-95. ' - '+lat_0=25 +lat_1=25 +units=m +no_defs'), - 100, - 200, - (-1000., -1500., 1000., 1500.), - ) - now = datetime(2018, 1, 1, 12, 0, 0) - ds = DataArray( - da.from_array(np.linspace(0., 1., 60000, dtype=np.float32).reshape((3, 200, 100)), chunks=50), - dims=('bands', 'y', 'x'), - coords={'bands': ['R', 'G', 'B']}, - attrs=dict( - name='test_ds', - platform_name='PLAT', - sensor='SENSOR', - units='1', - area=area_def, - start_time=now, - end_time=now + timedelta(minutes=20)) - ) + data = da.from_array(np.linspace(0., 1., 60000, dtype=np.float32).reshape((3, 200, 100)), chunks=50) + area_def = self._get_test_area() + ds = self._get_test_lcc_data(data, area_def) + ds = ds.rename(dict((old, new) for old, new in zip(ds.dims, ['bands', 'y', 'x']))) + ds.coords['bands'] = ['R', 'G', 'B'] + w.save_datasets([ds], sector_id='TEST', source_name="TESTS", tile_count=(3, 3)) chan_files = glob(os.path.join(self.base_dir, 'TESTS_AII*test_ds_R*.nc')) all_files = chan_files[:] @@ -522,8 +420,9 @@ def test_basic_numbered_tiles_rgb(self): assert len(chan_files) == 9 all_files.extend(chan_files) for fn in all_files: - ds = xr.open_dataset(fn, mask_and_scale=False) - check_required_common_attributes(ds) + unmasked_ds = xr.open_dataset(fn, mask_and_scale=False) + masked_ds = xr.open_dataset(fn, mask_and_scale=True) + check_required_properties(unmasked_ds, masked_ds) @pytest.mark.parametrize( "sector", @@ -542,34 +441,18 @@ def test_multivar_numbered_tiles_glm(self, sector, extra_kwargs): """Test creating a tiles with multiple variables.""" import xarray as xr from satpy.writers.awips_tiled import AWIPSTiledWriter - from xarray import DataArray - from pyresample.geometry import AreaDefinition - from pyresample.utils import proj4_str_to_dict import os os.environ['ORGANIZATION'] = '1' * 50 w = AWIPSTiledWriter(base_dir=self.base_dir, compress=True) - area_def = AreaDefinition( - 'test', - 'test', - 'test', - proj4_str_to_dict('+proj=lcc +datum=WGS84 +ellps=WGS84 +lon_0=-95. ' - '+lat_0=25 +lat_1=25 +units=m +no_defs'), - 100, - 200, - (-1000., -1500., 1000., 1500.), - ) - now = datetime(2018, 1, 1, 12, 0, 0) - end_time = now + timedelta(minutes=20) - ds1 = DataArray( - da.from_array(np.linspace(0., 1., 20000, dtype=np.float32).reshape((200, 100)), chunks=50), - attrs=dict( + data = self._get_test_data() + area_def = self._get_test_area() + ds1 = self._get_test_lcc_data(data, area_def) + ds1.attrs.update( + dict( name='total_energy', platform_name='GOES-17', sensor='SENSOR', units='1', - area=area_def, - start_time=now, - end_time=end_time, scan_mode='M3', scene_abbr=sector, platform_shortname="G17" @@ -598,12 +481,13 @@ def test_multivar_numbered_tiles_glm(self, sector, extra_kwargs): all_files = glob(os.path.join(self.base_dir, fn_glob)) assert len(all_files) == 9 for fn in all_files: - ds = xr.open_dataset(fn, mask_and_scale=False) - check_required_common_attributes(ds) + unmasked_ds = xr.open_dataset(fn, mask_and_scale=False) + masked_ds = xr.open_dataset(fn, mask_and_scale=True) + check_required_properties(unmasked_ds, masked_ds) if sector == 'C': - assert ds.attrs['time_coverage_end'] == end_time.strftime('%Y-%m-%dT%H:%M:%S.%fZ') + assert masked_ds.attrs['time_coverage_end'] == self.end_time.strftime('%Y-%m-%dT%H:%M:%S.%fZ') else: # 'F' - assert ds.attrs['time_coverage_end'] == end_time.strftime('%Y-%m-%dT%H:%M:%SZ') + assert masked_ds.attrs['time_coverage_end'] == self.end_time.strftime('%Y-%m-%dT%H:%M:%SZ') @staticmethod def _get_glm_glob_filename(extra_kwargs): diff --git a/satpy/writers/awips_tiled.py b/satpy/writers/awips_tiled.py index 29bed7a321..94a58c58ac 100644 --- a/satpy/writers/awips_tiled.py +++ b/satpy/writers/awips_tiled.py @@ -227,7 +227,7 @@ from trollsift.parser import StringFormatter, Parser import numpy as np -from pyproj import Proj +from pyproj import Proj, CRS, Transformer import dask import dask.array as da import xarray as xr @@ -303,7 +303,7 @@ def _get_tile_properties(self, tile_shape, tile_count): else: raise ValueError("Either 'tile_count' or 'tile_shape' must be provided") - # number of pixels per each tile + # number of pixels per each tile (rows, cols) self.tile_shape = tile_shape # number of tiles in each direction (rows, columns) self.tile_count = tile_count @@ -342,9 +342,7 @@ def _get_xy_arrays(self): new_extents, ) - x, y = imaginary_grid_def.get_proj_coords() - x = x[0].squeeze() # all rows should have the same coordinates - y = y[:, 0].squeeze() # all columns should have the same coordinates + x, y = imaginary_grid_def.get_proj_vectors() return x, y def _get_xy_scaling_parameters(self): @@ -352,7 +350,7 @@ def _get_xy_scaling_parameters(self): gd = self.area_definition bx = self.x.min() mx = gd.pixel_size_x - by = self.y.min() + by = self.y.max() my = -abs(gd.pixel_size_y) return mx, bx, my, by @@ -424,10 +422,20 @@ def __call__(self): class LetteredTileGenerator(NumberedTileGenerator): """Helper class to generate per-tile metadata for lettered tiles.""" - def __init__(self, area_definition, extents, + def __init__(self, area_definition, extents, sector_crs, cell_size=(2000000, 2000000), num_subtiles=None, use_sector_reference=False): - """Initialize tile information for later generation.""" + """Initialize tile information for later generation. + + Args: + area_definition (AreaDefinition): Area of the data being saved. + extents (tuple): Four element tuple of the configured lettered + area. + sector_crs (pyproj.CRS): CRS of the configured lettered sector + area. + cell_size (tuple): Two element tuple of resolution of each tile + in sector projection units (y, x). + """ # (row subtiles, col subtiles) self.num_subtiles = num_subtiles or (2, 2) self.cell_size = cell_size # (row tile height, col tile width) @@ -435,7 +443,8 @@ def __init__(self, area_definition, extents, self.ll_extents = extents[:2] # (x min, y min) self.ur_extents = extents[2:] # (x max, y max) self.use_sector_reference = use_sector_reference - super(LetteredTileGenerator, self).__init__(area_definition) + self._transformer = Transformer.from_crs(sector_crs, area_definition.crs) + super().__init__(area_definition) def _get_tile_properties(self, tile_shape, tile_count): """Calculate tile information for this particular sector/grid.""" @@ -447,8 +456,8 @@ def _get_tile_properties(self, tile_shape, tile_count): ad = self.area_definition x, y = ad.get_proj_vectors() - ll_xy = self.ll_extents - ur_xy = self.ur_extents + ll_xy = self._transformer.transform(*self.ll_extents) + ur_xy = self._transformer.transform(*self.ur_extents) cw = abs(ad.pixel_size_x) ch = abs(ad.pixel_size_y) st = self.num_subtiles @@ -864,7 +873,6 @@ def render(self, dataset_or_data_arrays, shared_attrs=None): for data_arr in data_arrays: new_var_name, new_data_arr = self._render_variable(data_arr) new_ds[new_var_name] = new_data_arr - new_coords = self._render_coordinates(new_ds) new_ds.coords.update(new_coords) # use first data array as "representative" for global attributes @@ -911,6 +919,11 @@ def _global_start_date_time(self, input_metadata): def _global_awips_id(self, input_metadata): return "AWIPS_" + input_metadata['name'] + def _global_physical_element(self, input_metadata): + var_config = self._var_tree.find_match(**input_metadata) + result = self._render_attrs(var_config["attributes"], input_metadata) + return result["physical_element"] + def _global_production_location(self, input_metadata): """Get default global production_location attribute.""" del input_metadata @@ -964,6 +977,7 @@ def _render_variable_encoding(self, var_config, input_data_arr): new_encoding['scale_factor'] = sf new_encoding['add_offset'] = ao new_encoding['_FillValue'] = fill + new_encoding['coordinates'] = ' '.join([ele for ele in input_data_arr.dims]) return new_encoding def _get_projection_attrs(self, area_def): @@ -1263,6 +1277,7 @@ def separate_init_kwargs(cls, kwargs): def _fill_sector_info(self): """Convert sector extents if needed.""" for sector_info in self.awips_sectors.values(): + sector_info['projection'] = CRS.from_user_input(sector_info['projection']) p = Proj(sector_info['projection']) if 'lower_left_xy' in sector_info: sector_info['lower_left_lonlat'] = p(*sector_info['lower_left_xy'], inverse=True) @@ -1299,6 +1314,7 @@ def _get_tile_generator(self, area_def, lettered_grid, sector_id, tile_gen = LetteredTileGenerator( area_def, sector_info['lower_left_xy'] + sector_info['upper_right_xy'], + sector_crs=sector_info['projection'], cell_size=sector_info['resolution'], num_subtiles=num_subtiles, use_sector_reference=use_sector_reference, @@ -1522,11 +1538,11 @@ def save_datasets(self, datasets, sector_id=None, grid's pixels. By default this is False meaning that the grid's tiles will be shifted to align with the data locations. If True, the data is shifted. At most the data will be shifted - by 0.5 pixels. See :mod:`satpy.writers.scmi` for more + by 0.5 pixels. See :mod:`satpy.writers.awips_tiled` for more information. template (str or dict): Name of the template configured in the writer YAML file. This can also be a dictionary with a full - template configuration. See the :mod:`satpy.writers.scmi` + template configuration. See the :mod:`satpy.writers.awips_tiled` documentation for more information on templates. Defaults to the 'polar' builtin template. check_categories (bool): Whether category and flag products should @@ -1691,7 +1707,6 @@ def _create_debug_array(sector_info, num_subtiles, font_path='Verdana.ttf'): img.save("test.png") - from pyresample.utils import proj4_str_to_dict new_extents = ( ll_extent[0], ur_extent[1] - 1001. * meters_ppy, @@ -1702,7 +1717,7 @@ def _create_debug_array(sector_info, num_subtiles, font_path='Verdana.ttf'): 'debug_grid', 'debug_grid', 'debug_grid', - proj4_str_to_dict(sector_info['projection']), + sector_info['projection'], 1000, 1000, new_extents