From 9a64b12b22998d913737e9466008b5eded806f59 Mon Sep 17 00:00:00 2001 From: Stephan Finkensieper Date: Thu, 29 Nov 2018 15:20:08 +0100 Subject: [PATCH 01/12] Add method to compute geostationary earth mask --- satpy/readers/utils.py | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/satpy/readers/utils.py b/satpy/readers/utils.py index 1ccfa9a277..6114d4e327 100644 --- a/satpy/readers/utils.py +++ b/satpy/readers/utils.py @@ -26,6 +26,7 @@ import logging from contextlib import closing +import dask.array as da import tempfile import bz2 import os @@ -79,6 +80,35 @@ def get_geostationary_angle_extent(geos_area): return xmax, ymax +def get_geostationary_mask(area, flip=False): + """Compute a mask of the earth's shape as seen by a geostationary satellite + + Args: + area (pyresample.geometry.AreaDefinition) : Corresponding area + definition + flip (bool) : Flip mask vertically to account for scanning direction + (north-south vs south-north) + + Returns: + Boolean mask, True inside the earth's shape, False outside. + """ + # Compute projection coordinates at the earth's limb + h = area.proj_dict['h'] + xmax, ymax = get_geostationary_angle_extent(area) + xmax *= h + ymax *= h + + # Compute projection coordinates at the centre of each pixel + x, y = area.get_proj_coords_dask() + if flip: + # TODO: It would probably be faster to flip the 1D vector in + # get_proj_vectors_dask() before meshgridding + y = da.flipud(y) + + # Compute mask of the earth's elliptical shape + return ((x / xmax) ** 2 + (y / ymax) ** 2) <= 1 + + def _lonlat_from_geos_angle(x, y, geos_area): """Get lons and lats from x, y in projection coordinates.""" h = (geos_area.proj_dict['h'] + geos_area.proj_dict['a']) / 1000 From f7e23a8056435fb20f5475f541c7d846140f082d Mon Sep 17 00:00:00 2001 From: Stephan Finkensieper Date: Thu, 29 Nov 2018 15:45:37 +0100 Subject: [PATCH 02/12] Add JMA HRIT readers for MTSAT-1R/2 - Add YAML definitions for MTSAT-1R/2 - Update YAML definition for Himawari-8 - Rename hrit_jma.yaml -> hrit_ahi.yaml - Remove radiance datasets because it is not supported - Add counts datasets - Modify existing JMA HRIT reader: - Add support for half disk scans - Extract platform name from the file (needs to be improved for Himawari-9) - Mask space pixels --- .../readers/{hrit_jma.yaml => hrit_ahi.yaml} | 102 +++++------ satpy/etc/readers/hrit_mtsat1.yaml | 122 +++++++++++++ satpy/etc/readers/hrit_mtsat2.yaml | 122 +++++++++++++ satpy/readers/hrit_jma.py | 161 ++++++++++++++++-- 4 files changed, 442 insertions(+), 65 deletions(-) rename satpy/etc/readers/{hrit_jma.yaml => hrit_ahi.yaml} (80%) create mode 100644 satpy/etc/readers/hrit_mtsat1.yaml create mode 100644 satpy/etc/readers/hrit_mtsat2.yaml diff --git a/satpy/etc/readers/hrit_jma.yaml b/satpy/etc/readers/hrit_ahi.yaml similarity index 80% rename from satpy/etc/readers/hrit_jma.yaml rename to satpy/etc/readers/hrit_ahi.yaml index c6f1d8ac00..2d1e79c611 100644 --- a/satpy/etc/readers/hrit_jma.yaml +++ b/satpy/etc/readers/hrit_ahi.yaml @@ -1,6 +1,10 @@ +# References: +# - http://www.data.jma.go.jp/mscweb/en/himawari89/space_segment/spsg_ahi.html +# - http://www.data.jma.go.jp/mscweb/en/himawari89/space_segment/sample_hrit.html + reader: description: JMA HRIT Reader - name: hrit_jma + name: hrit_ahi sensors: [ahi] default_channels: [] reader: !!python/name:satpy.readers.yaml_reader.FileYAMLReader @@ -112,9 +116,9 @@ datasets: reflectance: standard_name: toa_bidirectional_reflectance units: "%" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: 1 file_type: hrit_b01 B02: @@ -126,9 +130,9 @@ datasets: reflectance: standard_name: toa_bidirectional_reflectance units: "%" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: 1 file_type: hrit_b02 B03: @@ -140,9 +144,9 @@ datasets: reflectance: standard_name: toa_bidirectional_reflectance units: "%" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: 1 file_type: hrit_b03 B04: @@ -154,9 +158,9 @@ datasets: reflectance: standard_name: toa_bidirectional_reflectance units: "%" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: 1 file_type: hrit_b04 B05: @@ -168,9 +172,9 @@ datasets: reflectance: standard_name: toa_bidirectional_reflectance units: "%" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: 1 file_type: hrit_b05 B06: @@ -182,9 +186,9 @@ datasets: reflectance: standard_name: toa_bidirectional_reflectance units: "%" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: 1 file_type: hrit_b06 B07: @@ -196,9 +200,9 @@ datasets: brightness_temperature: standard_name: toa_brightness_temperature units: "K" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: 1 file_type: hrit_b07 B08: @@ -210,9 +214,9 @@ datasets: brightness_temperature: standard_name: toa_brightness_temperature units: "K" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: 1 file_type: hrit_b08 B09: @@ -224,9 +228,9 @@ datasets: brightness_temperature: standard_name: toa_brightness_temperature units: "K" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: 1 file_type: hrit_b09 B10: @@ -238,9 +242,9 @@ datasets: brightness_temperature: standard_name: toa_brightness_temperature units: "K" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: 1 file_type: hrit_b10 B11: @@ -252,9 +256,9 @@ datasets: brightness_temperature: standard_name: toa_brightness_temperature units: "K" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: 1 file_type: hrit_b11 B12: @@ -266,9 +270,9 @@ datasets: brightness_temperature: standard_name: toa_brightness_temperature units: "K" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: 1 file_type: hrit_b12 B13: @@ -280,9 +284,9 @@ datasets: brightness_temperature: standard_name: toa_brightness_temperature units: "K" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: 1 file_type: hrit_b13 B14: @@ -294,9 +298,9 @@ datasets: brightness_temperature: standard_name: toa_brightness_temperature units: "K" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: 1 file_type: hrit_b14 B15: @@ -308,9 +312,9 @@ datasets: brightness_temperature: standard_name: toa_brightness_temperature units: "K" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: 1 file_type: hrit_b15 B16: @@ -322,7 +326,7 @@ datasets: brightness_temperature: standard_name: toa_brightness_temperature units: "K" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: 1 file_type: hrit_b16 diff --git a/satpy/etc/readers/hrit_mtsat1.yaml b/satpy/etc/readers/hrit_mtsat1.yaml new file mode 100644 index 0000000000..3fc6edfe9b --- /dev/null +++ b/satpy/etc/readers/hrit_mtsat1.yaml @@ -0,0 +1,122 @@ +# References: +# - https://www.wmo-sat.info/oscar/instruments/view/236 +# - http://www.data.jma.go.jp/mscweb/notice/Himawari7_e.html +# +# Note that the there exist two versions of the dataset. A segmented (data +# split into multiple files) and a non-segmented version (all data in one +# file). + +reader: + description: Reader for MTSAT-1R JAMI data in JMA HRIT Format + name: hrit_mtsat1 + sensors: [jami] + default_channels: [] + reader: !!python/name:satpy.readers.yaml_reader.FileYAMLReader + +file_types: + hrit_vis: + file_reader: !!python/name:satpy.readers.hrit_jma.HRITJMAFileHandler + file_patterns: + - 'IMG_DK01VIS_{start_time:%Y%m%d%H%M}_{segment:03d}' + - 'IMG_DK01VIS_{start_time:%Y%m%d%H%M}' + - 'HRIT_MTSAT1_{start_time:%Y%m%d_%H%M}_DK01VIS' + + hrit_ir1: + file_reader: !!python/name:satpy.readers.hrit_jma.HRITJMAFileHandler + file_patterns: + - 'IMG_DK01IR1_{start_time:%Y%m%d%H%M}_{segment:03d}' + - 'IMG_DK01IR1_{start_time:%Y%m%d%H%M}' + - 'HRIT_MTSAT1_{start_time:%Y%m%d_%H%M}_DK{area:02d}IR1' + + hrit_ir2: + file_reader: !!python/name:satpy.readers.hrit_jma.HRITJMAFileHandler + file_patterns: + - 'IMG_DK01IR2_{start_time:%Y%m%d%H%M}_{segment:03d}' + - 'IMG_DK01IR2_{start_time:%Y%m%d%H%M}' + - 'HRIT_MTSAT1_{start_time:%Y%m%d_%H%M}_DK01IR2' + + hrit_ir3: + file_reader: !!python/name:satpy.readers.hrit_jma.HRITJMAFileHandler + file_patterns: + - 'IMG_DK01IR3_{start_time:%Y%m%d%H%M}_{segment:03d}' + - 'IMG_DK01IR3_{start_time:%Y%m%d%H%M}' + - 'HRIT_MTSAT1_{start_time:%Y%m%d_%H%M}_DK01IR3' + + + hrit_ir4: + file_reader: !!python/name:satpy.readers.hrit_jma.HRITJMAFileHandler + file_patterns: + - 'IMG_DK01IR4_{start_time:%Y%m%d%H%M}_{segment:03d}' + - 'IMG_DK01IR4_{start_time:%Y%m%d%H%M}' + - 'HRIT_MTSAT1_{start_time:%Y%m%d_%H%M}_DK01IR4' + +datasets: + VIS: + name: VIS + sensor: jami + wavelength: [0.55, 0.675, 0.90] + resolution: 1000 + calibration: + counts: + standard_name: counts + units: 1 + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + file_type: hrit_vis + + IR1: + name: IR1 + sensor: jami + wavelength: [10.3, 10.8, 11.3] + resolution: 4000 + calibration: + counts: + standard_name: counts + units: 1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + file_type: hrit_ir1 + + IR2: + name: IR2 + sensor: jami + wavelength: [11.5, 12.0, 12.5] + resolution: 4000 + calibration: + counts: + standard_name: counts + units: 1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + file_type: hrit_ir2 + + IR3: + name: IR3 + sensor: jami + wavelength: [6.5, 6.75, 7.0] + resolution: 4000 + calibration: + counts: + standard_name: counts + units: 1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + file_type: hrit_ir3 + + IR4: + name: IR4 + sensor: jami + wavelength: [3.5, 3.75, 4.0] + resolution: 4000 + calibration: + counts: + standard_name: counts + units: 1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + file_type: hrit_ir4 \ No newline at end of file diff --git a/satpy/etc/readers/hrit_mtsat2.yaml b/satpy/etc/readers/hrit_mtsat2.yaml new file mode 100644 index 0000000000..b6eaa9d8d8 --- /dev/null +++ b/satpy/etc/readers/hrit_mtsat2.yaml @@ -0,0 +1,122 @@ +# References: +# - https://www.wmo-sat.info/oscar/instruments/view/219 +# - http://www.data.jma.go.jp/mscweb/notice/Himawari7_e.html +# +# Note that the there exist two versions of the dataset. A segmented (data +# split into multiple files) and a non-segmented version (all data in one +# file). + +reader: + description: Reader for MTSAT-2 data in JMA HRIT Format + name: hrit_mtsat2 + sensors: [mtsat2_imager] + default_channels: [] + reader: !!python/name:satpy.readers.yaml_reader.FileYAMLReader + +file_types: + hrit_vis: + file_reader: !!python/name:satpy.readers.hrit_jma.HRITJMAFileHandler + file_patterns: + - 'IMG_DK01VIS_{start_time:%Y%m%d%H%M}_{segment:03d}' + - 'IMG_DK01VIS_{start_time:%Y%m%d%H%M}' + - 'HRIT_MTSAT2_{start_time:%Y%m%d_%H%M}_DK01VIS' + + hrit_ir1: + file_reader: !!python/name:satpy.readers.hrit_jma.HRITJMAFileHandler + file_patterns: + - 'IMG_DK01IR1_{start_time:%Y%m%d%H%M}_{segment:03d}' + - 'IMG_DK01IR1_{start_time:%Y%m%d%H%M}' + - 'HRIT_MTSAT2_{start_time:%Y%m%d_%H%M}_DK{area:02d}IR1' + + hrit_ir2: + file_reader: !!python/name:satpy.readers.hrit_jma.HRITJMAFileHandler + file_patterns: + - 'IMG_DK01IR2_{start_time:%Y%m%d%H%M}_{segment:03d}' + - 'IMG_DK01IR2_{start_time:%Y%m%d%H%M}' + - 'HRIT_MTSAT2_{start_time:%Y%m%d_%H%M}_DK01IR2' + + hrit_ir3: + file_reader: !!python/name:satpy.readers.hrit_jma.HRITJMAFileHandler + file_patterns: + - 'IMG_DK01IR3_{start_time:%Y%m%d%H%M}_{segment:03d}' + - 'IMG_DK01IR3_{start_time:%Y%m%d%H%M}' + - 'HRIT_MTSAT2_{start_time:%Y%m%d_%H%M}_DK01IR3' + + + hrit_ir4: + file_reader: !!python/name:satpy.readers.hrit_jma.HRITJMAFileHandler + file_patterns: + - 'IMG_DK01IR4_{start_time:%Y%m%d%H%M}_{segment:03d}' + - 'IMG_DK01IR4_{start_time:%Y%m%d%H%M}' + - 'HRIT_MTSAT2_{start_time:%Y%m%d_%H%M}_DK01IR4' + +datasets: + VIS: + name: VIS + sensor: mtsat2_imager + wavelength: [0.55, 0.675, 0.80] + resolution: 1000 + calibration: + counts: + standard_name: counts + units: 1 + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + file_type: hrit_vis + + IR1: + name: IR1 + sensor: mtsat2_imager + wavelength: [10.3, 10.8, 11.3] + resolution: 4000 + calibration: + counts: + standard_name: counts + units: 1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + file_type: hrit_ir1 + + IR2: + name: IR2 + sensor: mtsat2_imager + wavelength: [11.5, 12.0, 12.5] + resolution: 4000 + calibration: + counts: + standard_name: counts + units: 1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + file_type: hrit_ir2 + + IR3: + name: IR3 + sensor: mtsat2_imager + wavelength: [6.5, 6.75, 7.0] + resolution: 4000 + calibration: + counts: + standard_name: counts + units: 1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + file_type: hrit_ir3 + + IR4: + name: IR4 + sensor: mtsat2_imager + wavelength: [3.5, 3.75, 4.0] + resolution: 4000 + calibration: + counts: + standard_name: counts + units: 1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + file_type: hrit_ir4 \ No newline at end of file diff --git a/satpy/readers/hrit_jma.py b/satpy/readers/hrit_jma.py index b81740df10..4e74204f47 100644 --- a/satpy/readers/hrit_jma.py +++ b/satpy/readers/hrit_jma.py @@ -39,6 +39,7 @@ from satpy.readers.hrit_base import (HRITFileHandler, ancillary_text, annotation_header, base_hdr_map, image_data_function) +from .utils import get_geostationary_mask logger = logging.getLogger('hrit_jma') @@ -86,6 +87,31 @@ ('microseconds', '>u2'), ('nanoseconds', '>u2')]) +FULL_DISK = 1 +NORTH_HEMIS = 2 +SOUTH_HEMIS = 3 +UNKNOWN_AREA = -1 +AREA_NAMES = {FULL_DISK: 'Full Disk', + NORTH_HEMIS: 'Northern Hemisphere', + SOUTH_HEMIS: 'Southern Hemisphere', + UNKNOWN_AREA: 'Unknown Area'} + +MTSAT1R = 'MTSAT-1R' +MTSAT2 = 'MTSAT-2' +HIMAWARI8 = 'Himawari-8' +UNKNOWN_PLATFORM = 'Unknown Platform' +PLATFORMS = { + 'GEOS(140.00)': MTSAT1R, + 'GEOS(140.25)': MTSAT1R, + 'GEOS(140.70)': HIMAWARI8, + 'GEOS(145.00)': MTSAT2, +} +SENSORS = { + MTSAT1R: 'jami', + MTSAT2: 'mtsat2_imager', + HIMAWARI8: 'ahi' +} + class HRITJMAFileHandler(HRITFileHandler): """JMA HRIT format reader.""" @@ -123,13 +149,102 @@ def __init__(self, filename, filename_info, filetype_info): self.projection_name = self.mda['projection_name'].decode().strip() sublon = float(self.projection_name.split('(')[1][:-1]) self.mda['projection_parameters']['SSP_longitude'] = sublon + self.platform = self._get_platform() + self.is_segmented = self.mda['segment_sequence_number'] > 0 + self.area_id = filename_info.get('area', UNKNOWN_AREA) + self.area = self._get_area_def() + + def _get_platform(self): + """Get the platform name + + The platform is not specified explicitly in JMA HRIT files. For + segmented data it is not even specified in the filename. But it + can be derived indirectly from the projection name: + + GEOS(140.00): MTSAT-1R + GEOS(140.25): MTSAT-1R # TODO: Check if there is more... + GEOS(140.70): Himawari-8 + GEOS(145.00): MTSAT-2 + + See [MTSAT], section 3.1. Unfortunately Himawari-8 and 9 are not + distinguishable using that method at the moment. From [HIMAWARI]: + + "HRIT/LRIT files have the same file naming convention in the same + format in Himawari-8 and Himawari-9, so there is no particular + difference." + + TODO: Find another way to distinguish Himawari-8 and 9. + + References: + [MTSAT] http://www.data.jma.go.jp/mscweb/notice/Himawari7_e.html + [HIMAWARI] http://www.data.jma.go.jp/mscweb/en/himawari89/space_segment/sample_hrit.html + """ + try: + return PLATFORMS[self.projection_name] + except KeyError: + logger.error('Unable to determine platform: Unknown projection ' + 'name "{}"'.format(self.projection_name)) + return UNKNOWN_PLATFORM + + def _check_sensor_platform_consistency(self, sensor): + """Make sure sensor and platform are consistent + + Args: + sensor (str) : Sensor name from YAML dataset definition + + Raises: + ValueError if they don't match + """ + ref_sensor = SENSORS.get(self.platform, None) + if ref_sensor and not sensor == ref_sensor: + logger.error('Sensor-Platform mismatch: {} is not a payload ' + 'of {}. Did you choose the correct reader?' + .format(sensor, self.platform)) + + def _get_line_offset(self): + """Get line offset for the current segment + + Line offset must be set so that (line - loff) is negative in the + northern hemisphere and positive in the southern hemisphere (-> image + centre has projection coordinates (x,y)=(0,0)) + """ + # Get line offset from the file + nlines = int(self.mda['number_of_lines']) + loff = np.float32(self.mda['loff']) - def get_area_def(self, dsid): + # Adapt it to the current segment + if self.is_segmented: + # Segmented data + segment_number = self.mda['segment_sequence_number'] - 1 + loff -= segment_number * nlines + elif self.area_id in (NORTH_HEMIS, SOUTH_HEMIS): + # Non-segmented data (segment_sequence_number=0). Half disk images + # are cutoff at the northern/southern edge (compared to the full + # disk image), but still have half the lines so that they overlap + # the equator by scanlines. Example for IR channels: + # + # Northern hemisphere: + # - Ideal line offset 1375 + # - Offset in the file: ca 50 (due to cutoff at the northern + # edge) + # - Adapted line offset: ca 1325 -> last scanline a little bit + # below the equator + # Southern Hemisphere: + # - Ideal line offset 0 + # - Offset in the file: ca 1325 (= 1375 - ca. 50 due due to + # cutoff at the southern edge) + # - Adapted line offset: ca 50 -> first scanline a little bit + # above the equator + loff = nlines - loff + + return loff + + def _get_area_def(self): """Get the area definition of the band.""" cfac = np.int32(self.mda['cfac']) lfac = np.int32(self.mda['lfac']) coff = np.float32(self.mda['coff']) - loff = np.float32(self.mda['loff']) + loff = self._get_line_offset() a = self.mda['projection_parameters']['a'] b = self.mda['projection_parameters']['b'] @@ -139,10 +254,6 @@ def get_area_def(self, dsid): nlines = int(self.mda['number_of_lines']) ncols = int(self.mda['number_of_columns']) - segment_number = self.mda['segment_sequence_number'] - 1 - - loff -= segment_number * nlines - area_extent = self.get_area_extent((nlines, ncols), (loff, coff), (lfac, cfac), @@ -155,29 +266,47 @@ def get_area_def(self, dsid): 'proj': 'geos', 'units': 'm'} + area_name = '{} {}'.format(self.platform, AREA_NAMES[self.area_id]) area = geometry.AreaDefinition( - "FLDK", - "HRIT FLDK Area: {}".format(self.projection_name), - 'geosmsg', - proj_dict, - ncols, - nlines, - area_extent) + area_id=area_name.replace(' ', '-').lower(), + name=area_name, + proj_id='geosmsg', + proj_dict=proj_dict, + x_size=ncols, + y_size=nlines, + area_extent=area_extent) + return area + def get_area_def(self, dsid): + """Get the area definition of the band.""" + return self.area + def get_dataset(self, key, info): """Get the dataset designated by *key*.""" res = super(HRITJMAFileHandler, self).get_dataset(key, info) - res = self.calibrate(res, key.calibration) + # Filenames of segmented data is identical for MTSAT-1R, MTSAT-2 + # and Himawari-8/9. Make sure we have the correct reader for the data + # at hand. + self._check_sensor_platform_consistency(info['sensor']) + + res = self._mask_space(self.calibrate(res, key.calibration)) + res.attrs.update(info) - res.attrs['platform_name'] = 'Himawari-8' - res.attrs['sensor'] = 'ahi' + res.attrs['platform_name'] = self.platform res.attrs['satellite_longitude'] = float(self.mda['projection_parameters']['SSP_longitude']) res.attrs['satellite_latitude'] = 0. res.attrs['satellite_altitude'] = float(self.mda['projection_parameters']['h']) + return res + def _mask_space(self, data): + """Mask space pixels""" + geomask = get_geostationary_mask(area=self.area, + flip=self.is_segmented) + return data.where(geomask) + def calibrate(self, data, calibration): """Calibrate the data.""" tic = datetime.now() From 7eff02f6ebf44b85fa89cad9df88e42aa862974e Mon Sep 17 00:00:00 2001 From: Stephan Finkensieper Date: Thu, 29 Nov 2018 17:17:52 +0100 Subject: [PATCH 03/12] Fix area names in YAML definitions --- satpy/etc/readers/hrit_mtsat1.yaml | 28 ++++++++++++++-------------- satpy/etc/readers/hrit_mtsat2.yaml | 28 ++++++++++++++-------------- 2 files changed, 28 insertions(+), 28 deletions(-) diff --git a/satpy/etc/readers/hrit_mtsat1.yaml b/satpy/etc/readers/hrit_mtsat1.yaml index 3fc6edfe9b..f3f81bcf5b 100644 --- a/satpy/etc/readers/hrit_mtsat1.yaml +++ b/satpy/etc/readers/hrit_mtsat1.yaml @@ -17,38 +17,38 @@ file_types: hrit_vis: file_reader: !!python/name:satpy.readers.hrit_jma.HRITJMAFileHandler file_patterns: - - 'IMG_DK01VIS_{start_time:%Y%m%d%H%M}_{segment:03d}' - - 'IMG_DK01VIS_{start_time:%Y%m%d%H%M}' - - 'HRIT_MTSAT1_{start_time:%Y%m%d_%H%M}_DK01VIS' + - 'IMG_DK{area:02d}VIS_{start_time:%Y%m%d%H%M}_{segment:03d}' + - 'IMG_DK{area:02d}VIS_{start_time:%Y%m%d%H%M}' + - 'HRIT_MTSAT1_{start_time:%Y%m%d_%H%M}_DK{area:02d}VIS' hrit_ir1: file_reader: !!python/name:satpy.readers.hrit_jma.HRITJMAFileHandler file_patterns: - - 'IMG_DK01IR1_{start_time:%Y%m%d%H%M}_{segment:03d}' - - 'IMG_DK01IR1_{start_time:%Y%m%d%H%M}' + - 'IMG_DK{area:02d}IR1_{start_time:%Y%m%d%H%M}_{segment:03d}' + - 'IMG_DK{area:02d}IR1_{start_time:%Y%m%d%H%M}' - 'HRIT_MTSAT1_{start_time:%Y%m%d_%H%M}_DK{area:02d}IR1' hrit_ir2: file_reader: !!python/name:satpy.readers.hrit_jma.HRITJMAFileHandler file_patterns: - - 'IMG_DK01IR2_{start_time:%Y%m%d%H%M}_{segment:03d}' - - 'IMG_DK01IR2_{start_time:%Y%m%d%H%M}' - - 'HRIT_MTSAT1_{start_time:%Y%m%d_%H%M}_DK01IR2' + - 'IMG_DK{area:02d}IR2_{start_time:%Y%m%d%H%M}_{segment:03d}' + - 'IMG_DK{area:02d}IR2_{start_time:%Y%m%d%H%M}' + - 'HRIT_MTSAT1_{start_time:%Y%m%d_%H%M}_DK{area:02d}IR2' hrit_ir3: file_reader: !!python/name:satpy.readers.hrit_jma.HRITJMAFileHandler file_patterns: - - 'IMG_DK01IR3_{start_time:%Y%m%d%H%M}_{segment:03d}' - - 'IMG_DK01IR3_{start_time:%Y%m%d%H%M}' - - 'HRIT_MTSAT1_{start_time:%Y%m%d_%H%M}_DK01IR3' + - 'IMG_DK{area:02d}IR3_{start_time:%Y%m%d%H%M}_{segment:03d}' + - 'IMG_DK{area:02d}IR3_{start_time:%Y%m%d%H%M}' + - 'HRIT_MTSAT1_{start_time:%Y%m%d_%H%M}_DK{area:02d}IR3' hrit_ir4: file_reader: !!python/name:satpy.readers.hrit_jma.HRITJMAFileHandler file_patterns: - - 'IMG_DK01IR4_{start_time:%Y%m%d%H%M}_{segment:03d}' - - 'IMG_DK01IR4_{start_time:%Y%m%d%H%M}' - - 'HRIT_MTSAT1_{start_time:%Y%m%d_%H%M}_DK01IR4' + - 'IMG_DK{area:02d}IR4_{start_time:%Y%m%d%H%M}_{segment:03d}' + - 'IMG_DK{area:02d}IR4_{start_time:%Y%m%d%H%M}' + - 'HRIT_MTSAT1_{start_time:%Y%m%d_%H%M}_DK{area:02d}IR4' datasets: VIS: diff --git a/satpy/etc/readers/hrit_mtsat2.yaml b/satpy/etc/readers/hrit_mtsat2.yaml index b6eaa9d8d8..015bc81341 100644 --- a/satpy/etc/readers/hrit_mtsat2.yaml +++ b/satpy/etc/readers/hrit_mtsat2.yaml @@ -17,38 +17,38 @@ file_types: hrit_vis: file_reader: !!python/name:satpy.readers.hrit_jma.HRITJMAFileHandler file_patterns: - - 'IMG_DK01VIS_{start_time:%Y%m%d%H%M}_{segment:03d}' - - 'IMG_DK01VIS_{start_time:%Y%m%d%H%M}' - - 'HRIT_MTSAT2_{start_time:%Y%m%d_%H%M}_DK01VIS' + - 'IMG_DK{area:02d}VIS_{start_time:%Y%m%d%H%M}_{segment:03d}' + - 'IMG_DK{area:02d}VIS_{start_time:%Y%m%d%H%M}' + - 'HRIT_MTSAT2_{start_time:%Y%m%d_%H%M}_DK{area:02d}VIS' hrit_ir1: file_reader: !!python/name:satpy.readers.hrit_jma.HRITJMAFileHandler file_patterns: - - 'IMG_DK01IR1_{start_time:%Y%m%d%H%M}_{segment:03d}' - - 'IMG_DK01IR1_{start_time:%Y%m%d%H%M}' + - 'IMG_DK{area:02d}IR1_{start_time:%Y%m%d%H%M}_{segment:03d}' + - 'IMG_DK{area:02d}IR1_{start_time:%Y%m%d%H%M}' - 'HRIT_MTSAT2_{start_time:%Y%m%d_%H%M}_DK{area:02d}IR1' hrit_ir2: file_reader: !!python/name:satpy.readers.hrit_jma.HRITJMAFileHandler file_patterns: - - 'IMG_DK01IR2_{start_time:%Y%m%d%H%M}_{segment:03d}' - - 'IMG_DK01IR2_{start_time:%Y%m%d%H%M}' - - 'HRIT_MTSAT2_{start_time:%Y%m%d_%H%M}_DK01IR2' + - 'IMG_DK{area:02d}IR2_{start_time:%Y%m%d%H%M}_{segment:03d}' + - 'IMG_DK{area:02d}IR2_{start_time:%Y%m%d%H%M}' + - 'HRIT_MTSAT2_{start_time:%Y%m%d_%H%M}_DK{area:02d}IR2' hrit_ir3: file_reader: !!python/name:satpy.readers.hrit_jma.HRITJMAFileHandler file_patterns: - - 'IMG_DK01IR3_{start_time:%Y%m%d%H%M}_{segment:03d}' - - 'IMG_DK01IR3_{start_time:%Y%m%d%H%M}' - - 'HRIT_MTSAT2_{start_time:%Y%m%d_%H%M}_DK01IR3' + - 'IMG_DK{area:02d}IR3_{start_time:%Y%m%d%H%M}_{segment:03d}' + - 'IMG_DK{area:02d}IR3_{start_time:%Y%m%d%H%M}' + - 'HRIT_MTSAT2_{start_time:%Y%m%d_%H%M}_DK{area:02d}IR3' hrit_ir4: file_reader: !!python/name:satpy.readers.hrit_jma.HRITJMAFileHandler file_patterns: - - 'IMG_DK01IR4_{start_time:%Y%m%d%H%M}_{segment:03d}' - - 'IMG_DK01IR4_{start_time:%Y%m%d%H%M}' - - 'HRIT_MTSAT2_{start_time:%Y%m%d_%H%M}_DK01IR4' + - 'IMG_DK{area:02d}IR4_{start_time:%Y%m%d%H%M}_{segment:03d}' + - 'IMG_DK{area:02d}IR4_{start_time:%Y%m%d%H%M}' + - 'HRIT_MTSAT2_{start_time:%Y%m%d_%H%M}_DK{area:02d}IR4' datasets: VIS: From e57a4640208ff061071ee58874fdbf2939778113 Mon Sep 17 00:00:00 2001 From: Stephan Finkensieper Date: Thu, 29 Nov 2018 17:18:22 +0100 Subject: [PATCH 04/12] Make area IDs short --- satpy/readers/hrit_jma.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/satpy/readers/hrit_jma.py b/satpy/readers/hrit_jma.py index 4e74204f47..2a3ae62137 100644 --- a/satpy/readers/hrit_jma.py +++ b/satpy/readers/hrit_jma.py @@ -91,10 +91,10 @@ NORTH_HEMIS = 2 SOUTH_HEMIS = 3 UNKNOWN_AREA = -1 -AREA_NAMES = {FULL_DISK: 'Full Disk', - NORTH_HEMIS: 'Northern Hemisphere', - SOUTH_HEMIS: 'Southern Hemisphere', - UNKNOWN_AREA: 'Unknown Area'} +AREA_NAMES = {FULL_DISK: {'short': 'FLDK', 'long': 'Full Disk'}, + NORTH_HEMIS: {'short': 'NH', 'long': 'Northern Hemisphere'}, + SOUTH_HEMIS: {'short': 'SH', 'long': 'Southern Hemisphere'}, + UNKNOWN_AREA: {'short': 'UNKNOWN', 'long': 'Unknown Area'}} MTSAT1R = 'MTSAT-1R' MTSAT2 = 'MTSAT-2' @@ -266,10 +266,9 @@ def _get_area_def(self): 'proj': 'geos', 'units': 'm'} - area_name = '{} {}'.format(self.platform, AREA_NAMES[self.area_id]) area = geometry.AreaDefinition( - area_id=area_name.replace(' ', '-').lower(), - name=area_name, + area_id=AREA_NAMES[self.area_id]['short'], + name=AREA_NAMES[self.area_id]['long'], proj_id='geosmsg', proj_dict=proj_dict, x_size=ncols, From 8a53cf8bde53cf373842c5f78047f7bbd0a7c8f3 Mon Sep 17 00:00:00 2001 From: Stephan Finkensieper Date: Fri, 30 Nov 2018 15:56:01 +0100 Subject: [PATCH 05/12] Remove flip workaround --- satpy/readers/utils.py | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/satpy/readers/utils.py b/satpy/readers/utils.py index 6114d4e327..ebbcd6cce5 100644 --- a/satpy/readers/utils.py +++ b/satpy/readers/utils.py @@ -80,14 +80,12 @@ def get_geostationary_angle_extent(geos_area): return xmax, ymax -def get_geostationary_mask(area, flip=False): +def get_geostationary_mask(area): """Compute a mask of the earth's shape as seen by a geostationary satellite Args: area (pyresample.geometry.AreaDefinition) : Corresponding area definition - flip (bool) : Flip mask vertically to account for scanning direction - (north-south vs south-north) Returns: Boolean mask, True inside the earth's shape, False outside. @@ -100,10 +98,6 @@ def get_geostationary_mask(area, flip=False): # Compute projection coordinates at the centre of each pixel x, y = area.get_proj_coords_dask() - if flip: - # TODO: It would probably be faster to flip the 1D vector in - # get_proj_vectors_dask() before meshgridding - y = da.flipud(y) # Compute mask of the earth's elliptical shape return ((x / xmax) ** 2 + (y / ymax) ** 2) <= 1 From fca455e0fe934729556fd7d938edeb5f456e5b46 Mon Sep 17 00:00:00 2001 From: Stephan Finkensieper Date: Fri, 30 Nov 2018 18:23:56 +0100 Subject: [PATCH 06/12] Catch unknown area ID --- satpy/readers/hrit_jma.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/satpy/readers/hrit_jma.py b/satpy/readers/hrit_jma.py index 2a3ae62137..5b00ddf3f5 100644 --- a/satpy/readers/hrit_jma.py +++ b/satpy/readers/hrit_jma.py @@ -152,6 +152,8 @@ def __init__(self, filename, filename_info, filetype_info): self.platform = self._get_platform() self.is_segmented = self.mda['segment_sequence_number'] > 0 self.area_id = filename_info.get('area', UNKNOWN_AREA) + if self.area_id not in AREA_NAMES: + self.area_id = UNKNOWN_AREA self.area = self._get_area_def() def _get_platform(self): From 4ff994779308a934965c28da24f67e4f7d06b177 Mon Sep 17 00:00:00 2001 From: Stephan Finkensieper Date: Fri, 30 Nov 2018 18:24:48 +0100 Subject: [PATCH 07/12] Fix area definitions and space masking --- satpy/readers/hrit_jma.py | 40 ++++++++++++++++----------------------- 1 file changed, 16 insertions(+), 24 deletions(-) diff --git a/satpy/readers/hrit_jma.py b/satpy/readers/hrit_jma.py index 5b00ddf3f5..a35025c063 100644 --- a/satpy/readers/hrit_jma.py +++ b/satpy/readers/hrit_jma.py @@ -206,9 +206,12 @@ def _check_sensor_platform_consistency(self, sensor): def _get_line_offset(self): """Get line offset for the current segment - Line offset must be set so that (line - loff) is negative in the - northern hemisphere and positive in the southern hemisphere (-> image - centre has projection coordinates (x,y)=(0,0)) + Read line offset from the file and adapt it to the current segment + or half disk scan so that + + y(l) ~ l - loff + + because this is what get_geostationary_area_extent() expects. """ # Get line offset from the file nlines = int(self.mda['number_of_lines']) @@ -216,28 +219,16 @@ def _get_line_offset(self): # Adapt it to the current segment if self.is_segmented: - # Segmented data + # loff in the file specifies the offset of the full disk image + # centre (1375/2750 for VIS/IR) segment_number = self.mda['segment_sequence_number'] - 1 - loff -= segment_number * nlines + loff -= (self.mda['total_no_image_segm'] - segment_number - 1) * nlines elif self.area_id in (NORTH_HEMIS, SOUTH_HEMIS): - # Non-segmented data (segment_sequence_number=0). Half disk images - # are cutoff at the northern/southern edge (compared to the full - # disk image), but still have half the lines so that they overlap - # the equator by scanlines. Example for IR channels: - # - # Northern hemisphere: - # - Ideal line offset 1375 - # - Offset in the file: ca 50 (due to cutoff at the northern - # edge) - # - Adapted line offset: ca 1325 -> last scanline a little bit - # below the equator - # Southern Hemisphere: - # - Ideal line offset 0 - # - Offset in the file: ca 1325 (= 1375 - ca. 50 due due to - # cutoff at the southern edge) - # - Adapted line offset: ca 50 -> first scanline a little bit - # above the equator + # loff in the file specifies the start line of the half disk image + # in the full disk image loff = nlines - loff + elif self.area_id == UNKNOWN_AREA: + logger.error('Cannot compute line offset for unknown area') return loff @@ -292,8 +283,10 @@ def get_dataset(self, key, info): # at hand. self._check_sensor_platform_consistency(info['sensor']) + # Calibrate and mask space pixels res = self._mask_space(self.calibrate(res, key.calibration)) + # Update attributes res.attrs.update(info) res.attrs['platform_name'] = self.platform res.attrs['satellite_longitude'] = float(self.mda['projection_parameters']['SSP_longitude']) @@ -304,8 +297,7 @@ def get_dataset(self, key, info): def _mask_space(self, data): """Mask space pixels""" - geomask = get_geostationary_mask(area=self.area, - flip=self.is_segmented) + geomask = get_geostationary_mask(area=self.area) return data.where(geomask) def calibrate(self, data, calibration): From 2ead428b48665e14503f6242ac9ce3ade8dc8ef8 Mon Sep 17 00:00:00 2001 From: Stephan Finkensieper Date: Fri, 30 Nov 2018 18:26:19 +0100 Subject: [PATCH 08/12] Add/update tests for hrit_jma --- satpy/tests/reader_tests/test_hrit_jma.py | 285 +++++++++++++++++----- 1 file changed, 224 insertions(+), 61 deletions(-) diff --git a/satpy/tests/reader_tests/test_hrit_jma.py b/satpy/tests/reader_tests/test_hrit_jma.py index b8f69be7c2..10a39ec8af 100644 --- a/satpy/tests/reader_tests/test_hrit_jma.py +++ b/satpy/tests/reader_tests/test_hrit_jma.py @@ -37,81 +37,244 @@ class TestHRITJMAFileHandler(unittest.TestCase): """Test the HRITJMAFileHandler.""" @mock.patch('satpy.readers.hrit_jma.HRITFileHandler.__init__') - def setUp(self, new_fh_init): - """Setup the hrit file handler for testing.""" + def _get_reader(self, mocked_init, mda, filename_info=None): from satpy.readers.hrit_jma import HRITJMAFileHandler - mda = { - 'projection_parameters': { - 'a': 6378169.00, - 'b': 6356583.80, - 'h': 35785831.00, - }, - 'image_data_function': b'$HALFTONE:=16\r_NAME:=VISIBLE\r_UNIT:=ALBEDO(%)\r' - b'0:=-0.10\r1023:=100.00\r65535:=100.00\r', - 'image_segm_seq_no': 0, - 'total_no_image_segm': 1, - 'projection_name': b'GEOS(140.70) ', - 'cfac': 40932549, - 'lfac': 40932549, - 'coff': 5500, - 'loff': 5500, - 'number_of_columns': 11000, - 'number_of_lines': 11000, - } + if not filename_info: + filename_info = {} HRITJMAFileHandler.filename = 'filename' HRITJMAFileHandler.mda = mda - self.reader = HRITJMAFileHandler('filename', {}, {}) + return HRITJMAFileHandler('filename', filename_info, {}) + + def _get_mda(self, loff=5500.0, coff=5500.0, nlines=11000, ncols=11000, + segno=0, numseg=1, vis=True): + """Create metadata dict like HRITFileHandler would do it""" + if vis: + idf = b'$HALFTONE:=16\r_NAME:=VISIBLE\r_UNIT:=ALBEDO(%)\r' \ + b'0:=-0.10\r1023:=100.00\r65535:=100.00\r' + else: + idf = b'$HALFTONE:=16\r_NAME:=INFRARED\r_UNIT:=KELVIN\r' \ + b'0:=329.98\r1023:=130.02\r65535:=130.02\r' + + return {'image_segm_seq_no': segno, + 'total_no_image_segm': numseg, + 'projection_name': b'GEOS(140.70) ', + 'projection_parameters': { + 'a': 6378169.00, + 'b': 6356583.80, + 'h': 35785831.00, + }, + 'cfac': 10233128, + 'lfac': 10233128, + 'coff': coff, + 'loff': loff, + 'number_of_columns': ncols, + 'number_of_lines': nlines, + 'image_data_function': idf + } def test_init(self): """Test creating the file handler.""" - mda = { - 'image_segm_seq_no': 0, - 'planned_end_segment_number': 1, - 'planned_start_segment_number': 1, - 'segment_sequence_number': 0, - 'total_no_image_segm': 1, - 'unit': 'ALBEDO(%)', - 'projection_name': b'GEOS(140.70) ', - 'projection_parameters': { - 'a': 6378169.00, - 'b': 6356583.80, - 'h': 35785831.00, - 'SSP_longitude': 140.7 - }, - 'cfac': 40932549, - 'lfac': 40932549, - 'coff': 5500, - 'loff': 5500, - 'number_of_columns': 11000, - 'number_of_lines': 11000, - 'image_data_function': b'$HALFTONE:=16\r_NAME:=VISIBLE\r_UNIT:=ALBEDO(%)\r' - b'0:=-0.10\r1023:=100.00\r65535:=100.00\r', - } - self.assertEqual(self.reader.mda, mda) + from satpy.readers.hrit_jma import UNKNOWN_AREA, HIMAWARI8 + + # Test addition of extra metadata + mda = self._get_mda() + mda_expected = mda.copy() + mda_expected.update( + {'planned_end_segment_number': 1, + 'planned_start_segment_number': 1, + 'segment_sequence_number': 0, + 'unit': 'ALBEDO(%)'}) + mda_expected['projection_parameters']['SSP_longitude'] = 140.7 + reader = self._get_reader(mda=mda) + self.assertEqual(reader.mda, mda_expected) + + # Check projection name + self.assertEqual(reader.projection_name, 'GEOS(140.70)') + + # Check calibration table + cal_expected = np.array([[0, -0.1], + [1023, 100], + [65535, 100]]) + self.assertTrue(np.all(reader.calibration_table == cal_expected)) + + # Check platform + self.assertEqual(reader.platform, HIMAWARI8) + + # Check is_segmented attribute + expected = {0: False, 1: True, 8: True} + for segno, is_segmented in expected.items(): + mda = self._get_mda(segno=segno) + reader = self._get_reader(mda=mda) + self.assertEqual(reader.is_segmented, is_segmented) + + # Check area IDs + expected = [ + ({'area': 1}, 1), + ({'area': 1234}, UNKNOWN_AREA), + ({}, UNKNOWN_AREA) + ] + mda = self._get_mda() + for filename_info, area_id in expected: + reader = self._get_reader(mda=mda, filename_info=filename_info) + self.assertEqual(reader.area_id, area_id) + + @mock.patch('satpy.readers.hrit_jma.HRITJMAFileHandler.__init__') + def test_get_platform(self, mocked_init): + """Test platform identification""" + from satpy.readers.hrit_jma import HRITJMAFileHandler + from satpy.readers.hrit_jma import PLATFORMS, UNKNOWN_PLATFORM + + mocked_init.return_value = None + reader = HRITJMAFileHandler() + + for proj_name, platform in PLATFORMS.items(): + reader.projection_name = proj_name + self.assertEqual(reader._get_platform(), platform) + + with mock.patch('logging.Logger.error') as mocked_log: + reader.projection_name = 'invalid' + self.assertEqual(reader._get_platform(), UNKNOWN_PLATFORM) + mocked_log.assert_called() + + def test_get_area_def(self): + """Test getting an AreaDefinition.""" + from satpy.readers.hrit_jma import FULL_DISK, NORTH_HEMIS, SOUTH_HEMIS + + cases = [ + # Non-segmented, full disk + {'loff': 1375.0, 'coff': 1375.0, + 'nlines': 2750, 'ncols': 2750, + 'segno': 0, 'numseg': 1, + 'area': FULL_DISK, + 'extent': (-5498000.088960204, -5498000.088960204, + 5502000.089024927, 5502000.089024927)}, + # Non-segmented, northern hemisphere + {'loff': 1325.0, 'coff': 1375.0, + 'nlines': 1375, 'ncols': 2750, + 'segno': 0, 'numseg': 1, + 'area': NORTH_HEMIS, + 'extent': (-5498000.088960204, -198000.00320373234, + 5502000.089024927, 5302000.085788833)}, + # Non-segmented, southern hemisphere + {'loff': 50, 'coff': 1375.0, + 'nlines': 1375, 'ncols': 2750, + 'segno': 0, 'numseg': 1, + 'area': SOUTH_HEMIS, + 'extent': (-5498000.088960204, -5298000.085724112, + 5502000.089024927, 202000.0032684542)}, + # Segmented, segment #1 + {'loff': 1375.0, 'coff': 1375.0, + 'nlines': 275, 'ncols': 2750, + 'segno': 1, 'numseg': 10, + 'area': FULL_DISK, + 'extent': (-5498000.088960204, 4402000.071226413, + 5502000.089024927, 5502000.089024927)}, + # Segmented, segment #7 + {'loff': 1375.0, 'coff': 1375.0, + 'nlines': 275, 'ncols': 2750, + 'segno': 7, 'numseg': 10, + 'area': FULL_DISK, + 'extent': (-5498000.088960204, -2198000.035564665, + 5502000.089024927, -1098000.0177661523)}, + ] + for case in cases: + mda = self._get_mda(loff=case['loff'], coff=case['coff'], + nlines=case['nlines'], ncols=case['ncols'], + segno=case['segno'], numseg=case['numseg']) + reader = self._get_reader(mda=mda, + filename_info={'area': case['area']}) + + self.assertTupleEqual(reader._get_area_def().area_extent, + case['extent']) + + def test_calibrate(self): + """Test calibration""" + # Generate test data + counts = DataArray(da.linspace(0, 1200, 25, chunks=5).reshape(5, 5)) + refl = np.array( + [[np.nan, 4.79247312, 9.68494624, 14.57741935, 19.46989247], + [24.36236559, 29.25483871, 34.14731183, 39.03978495, 43.93225806], + [48.82473118, 53.7172043, 58.60967742, 63.50215054, 68.39462366], + [73.28709677, 78.17956989, 83.07204301, 87.96451613, 92.85698925], + [97.74946237, 100., 100., 100., 100.]] + ) + bt = np.array( + [[np.nan, 320.20678397, 310.43356794, 300.66035191, 290.88713587], + [281.11391984, 271.34070381, 261.56748778, 251.79427175, 242.02105572], + [232.24783969, 222.47462366, 212.70140762, 202.92819159, 193.15497556], + [183.38175953, 173.6085435, 163.83532747, 154.06211144, 144.28889541], + [134.51567937, 130.02, 130.02, 130.02, 130.02]] + ) + + # Choose an area near the subsatellite point to avoid masking + # of space pixels + mda = self._get_mda(nlines=5, ncols=5, loff=1375.0, coff=1375.0, + segno=0) + reader = self._get_reader(mda=mda) + + # 1. Counts + res = reader.calibrate(data=counts, calibration='counts') + self.assertTrue(np.all(counts.values == res.values)) + + # 2. Reflectance + res = reader.calibrate(data=counts, calibration='reflectance') + np.testing.assert_allclose(refl, res.values) # also compares NaN + + # 3. Brightness temperature + mda_bt = self._get_mda(nlines=5, ncols=5, loff=1375.0, coff=1375.0, + segno=0, vis=False) + reader_bt = self._get_reader(mda=mda_bt) + res = reader_bt.calibrate(data=counts, + calibration='brightness_temperature') + np.testing.assert_allclose(bt, res.values) # also compares NaN + + def test_mask_space(self): + """Test masking of space pixels""" + mda = self._get_mda(loff=1375.0, coff=1375.0, nlines=275, ncols=1375, + segno=1, numseg=10) + reader = self._get_reader(mda=mda) + data = DataArray(da.ones((275, 1375), chunks=1024)) + masked = reader._mask_space(data) + + # First line of the segment should be space, in the middle of the + # last line there should be some valid pixels + np.testing.assert_allclose(masked.values[0, :], np.nan) + self.assertTrue(np.all(masked.values[-1, 588:788] == 1)) @mock.patch('satpy.readers.hrit_jma.HRITFileHandler.get_dataset') def test_get_dataset(self, base_get_dataset): - """Test getting a reflectance DataArray.""" + """Test getting a dataset""" + from satpy.readers.hrit_jma import HIMAWARI8 + + mda = self._get_mda(loff=1375.0, coff=1375.0, nlines=275, ncols=1375, + segno=1, numseg=10) + reader = self._get_reader(mda=mda) + key = mock.MagicMock() key.calibration = 'reflectance' - base_get_dataset.return_value = DataArray(da.arange(25, chunks=5).reshape(5, 5)) - res = self.reader.get_dataset(key, {'units': '%'}) - expected = np.array([ - [np.nan, -2.15053763e-03, 9.56989247e-02, 1.93548387e-01, 2.91397849e-01], - [3.89247312e-01, 4.87096774e-01, 5.84946237e-01, 6.82795699e-01, 7.80645161e-01], - [8.78494624e-01, 9.76344086e-01, 1.07419355e+00, 1.17204301e+00, 1.26989247e+00], - [1.36774194e+00, 1.46559140e+00, 1.56344086e+00, 1.66129032e+00, 1.75913978e+00], - [1.85698925e+00, 1.95483871e+00, 2.05268817e+00, 2.15053763e+00, 2.24838710e+00]]) - np.testing.assert_allclose(res.values, expected) + + base_get_dataset.return_value = DataArray(da.ones((275, 1375), + chunks=1024)) + + # Check attributes + res = reader.get_dataset(key, {'units': '%', 'sensor': 'ahi'}) self.assertEqual(res.attrs['units'], '%') + self.assertEqual(res.attrs['sensor'], 'ahi') + self.assertEqual(res.attrs['platform_name'], HIMAWARI8) self.assertEqual(res.attrs['satellite_longitude'], 140.7) + self.assertEqual(res.attrs['satellite_longitude'], 140.7) + self.assertEqual(res.attrs['satellite_altitude'], 35785831.0) - def test_get_area_def(self): - """Test getting an AreaDefinition.""" - from satpy import DatasetID - area_def = self.reader.get_area_def(DatasetID(name='B03', calibration='reflectance')) - self.assertTupleEqual(area_def.area_extent, - (-5499495.117842725, -16499485.352640428, 5500495.116954979, -5499495.117842725)) + # Check called methods + with mock.patch.object(reader, '_mask_space') as mask_space: + with mock.patch.object(reader, 'calibrate') as calibrate: + reader.get_dataset(key, {'units': '%', 'sensor': 'ahi'}) + mask_space.assert_called() + calibrate.assert_called() + + with mock.patch('logging.Logger.error') as log_mock: + reader.get_dataset(key, {'units': '%', 'sensor': 'jami'}) + log_mock.assert_called() def suite(): From c8e08085faa19c8aee28d9540e78a5b752ac51b0 Mon Sep 17 00:00:00 2001 From: Stephan Finkensieper Date: Fri, 30 Nov 2018 19:27:13 +0100 Subject: [PATCH 09/12] Add test for get_geostationary_mask --- satpy/tests/reader_tests/test_utils.py | 42 ++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/satpy/tests/reader_tests/test_utils.py b/satpy/tests/reader_tests/test_utils.py index ada95beaeb..06df560e61 100644 --- a/satpy/tests/reader_tests/test_utils.py +++ b/satpy/tests/reader_tests/test_utils.py @@ -30,6 +30,7 @@ import mock import numpy as np +import pyresample.geometry from satpy.readers import utils as hf @@ -316,6 +317,47 @@ def test_get_geostationary_angle_extent(self): np.testing.assert_allclose(expected, hf.get_geostationary_angle_extent(geos_area)) + def test_geostationary_mask(self): + """Test geostationary mask""" + # Compute mask of a very elliptical earth + area = pyresample.geometry.AreaDefinition( + area_id='FLDK', + name='Full Disk', + proj_id='geos', + proj_dict={'a': '6378169.0', + 'b': '3000000.0', + 'h': '35785831.0', + 'lon_0': '145.0', + 'proj': 'geos', + 'units': 'm'}, + x_size=101, + y_size=101, + area_extent=(-6498000.088960204, -6498000.088960204, + 6502000.089024927, 6502000.089024927)) + + mask = hf.get_geostationary_mask(area).astype(np.int).compute() + + # Check results along a couple of lines + # a) Horizontal + self.assertTrue(np.all(mask[50, :8] == 0)) + self.assertTrue(np.all(mask[50, 8:93] == 1)) + self.assertTrue(np.all(mask[50, 93:] == 0)) + + # b) Vertical + self.assertTrue(np.all(mask[:31, 50] == 0)) + self.assertTrue(np.all(mask[31:70, 50] == 1)) + self.assertTrue(np.all(mask[70:, 50] == 0)) + + # c) Top left to bottom right + self.assertTrue(np.all(mask[range(33), range(33)] == 0)) + self.assertTrue(np.all(mask[range(33, 68), range(33, 68)] == 1)) + self.assertTrue(np.all(mask[range(68, 101), range(68, 101)] == 0)) + + # d) Bottom left to top right + self.assertTrue(np.all(mask[range(101-1, 68-1, -1), range(33)] == 0)) + self.assertTrue(np.all(mask[range(68-1, 33-1, -1), range(33, 68)] == 1)) + self.assertTrue(np.all(mask[range(33-1, -1, -1), range(68, 101)] == 0)) + @mock.patch('satpy.readers.utils.AreaDefinition') def test_sub_area(self, adef): """Sub area slicing.""" From 0c4be48bc6f2dcc86864aa8fff59313e704b7641 Mon Sep 17 00:00:00 2001 From: Stephan Finkensieper Date: Fri, 30 Nov 2018 19:33:00 +0100 Subject: [PATCH 10/12] Fix stickler suggestions --- satpy/readers/utils.py | 1 - satpy/tests/reader_tests/test_hrit_jma.py | 3 +-- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/satpy/readers/utils.py b/satpy/readers/utils.py index ebbcd6cce5..f0261fb76b 100644 --- a/satpy/readers/utils.py +++ b/satpy/readers/utils.py @@ -26,7 +26,6 @@ import logging from contextlib import closing -import dask.array as da import tempfile import bz2 import os diff --git a/satpy/tests/reader_tests/test_hrit_jma.py b/satpy/tests/reader_tests/test_hrit_jma.py index 10a39ec8af..302a119adb 100644 --- a/satpy/tests/reader_tests/test_hrit_jma.py +++ b/satpy/tests/reader_tests/test_hrit_jma.py @@ -69,8 +69,7 @@ def _get_mda(self, loff=5500.0, coff=5500.0, nlines=11000, ncols=11000, 'loff': loff, 'number_of_columns': ncols, 'number_of_lines': nlines, - 'image_data_function': idf - } + 'image_data_function': idf} def test_init(self): """Test creating the file handler.""" From 20e8d324932e7122d32cd2af0df7c28ba75ed0f8 Mon Sep 17 00:00:00 2001 From: Stephan Finkensieper Date: Mon, 3 Dec 2018 10:58:31 +0000 Subject: [PATCH 11/12] Add hrit_mtsat* to the documentation --- doc/source/index.rst | 3 +++ 1 file changed, 3 insertions(+) diff --git a/doc/source/index.rst b/doc/source/index.rst index 94baf897bb..77e756325b 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -66,6 +66,9 @@ installation. * - Himawari 8 and 9 AHI data in HRIT format - `hrit_jma` - Nominal + * - MTSAT-1R/2 data in JMA HRIT format + - `hrit_mtsat1`, `hrit_mtsat2` + - Beta * - GOES 16 imager data in netcdf format - `abi_l1b` - Nominal From 19a304cc40a17d6983c345e0a3f40ff32d1e730b Mon Sep 17 00:00:00 2001 From: Stephan Finkensieper Date: Mon, 3 Dec 2018 16:25:41 +0000 Subject: [PATCH 12/12] Rename readers according to #527 --- doc/source/index.rst | 7 +++++-- satpy/etc/readers/{hrit_ahi.yaml => ahi_hrit.yaml} | 2 +- satpy/etc/readers/{hrit_mtsat1.yaml => jami_hrit.yaml} | 2 +- .../readers/{hrit_mtsat2.yaml => mtsat2-imager_hrit.yaml} | 2 +- 4 files changed, 8 insertions(+), 5 deletions(-) rename satpy/etc/readers/{hrit_ahi.yaml => ahi_hrit.yaml} (99%) rename satpy/etc/readers/{hrit_mtsat1.yaml => jami_hrit.yaml} (99%) rename satpy/etc/readers/{hrit_mtsat2.yaml => mtsat2-imager_hrit.yaml} (99%) diff --git a/doc/source/index.rst b/doc/source/index.rst index 77e756325b..4569c31b17 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -66,8 +66,11 @@ installation. * - Himawari 8 and 9 AHI data in HRIT format - `hrit_jma` - Nominal - * - MTSAT-1R/2 data in JMA HRIT format - - `hrit_mtsat1`, `hrit_mtsat2` + * - MTSAT-1R JAMI data in JMA HRIT format + - `jami_hrit` + - Beta + * - MTSAT-2 Imager data in JMA HRIT format + - `mtsat2-imager_hrit` - Beta * - GOES 16 imager data in netcdf format - `abi_l1b` diff --git a/satpy/etc/readers/hrit_ahi.yaml b/satpy/etc/readers/ahi_hrit.yaml similarity index 99% rename from satpy/etc/readers/hrit_ahi.yaml rename to satpy/etc/readers/ahi_hrit.yaml index 2d1e79c611..b27388d822 100644 --- a/satpy/etc/readers/hrit_ahi.yaml +++ b/satpy/etc/readers/ahi_hrit.yaml @@ -4,7 +4,7 @@ reader: description: JMA HRIT Reader - name: hrit_ahi + name: ahi_hrit sensors: [ahi] default_channels: [] reader: !!python/name:satpy.readers.yaml_reader.FileYAMLReader diff --git a/satpy/etc/readers/hrit_mtsat1.yaml b/satpy/etc/readers/jami_hrit.yaml similarity index 99% rename from satpy/etc/readers/hrit_mtsat1.yaml rename to satpy/etc/readers/jami_hrit.yaml index f3f81bcf5b..83279cd1bc 100644 --- a/satpy/etc/readers/hrit_mtsat1.yaml +++ b/satpy/etc/readers/jami_hrit.yaml @@ -8,7 +8,7 @@ reader: description: Reader for MTSAT-1R JAMI data in JMA HRIT Format - name: hrit_mtsat1 + name: jami_hrit sensors: [jami] default_channels: [] reader: !!python/name:satpy.readers.yaml_reader.FileYAMLReader diff --git a/satpy/etc/readers/hrit_mtsat2.yaml b/satpy/etc/readers/mtsat2-imager_hrit.yaml similarity index 99% rename from satpy/etc/readers/hrit_mtsat2.yaml rename to satpy/etc/readers/mtsat2-imager_hrit.yaml index 015bc81341..6d06161dca 100644 --- a/satpy/etc/readers/hrit_mtsat2.yaml +++ b/satpy/etc/readers/mtsat2-imager_hrit.yaml @@ -8,7 +8,7 @@ reader: description: Reader for MTSAT-2 data in JMA HRIT Format - name: hrit_mtsat2 + name: mtsat2-imager_hrit sensors: [mtsat2_imager] default_channels: [] reader: !!python/name:satpy.readers.yaml_reader.FileYAMLReader