diff --git a/doc/source/index.rst b/doc/source/index.rst index 2b66328518..a25a0396e3 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -215,6 +215,9 @@ the base Satpy installation. - `hsaf_grib` - | Beta | Only the h03, h03b, h05 and h05B products are supported at-present + * - GEO-KOMPSAT-2 AMI L1B data in NetCDF4 format + - `ami_l1b` + - Beta Indices and tables diff --git a/satpy/composites/__init__.py b/satpy/composites/__init__.py index f5a94ccc50..33f9464f08 100644 --- a/satpy/composites/__init__.py +++ b/satpy/composites/__init__.py @@ -761,6 +761,7 @@ def __init__(self, name, common_channel_mask=True, **kwargs): Args: common_channel_mask (bool): If True, mask all the channels with a mask that combines all the invalid areas of the given data. + """ self.common_channel_mask = common_channel_mask super(GenericCompositor, self).__init__(name, **kwargs) @@ -974,6 +975,7 @@ def __init__(self, name, lim_low=85., lim_high=95., **kwargs): blending of the given channels lim_high (float): upper limit of Sun zenith angle for the blending of the given channels + """ self.lim_low = lim_low self.lim_high = lim_high @@ -1172,7 +1174,7 @@ class RatioSharpenedRGB(GenericCompositor): footprint. Note that the input data to this compositor must already be resampled so all data arrays are the same shape. - Example: + Example:: R_lo - 1000m resolution - shape=(2000, 2000) G - 1000m resolution - shape=(2000, 2000) @@ -1293,7 +1295,7 @@ def _mean4(data, offset=(0, 0), block_id=None): class SelfSharpenedRGB(RatioSharpenedRGB): """Sharpen RGB with ratio of a band with a strided-version of itself. - Example: + Example:: R - 500m resolution - shape=(4000, 4000) G - 1000m resolution - shape=(2000, 2000) @@ -1402,6 +1404,7 @@ def __init__(self, name, filename=None, area=None, **kwargs): filename (str): Filename of the image to load area (str): Name of area definition for the image. Optional for images with built-in area definitions (geotiff) + """ if filename is None: raise ValueError("No image configured for static image compositor") @@ -1417,9 +1420,9 @@ def __call__(self, *args, **kwargs): """Call the compositor.""" from satpy import Scene # Check if filename exists, if not then try from SATPY_ANCPATH - if (not os.path.isfile(self.filename)): + if not os.path.isfile(self.filename): tmp_filename = os.path.join(get_environ_ancpath(), self.filename) - if (os.path.isfile(tmp_filename)): + if os.path.isfile(tmp_filename): self.filename = tmp_filename scn = Scene(reader='generic_image', filenames=[self.filename]) scn.load(['image']) diff --git a/satpy/composites/ahi.py b/satpy/composites/ahi.py index ab069609b2..c170eb0543 100644 --- a/satpy/composites/ahi.py +++ b/satpy/composites/ahi.py @@ -15,8 +15,7 @@ # # You should have received a copy of the GNU General Public License along with # satpy. If not, see . -"""Composite classes for the AHI instrument. -""" +"""Composite classes for the AHI instrument.""" import logging @@ -26,17 +25,19 @@ class GreenCorrector(GenericCompositor): - """Corrector of the AHI green band to compensate for the deficit of - chlorophyl signal. - """ + """Corrector of the AHI green band to compensate for the deficit of chlorophyll signal.""" + + def __init__(self, *args, **kwargs): + """Set default keyword argument values.""" + # XXX: Should this be 0.93 and 0.07 + self.fractions = kwargs.pop('fractions', [0.85, 0.15]) + super(GreenCorrector, self).__init__(*args, **kwargs) def __call__(self, projectables, optional_datasets=None, **attrs): """Boost vegetation effect thanks to NIR (0.8µm) band.""" - green, nir = self.match_data_arrays(projectables) LOG.info('Boosting vegetation on green band') - # XXX: Should this be 0.93 and 0.07 - new_green = green * 0.85 + nir * 0.15 + new_green = green * self.fractions[0] + nir * self.fractions[1] new_green.attrs = green.attrs.copy() return super(GreenCorrector, self).__call__((new_green,), **attrs) diff --git a/satpy/etc/composites/ami.yaml b/satpy/etc/composites/ami.yaml new file mode 100644 index 0000000000..ae670df403 --- /dev/null +++ b/satpy/etc/composites/ami.yaml @@ -0,0 +1,153 @@ +sensor_name: visir/ami + +composites: + green_raw: + compositor: !!python/name:satpy.composites.ahi.GreenCorrector + prerequisites: + - name: VI005 + modifiers: [sunz_corrected] + - name: VI008 + modifiers: [sunz_corrected] + standard_name: toa_bidirectional_reflectance + fractions: [0.85, 0.15] + + green: + compositor: !!python/name:satpy.composites.ahi.GreenCorrector + prerequisites: + - name: VI005 + modifiers: [sunz_corrected, rayleigh_corrected] + - name: VI008 + modifiers: [sunz_corrected] + standard_name: toa_bidirectional_reflectance + fractions: [0.85, 0.15] + + true_color_raw: + compositor: !!python/name:satpy.composites.SelfSharpenedRGB + prerequisites: + - name: VI006 + modifiers: [sunz_corrected] + - name: green_raw + - name: VI004 + modifiers: [sunz_corrected] + standard_name: true_color + + true_color: + compositor: !!python/name:satpy.composites.SelfSharpenedRGB + prerequisites: + - name: VI006 + modifiers: [sunz_corrected, rayleigh_corrected] + - name: green + - name: VI004 + modifiers: [sunz_corrected, rayleigh_corrected] + standard_name: true_color + + overview: + compositor: !!python/name:satpy.composites.GenericCompositor + prerequisites: + - 0.65 + - 0.85 + - 10.4 + standard_name: overview + + natural_color: + compositor: !!python/name:satpy.composites.SelfSharpenedRGB + prerequisites: + - name: NR016 + modifiers: [sunz_corrected] #, rayleigh_corrected] + - name: VI008 + modifiers: [sunz_corrected] #, rayleigh_corrected] + - name: VI006 + modifiers: [sunz_corrected] #, rayleigh_corrected] + high_resolution_band: blue + standard_name: natural_color + + day_microphysics_eum: + compositor: !!python/name:satpy.composites.GenericCompositor + prerequisites: + - wavelength: 0.86 + - wavelength: 3.9 + modifiers: [nir_reflectance] + - wavelength: 10.4 + standard_name: day_microphysics + + cloud_phase_distinction: + compositor: !!python/name:satpy.composites.GenericCompositor + prerequisites: + - wavelength: 10.4 + - wavelength: 0.64 + - wavelength: 1.6 + standard_name: cloud_phase_distinction + + water_vapors1: + compositor: !!python/name:satpy.composites.GenericCompositor + prerequisites: + - wavelength: 10.4 + - wavelength: 6.2 + - wavelength: 7.3 + standard_name: water_vapors1 + + mid_vapor: + compositor: !!python/name:satpy.composites.DifferenceCompositor + prerequisites: + - wavelength: 7.3 + - wavelength: 6.2 + standard_name: mid_vapor + + water_vapors2: + compositor: !!python/name:satpy.composites.GenericCompositor + prerequisites: + - name: mid_vapor + - wavelength: 7.3 + - wavelength: 6.2 + standard_name: water_vapors2 + + convection: + compositor: !!python/name:satpy.composites.GenericCompositor + prerequisites: + - compositor: !!python/name:satpy.composites.DifferenceCompositor + prerequisites: + - WV069 + - WV073 + - compositor: !!python/name:satpy.composites.DifferenceCompositor + prerequisites: + - SW038 + - IR105 + - compositor: !!python/name:satpy.composites.DifferenceCompositor + prerequisites: + - NR016 + - VI006 + standard_name: convection + + ir_cloud_day: + standard_name: ir_cloud_day + compositor: !!python/name:satpy.composites.CloudCompositor + prerequisites: + - name: IR112 + + airmass: + compositor: !!python/name:satpy.composites.GenericCompositor + prerequisites: + - compositor: !!python/name:satpy.composites.DifferenceCompositor + prerequisites: + - name: WV063 + - name: WV073 + - compositor: !!python/name:satpy.composites.DifferenceCompositor + prerequisites: + - name: IR096 + - name: IR105 + - name: WV063 + standard_name: airmass + + ash: + compositor: !!python/name:satpy.composites.GenericCompositor + prerequisites: + - compositor: !!python/name:satpy.composites.DifferenceCompositor + prerequisites: + - IR123 + - IR112 + - compositor: !!python/name:satpy.composites.DifferenceCompositor + prerequisites: + - IR112 + - IR087 + - IR112 + standard_name: ash diff --git a/satpy/etc/readers/ami_l1b.yaml b/satpy/etc/readers/ami_l1b.yaml new file mode 100644 index 0000000000..b0311e7b35 --- /dev/null +++ b/satpy/etc/readers/ami_l1b.yaml @@ -0,0 +1,341 @@ +reader: + name: ami_l1b + short_name: AMI L1b + long_name: GEO-KOMPSAT-2 AMI Level 1b + description: > + GEO-KOMPSAT-2 AMI Level 1b data reader in the NetCDF4 format. The file format and + instrument are described on KMA's website + `here `_. + sensors: [ami] + default_channels: + reader: !!python/name:satpy.readers.yaml_reader.FileYAMLReader + # file pattern keys to sort files by with 'satpy.utils.group_files' + group_keys: ['start_time', 'platform_shortname', 'sensor', 'sector_info'] + +file_types: + # Example: gk2a_ami_le1b_ir087_fd020ge_201901260310.nc + # Below list is alphabetical + ir087: + file_reader: !!python/name:satpy.readers.ami_l1b.AMIL1bNetCDF + file_patterns: ['{platform_shortname:4s}_{sensor:3s}_le1b_ir087_{sector_info:s}_{start_time:%Y%m%d%H%M}.nc'] + ir096: + file_reader: !!python/name:satpy.readers.ami_l1b.AMIL1bNetCDF + file_patterns: ['{platform_shortname:4s}_{sensor:3s}_le1b_ir096_{sector_info:s}_{start_time:%Y%m%d%H%M}.nc'] + ir105: + file_reader: !!python/name:satpy.readers.ami_l1b.AMIL1bNetCDF + file_patterns: ['{platform_shortname:4s}_{sensor:3s}_le1b_ir105_{sector_info:s}_{start_time:%Y%m%d%H%M}.nc'] + ir112: + file_reader: !!python/name:satpy.readers.ami_l1b.AMIL1bNetCDF + file_patterns: ['{platform_shortname:4s}_{sensor:3s}_le1b_ir112_{sector_info:s}_{start_time:%Y%m%d%H%M}.nc'] + ir123: + file_reader: !!python/name:satpy.readers.ami_l1b.AMIL1bNetCDF + file_patterns: ['{platform_shortname:4s}_{sensor:3s}_le1b_ir123_{sector_info:s}_{start_time:%Y%m%d%H%M}.nc'] + ir133: + file_reader: !!python/name:satpy.readers.ami_l1b.AMIL1bNetCDF + file_patterns: ['{platform_shortname:4s}_{sensor:3s}_le1b_ir133_{sector_info:s}_{start_time:%Y%m%d%H%M}.nc'] + nr013: + file_reader: !!python/name:satpy.readers.ami_l1b.AMIL1bNetCDF + file_patterns: ['{platform_shortname:4s}_{sensor:3s}_le1b_nr013_{sector_info:s}_{start_time:%Y%m%d%H%M}.nc'] + nr016: + file_reader: !!python/name:satpy.readers.ami_l1b.AMIL1bNetCDF + file_patterns: ['{platform_shortname:4s}_{sensor:3s}_le1b_nr016_{sector_info:s}_{start_time:%Y%m%d%H%M}.nc'] + sw038: + file_reader: !!python/name:satpy.readers.ami_l1b.AMIL1bNetCDF + file_patterns: ['{platform_shortname:4s}_{sensor:3s}_le1b_sw038_{sector_info:s}_{start_time:%Y%m%d%H%M}.nc'] + vi004: + file_reader: !!python/name:satpy.readers.ami_l1b.AMIL1bNetCDF + file_patterns: ['{platform_shortname:4s}_{sensor:3s}_le1b_vi004_{sector_info:s}_{start_time:%Y%m%d%H%M}.nc'] + vi005: + file_reader: !!python/name:satpy.readers.ami_l1b.AMIL1bNetCDF + file_patterns: ['{platform_shortname:4s}_{sensor:3s}_le1b_vi005_{sector_info:s}_{start_time:%Y%m%d%H%M}.nc'] + vi006: + file_reader: !!python/name:satpy.readers.ami_l1b.AMIL1bNetCDF + file_patterns: ['{platform_shortname:4s}_{sensor:3s}_le1b_vi006_{sector_info:s}_{start_time:%Y%m%d%H%M}.nc'] + vi008: + file_reader: !!python/name:satpy.readers.ami_l1b.AMIL1bNetCDF + file_patterns: ['{platform_shortname:4s}_{sensor:3s}_le1b_vi008_{sector_info:s}_{start_time:%Y%m%d%H%M}.nc'] + wv063: + file_reader: !!python/name:satpy.readers.ami_l1b.AMIL1bNetCDF + file_patterns: ['{platform_shortname:4s}_{sensor:3s}_le1b_wv063_{sector_info:s}_{start_time:%Y%m%d%H%M}.nc'] + wv069: + file_reader: !!python/name:satpy.readers.ami_l1b.AMIL1bNetCDF + file_patterns: ['{platform_shortname:4s}_{sensor:3s}_le1b_wv069_{sector_info:s}_{start_time:%Y%m%d%H%M}.nc'] + wv073: + file_reader: !!python/name:satpy.readers.ami_l1b.AMIL1bNetCDF + file_patterns: ['{platform_shortname:4s}_{sensor:3s}_le1b_wv073_{sector_info:s}_{start_time:%Y%m%d%H%M}.nc'] + +datasets: + # Below list is ordered the same as the table: + # https://directory.eoportal.org/web/eoportal/satellite-missions/content/-/article/geo-kompsat-2 + C01: + name: VI004 + wavelength: [0.450, 0.470, 0.490] + resolution: 1000 + calibration: + counts: + standard_name: counts + units: 1 + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavelength + units: W m-2 um-1 sr-1 + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + file_type: vi004 + file_key: image_pixel_values + + C02: + name: VI005 + wavelength: [0.495, 0.509, 0.523] + resolution: 1000 + calibration: + counts: + standard_name: counts + units: 1 + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavelength + units: W m-2 um-1 sr-1 + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + file_type: vi005 + file_key: image_pixel_values + + C03: + name: VI006 + wavelength: [0.599, 0.639, 0.679] + resolution: 500 + calibration: + counts: + standard_name: counts + units: 1 + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavelength + units: W m-2 um-1 sr-1 + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + file_type: vi006 + file_key: image_pixel_values + + C04: + name: VI008 + wavelength: [0.846, 0.863, 0.880] + resolution: 1000 + calibration: + counts: + standard_name: counts + units: 1 + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavelength + units: W m-2 um-1 sr-1 + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + file_type: vi008 + file_key: image_pixel_values + + C05: + name: NR013 + wavelength: [1.363, 1.37, 1.377] + resolution: 2000 + calibration: + counts: + standard_name: counts + units: 1 + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavelength + units: W m-2 um-1 sr-1 + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + file_type: nr013 + file_key: image_pixel_values + + C06: + name: NR016 + wavelength: [1.590, 1.61, 1.630] + resolution: 2000 + calibration: + counts: + standard_name: counts + units: 1 + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavelength + units: W m-2 um-1 sr-1 + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + file_type: nr016 + file_key: image_pixel_values + + C07: + name: SW038 + wavelength: [3.74, 3.83, 3.92] + resolution: 2000 + calibration: + counts: + standard_name: counts + units: 1 + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: K + file_type: sw038 + file_key: image_pixel_values + + C08: + name: WV063 + wavelength: [5.79, 6.21, 6.63] + resolution: 2000 + calibration: + counts: + standard_name: counts + units: 1 + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: K + file_type: wv063 + file_key: image_pixel_values + + C09: + name: WV069 + wavelength: [6.74, 6.94, 7.14] + resolution: 2000 + calibration: + counts: + standard_name: counts + units: 1 + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: K + file_type: wv069 + file_key: image_pixel_values + + C10: + name: WV073 + wavelength: [7.24, 7.33, 7.42] + resolution: 2000 + calibration: + counts: + standard_name: counts + units: 1 + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: K + file_type: wv073 + file_key: image_pixel_values + + C11: + name: IR087 + wavelength: [8.415, 8.59, 8.765] + resolution: 2000 + calibration: + counts: + standard_name: counts + units: 1 + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: K + file_type: ir087 + file_key: image_pixel_values + + C12: + name: IR096 + wavelength: [9.43, 9.62, 9.81] + resolution: 2000 + calibration: + counts: + standard_name: counts + units: 1 + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: K + file_type: ir096 + file_key: image_pixel_values + + C13: + name: IR105 + wavelength: [10.115, 10.35, 10.585] + resolution: 2000 + calibration: + counts: + standard_name: counts + units: 1 + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: K + file_type: ir105 + file_key: image_pixel_values + + C14: + name: IR112 + wavelength: [10.90, 11.23, 11.56] + resolution: 2000 + calibration: + counts: + standard_name: counts + units: 1 + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: K + file_type: ir112 + file_key: image_pixel_values + + C15: + name: IR123 + wavelength: [11.805, 12.36, 12.915] + resolution: 2000 + calibration: + counts: + standard_name: counts + units: 1 + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: K + file_type: ir123 + file_key: image_pixel_values + + C16: + name: IR133 + wavelength: [13.005, 13.29, 13.575] + resolution: 2000 + calibration: + counts: + standard_name: counts + units: 1 + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: K + file_type: ir133 + file_key: image_pixel_values + diff --git a/satpy/readers/ami_l1b.py b/satpy/readers/ami_l1b.py new file mode 100644 index 0000000000..7867d679d2 --- /dev/null +++ b/satpy/readers/ami_l1b.py @@ -0,0 +1,230 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Copyright (c) 2019 Satpy developers +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +"""Advanced Meteorological Imager reader for the Level 1b NetCDF4 format.""" + +import logging +from datetime import datetime, timedelta + +import numpy as np +import xarray as xr +import dask.array as da +import pyproj + +from pyresample import geometry +from pyspectral.blackbody import blackbody_wn_rad2temp as rad2temp +from satpy.readers.file_handlers import BaseFileHandler +from satpy import CHUNK_SIZE + +logger = logging.getLogger(__name__) + +PLATFORM_NAMES = { + 'GK-2A': 'GEO-KOMPSAT-2A', + 'GK-2B': 'GEO-KOMPSAT-2B', +} + + +class AMIL1bNetCDF(BaseFileHandler): + """Base reader for AMI L1B NetCDF4 files.""" + + def __init__(self, filename, filename_info, filetype_info, + calib_mode='PYSPECTRAL', allow_conditional_pixels=False): + """Open the NetCDF file with xarray and prepare the Dataset for reading.""" + super(AMIL1bNetCDF, self).__init__(filename, filename_info, filetype_info) + self.nc = xr.open_dataset(self.filename, + decode_cf=True, + mask_and_scale=False, + chunks={'dim_image_x': CHUNK_SIZE, 'dim_image_y': CHUNK_SIZE}) + self.nc = self.nc.rename({'dim_image_x': 'x', 'dim_image_y': 'y'}) + + platform_shortname = self.nc.attrs['satellite_name'] + self.platform_name = PLATFORM_NAMES.get(platform_shortname) + self.sensor = 'ami' + self.allow_conditional_pixels = allow_conditional_pixels + calib_mode_choices = ('FILE', 'PYSPECTRAL') + if calib_mode.upper() not in calib_mode_choices: + raise ValueError('Invalid calibration mode: {}. Choose one of {}'.format( + calib_mode, calib_mode_choices)) + self.calib_mode = calib_mode.upper() + + @property + def start_time(self): + """Get observation start time.""" + base = datetime(2000, 1, 1, 12, 0, 0) + return base + timedelta(seconds=self.nc.attrs['observation_start_time']) + + @property + def end_time(self): + """Get observation end time.""" + base = datetime(2000, 1, 1, 12, 0, 0) + return base + timedelta(seconds=self.nc.attrs['observation_end_time']) + + def get_area_def(self, dsid): + """Get area definition for this file.""" + a = self.nc.attrs['earth_equatorial_radius'] + b = self.nc.attrs['earth_polar_radius'] + h = self.nc.attrs['nominal_satellite_height'] - a + lon_0 = self.nc.attrs['sub_longitude'] * 180 / np.pi # it's in radians? + cols = self.nc.attrs['number_of_columns'] + rows = self.nc.attrs['number_of_lines'] + obs_mode = self.nc.attrs['observation_mode'] + resolution = self.nc.attrs['channel_spatial_resolution'] + + cfac = self.nc.attrs['cfac'] + coff = self.nc.attrs['coff'] + lfac = self.nc.attrs['lfac'] + loff = self.nc.attrs['loff'] + bit_shift = 2**16 + area_extent = ( + h * np.deg2rad((0 - coff - 0.5) * bit_shift / cfac), + -h * np.deg2rad((0 - loff - 0.5) * bit_shift / lfac), + h * np.deg2rad((cols - coff + 0.5) * bit_shift / cfac), + -h * np.deg2rad((rows - loff + 0.5) * bit_shift / lfac), + ) + + proj_dict = { + 'proj': 'geos', + 'lon_0': float(lon_0), + 'a': float(a), + 'b': float(b), + 'h': h, + 'units': 'm' + } + + fg_area_def = geometry.AreaDefinition( + 'ami_geos_{}'.format(obs_mode.lower()), + 'AMI {} Area at {} resolution'.format(obs_mode, resolution), + 'ami_fixed_grid', + proj_dict, + cols, + rows, + np.asarray(area_extent)) + + return fg_area_def + + def get_orbital_parameters(self): + """Collect orbital parameters for this file.""" + a = float(self.nc.attrs['earth_equatorial_radius']) + b = float(self.nc.attrs['earth_polar_radius']) + # nominal_satellite_height seems to be from the center of the earth + h = float(self.nc.attrs['nominal_satellite_height']) - a + lon_0 = self.nc.attrs['sub_longitude'] * 180 / np.pi # it's in radians? + sc_position = self.nc['sc_position'].attrs['sc_position_center_pixel'] + + # convert ECEF coordinates to lon, lat, alt + ecef = pyproj.Proj(proj='geocent', a=a, b=b) + lla = pyproj.Proj(proj='latlong', a=a, b=b) + sc_position = pyproj.transform( + ecef, lla, sc_position[0], sc_position[1], sc_position[2]) + + orbital_parameters = { + 'projection_longitude': float(lon_0), + 'projection_latitude': 0.0, + 'projection_altitude': h, + 'satellite_actual_longitude': sc_position[0], + 'satellite_actual_latitude': sc_position[1], + 'satellite_actual_altitude': sc_position[2] / 1000.0, # km + } + return orbital_parameters + + def get_dataset(self, dataset_id, ds_info): + """Load a dataset as a xarray DataArray.""" + file_key = ds_info.get('file_key', dataset_id.name) + data = self.nc[file_key] + # hold on to attributes for later + attrs = data.attrs + # highest 2 bits are data quality flags + # 00=no error + # 01=available under conditions + # 10=outside the viewing area + # 11=Error exists + if self.allow_conditional_pixels: + qf = data & 0b1000000000000000 + else: + qf = data & 0b1100000000000000 + + # mask DQF bits + bits = attrs['number_of_valid_bits_per_pixel'] + data &= 2**bits - 1 + # only take "no error" pixels as valid + data = data.where(qf == 0) + + # Calibration values from file, fall back to built-in if unavailable + gain = self.nc.attrs['DN_to_Radiance_Gain'] + offset = self.nc.attrs['DN_to_Radiance_Offset'] + + if dataset_id.calibration in ('radiance', 'reflectance', 'brightness_temperature'): + data = gain * data + offset + if dataset_id.calibration == 'reflectance': + # depends on the radiance calibration above + rad_to_alb = self.nc.attrs['Radiance_to_Albedo_c'] + if ds_info.get('units') == '%': + rad_to_alb *= 100 + data = data * rad_to_alb + elif dataset_id.calibration == 'brightness_temperature': + data = self._calibrate_ir(dataset_id, data) + elif dataset_id.calibration not in ('counts', 'radiance'): + raise ValueError("Unknown calibration: '{}'".format(dataset_id.calibration)) + + for attr_name in ('standard_name', 'units'): + attrs[attr_name] = ds_info[attr_name] + attrs.update(dataset_id.to_dict()) + attrs['orbital_parameters'] = self.get_orbital_parameters() + attrs['platform_name'] = self.platform_name + attrs['sensor'] = self.sensor + data.attrs = attrs + return data + + def _calibrate_ir(self, dataset_id, data): + """Calibrate radiance data to BTs using either pyspectral or in-file coefficients.""" + if self.calib_mode == 'PYSPECTRAL': + # depends on the radiance calibration above + # Convert um to m^-1 (SI units for pyspectral) + wn = 1 / (dataset_id.wavelength[1] / 1e6) + # Convert cm^-1 (wavenumbers) and (mW/m^2)/(str/cm^-1) (radiance data) + # to SI units m^-1, mW*m^-3*str^-1. + bt_data = rad2temp(wn, data.data * 1e-5) + if isinstance(bt_data, np.ndarray): + # old versions of pyspectral produce numpy arrays + data.data = da.from_array(bt_data, chunks=data.data.chunks) + else: + # new versions of pyspectral can do dask arrays + data.data = bt_data + else: + # IR coefficients from the file + # Channel specific + c0 = self.nc.attrs['Teff_to_Tbb_c0'] + c1 = self.nc.attrs['Teff_to_Tbb_c1'] + c2 = self.nc.attrs['Teff_to_Tbb_c2'] + + # These should be fixed, but load anyway + cval = self.nc.attrs['light_speed'] + kval = self.nc.attrs['Boltzmann_constant_k'] + hval = self.nc.attrs['Plank_constant_h'] + + # Compute wavenumber as cm-1 + wn = (10000 / dataset_id.wavelength[1]) * 100 + + # Convert radiance to effective brightness temperature + e1 = (2 * hval * cval * cval) * np.power(wn, 3) + e2 = (data.data * 1e-5) + t_eff = ((hval * cval / kval) * wn) / np.log((e1 / e2) + 1) + + # Now convert to actual brightness temperature + bt_data = c0 + c1 * t_eff + c2 * t_eff * t_eff + data.data = bt_data + return data diff --git a/satpy/tests/reader_tests/__init__.py b/satpy/tests/reader_tests/__init__.py index ccb1062ac0..ce2dcde3d1 100644 --- a/satpy/tests/reader_tests/__init__.py +++ b/satpy/tests/reader_tests/__init__.py @@ -38,7 +38,8 @@ test_electrol_hrit, test_mersi2_l1b, test_avhrr_l1b_gaclac, test_vaisala_gld360, test_fci_l1c_fdhsi, test_tropomi_l2, - test_hsaf_grib, test_abi_l2_nc, test_eum_base) + test_hsaf_grib, test_abi_l2_nc, test_eum_base, + test_ami_l1b) if sys.version_info < (2, 7): import unittest2 as unittest @@ -95,5 +96,6 @@ def suite(): mysuite.addTests(test_tropomi_l2.suite()) mysuite.addTests(test_hsaf_grib.suite()) mysuite.addTests(test_eum_base.suite()) + mysuite.addTests(test_ami_l1b.suite()) return mysuite diff --git a/satpy/tests/reader_tests/test_ami_l1b.py b/satpy/tests/reader_tests/test_ami_l1b.py new file mode 100644 index 0000000000..ff2667fb90 --- /dev/null +++ b/satpy/tests/reader_tests/test_ami_l1b.py @@ -0,0 +1,313 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# Copyright (c) 2019 Satpy developers +# +# This file is part of satpy. +# +# satpy is free software: you can redistribute it and/or modify it under the +# terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. +# +# satpy is distributed in the hope that it will be useful, but WITHOUT ANY +# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +# A PARTICULAR PURPOSE. See the GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License along with +# satpy. If not, see . +"""The ami_l1b reader tests package.""" + +import numpy as np +import xarray as xr +import dask.array as da + +import unittest + +try: + from unittest import mock +except ImportError: + import mock + + +class FakeDataset(object): + """Mimic xarray Dataset object.""" + + def __init__(self, info, attrs): + """Initialize test data.""" + for var_name, var_data in list(info.items()): + if isinstance(var_data, np.ndarray): + info[var_name] = xr.DataArray(var_data) + self.info = info + self.attrs = attrs + + def __getitem__(self, key): + """Mimic getitem method.""" + return self.info[key] + + def __contains__(self, key): + """Mimic contains method.""" + return key in self.info + + def rename(self, *args, **kwargs): + """Mimic rename method.""" + return self + + def close(self): + """Act like close method.""" + return + + +class TestAMIL1bNetCDFBase(unittest.TestCase): + """Common setup for NC_ABI_L1B tests.""" + + @mock.patch('satpy.readers.ami_l1b.xr') + def setUp(self, xr_, counts=None): + """Create a fake dataset using the given counts data.""" + from satpy.readers.ami_l1b import AMIL1bNetCDF + + if counts is None: + rad_data = (np.arange(10.).reshape((2, 5)) + 1.) * 50. + rad_data = (rad_data + 1.) / 0.5 + rad_data = rad_data.astype(np.int16) + counts = xr.DataArray( + da.from_array(rad_data, chunks='auto'), + dims=('y', 'x'), + attrs={ + 'channel_name': "VI006", + 'detector_side': 2, + 'number_of_total_pixels': 484000000, + 'number_of_error_pixels': 113892451, + 'max_pixel_value': 32768, + 'min_pixel_value': 6, + 'average_pixel_value': 8228.98770845248, + 'stddev_pixel_value': 13621.130386551, + 'number_of_total_bits_per_pixel': 16, + 'number_of_data_quality_flag_bits_per_pixel': 2, + 'number_of_valid_bits_per_pixel': 12, + 'data_quality_flag_meaning': + "0:good_pixel, 1:conditionally_usable_pixel, 2:out_of_scan_area_pixel, 3:error_pixel", + 'ground_sample_distance_ew': 1.4e-05, + 'ground_sample_distance_ns': 1.4e-05, + } + ) + sc_position = xr.DataArray(0., attrs={ + 'sc_position_center_pixel': [-26113466.1974016, 33100139.1630508, 3943.75470244799], + }) + xr_.open_dataset.return_value = FakeDataset( + { + 'image_pixel_values': counts, + 'sc_position': sc_position, + }, + { + "satellite_name": "GK-2A", + "observation_start_time": 623084431.957882, + "observation_end_time": 623084975.606133, + "projection_type": "GEOS", + "sub_longitude": 2.23751210105673, + "cfac": 81701355.6133574, + "lfac": -81701355.6133574, + "coff": 11000.5, + "loff": 11000.5, + "nominal_satellite_height": 42164000., + "earth_equatorial_radius": 6378137., + "earth_polar_radius": 6356752.3, + "number_of_columns": 22000, + "number_of_lines": 22000, + "observation_mode": "FD", + "channel_spatial_resolution": "0.5", + "Radiance_to_Albedo_c": 1, + "DN_to_Radiance_Gain": -0.0144806550815701, + "DN_to_Radiance_Offset": 118.050903320312, + "Teff_to_Tbb_c0": -0.141418528203155, + "Teff_to_Tbb_c1": 1.00052232906885, + "Teff_to_Tbb_c2": -0.00000036287276076109, + "light_speed": 2.9979245800E+08, + "Boltzmann_constant_k": 1.3806488000E-23, + "Plank_constant_h": 6.6260695700E-34, + } + ) + + self.reader = AMIL1bNetCDF('filename', + {'platform_shortname': 'gk2a'}, + {'filetype': 'info'}) + + +class TestAMIL1bNetCDF(TestAMIL1bNetCDFBase): + """Test the AMI L1b reader.""" + + def _check_orbital_parameters(self, orb_params): + """Check that orbital parameters match expected values.""" + exp_params = { + 'projection_altitude': 35785863.0, + 'projection_latitude': 0.0, + 'projection_longitude': 128.2, + 'satellite_actual_altitude': 35782.65456070405, + 'satellite_actual_latitude': 0.005364927, + 'satellite_actual_longitude': 128.2707, + } + for key, val in exp_params.items(): + self.assertAlmostEqual(val, orb_params[key], places=3) + + def test_basic_attributes(self): + """Test getting basic file attributes.""" + from datetime import datetime + self.assertEqual(self.reader.start_time, + datetime(2019, 9, 30, 3, 0, 31, 957882)) + self.assertEqual(self.reader.end_time, + datetime(2019, 9, 30, 3, 9, 35, 606133)) + + def test_get_dataset(self): + """Test gettting radiance data.""" + from satpy import DatasetID + key = DatasetID(name='VI006', calibration='radiance') + res = self.reader.get_dataset(key, { + 'file_key': 'image_pixel_values', + 'standard_name': 'toa_outgoing_radiance_per_unit_wavelength', + 'units': 'W m-2 um-1 sr-1', + }) + exp = {'calibration': 'radiance', + 'modifiers': (), + 'platform_name': 'GEO-KOMPSAT-2A', + 'sensor': 'ami', + 'units': 'W m-2 um-1 sr-1'} + for key, val in exp.items(): + self.assertEqual(val, res.attrs[key]) + self._check_orbital_parameters(res.attrs['orbital_parameters']) + + def test_bad_calibration(self): + """Test that asking for a bad calibration fails.""" + from satpy import DatasetID + self.assertRaises(ValueError, self.reader.get_dataset, + DatasetID(name='VI006', calibration='_bad_'), + {'file_key': 'image_pixel_values', + 'standard_name': 'toa_outgoing_radiance_per_unit_wavelength', + 'units': 'W m-2 um-1 sr-1', + }) + + @mock.patch('satpy.readers.abi_base.geometry.AreaDefinition') + def test_get_area_def(self, adef): + """Test the area generation.""" + self.reader.get_area_def(None) + + self.assertEqual(adef.call_count, 1) + call_args = tuple(adef.call_args)[0] + exp = {'a': 6378137.0, 'b': 6356752.3, 'h': 35785863.0, + 'lon_0': 128.2, 'proj': 'geos', 'units': 'm'} + for key, val in exp.items(): + self.assertIn(key, call_args[3]) + self.assertAlmostEqual(val, call_args[3][key]) + self.assertEqual(call_args[4], self.reader.nc.attrs['number_of_columns']) + self.assertEqual(call_args[5], self.reader.nc.attrs['number_of_lines']) + np.testing.assert_allclose(call_args[6], + [-5511523.904082, -5511523.904082, 5511022.902, 5511022.902]) + + def test_get_dataset_vis(self): + """Test get visible calibrated data.""" + from satpy import DatasetID + key = DatasetID(name='VI006', calibration='reflectance') + res = self.reader.get_dataset(key, { + 'file_key': 'image_pixel_values', + 'standard_name': 'toa_bidirectional_reflectance', + 'units': '%', + }) + exp = {'calibration': 'reflectance', + 'modifiers': (), + 'platform_name': 'GEO-KOMPSAT-2A', + 'sensor': 'ami', + 'units': '%'} + for key, val in exp.items(): + self.assertEqual(val, res.attrs[key]) + self._check_orbital_parameters(res.attrs['orbital_parameters']) + + def test_get_dataset_counts(self): + """Test get counts data.""" + from satpy import DatasetID + key = DatasetID(name='VI006', calibration='counts') + res = self.reader.get_dataset(key, { + 'file_key': 'image_pixel_values', + 'standard_name': 'counts', + 'units': '1', + }) + exp = {'calibration': 'counts', + 'modifiers': (), + 'platform_name': 'GEO-KOMPSAT-2A', + 'sensor': 'ami', + 'units': '1'} + for key, val in exp.items(): + self.assertEqual(val, res.attrs[key]) + self._check_orbital_parameters(res.attrs['orbital_parameters']) + + +class TestAMIL1bNetCDFIRCal(TestAMIL1bNetCDFBase): + """Test IR specific things about the AMI reader.""" + + def setUp(self): + """Create test data for IR calibration tests.""" + count_data = (np.arange(10).reshape((2, 5))) + 7000 + count_data = count_data.astype(np.uint16) + count = xr.DataArray( + da.from_array(count_data, chunks='auto'), + dims=('y', 'x'), + attrs={ + 'channel_name': "IR087", + 'detector_side': 2, + 'number_of_total_pixels': 484000000, + 'number_of_error_pixels': 113892451, + 'max_pixel_value': 32768, + 'min_pixel_value': 6, + 'average_pixel_value': 8228.98770845248, + 'stddev_pixel_value': 13621.130386551, + 'number_of_total_bits_per_pixel': 16, + 'number_of_data_quality_flag_bits_per_pixel': 2, + 'number_of_valid_bits_per_pixel': 13, + 'data_quality_flag_meaning': + "0:good_pixel, 1:conditionally_usable_pixel, 2:out_of_scan_area_pixel, 3:error_pixel", + 'ground_sample_distance_ew': 1.4e-05, + 'ground_sample_distance_ns': 1.4e-05, + } + ) + super(TestAMIL1bNetCDFIRCal, self).setUp(counts=count) + + def test_ir_calibrate(self): + """Test IR calibration.""" + from satpy import DatasetID + from satpy.readers.ami_l1b import rad2temp + ds_id = DatasetID(name='IR087', wavelength=[8.415, 8.59, 8.765], + calibration='brightness_temperature') + ds_info = { + 'file_key': 'image_pixel_values', + 'wavelength': [8.415, 8.59, 8.765], + 'standard_name': 'toa_brightness_temperature', + 'units': 'K', + } + with mock.patch('satpy.readers.ami_l1b.rad2temp', wraps=rad2temp) as r2t_mock: + res = self.reader.get_dataset(ds_id, ds_info) + r2t_mock.assert_called_once() + expected = np.array([[238.34385135, 238.31443527, 238.28500087, 238.25554813, 238.22607701], + [238.1965875, 238.16707956, 238.13755317, 238.10800829, 238.07844489]]) + np.testing.assert_allclose(res.data.compute(), expected, equal_nan=True) + # make sure the attributes from the file are in the data array + self.assertEqual(res.attrs['standard_name'], 'toa_brightness_temperature') + + # test builtin coefficients + self.reader.calib_mode = 'FILE' + with mock.patch('satpy.readers.ami_l1b.rad2temp', wraps=rad2temp) as r2t_mock: + res = self.reader.get_dataset(ds_id, ds_info) + r2t_mock.assert_not_called() + # file coefficients are pretty close, give some wiggle room + np.testing.assert_allclose(res.data.compute(), expected, equal_nan=True, atol=0.04) + # make sure the attributes from the file are in the data array + self.assertEqual(res.attrs['standard_name'], 'toa_brightness_temperature') + + +def suite(): + """Create the test suite for test_scene.""" + loader = unittest.TestLoader() + mysuite = unittest.TestSuite() + mysuite.addTest(loader.loadTestsFromTestCase(TestAMIL1bNetCDF)) + mysuite.addTest(loader.loadTestsFromTestCase(TestAMIL1bNetCDFIRCal)) + return mysuite + + +if __name__ == '__main__': + unittest.main()