From dca2e1bf36a8ef1662b75db130ea0dcdaa6dff0f Mon Sep 17 00:00:00 2001 From: David Hoese Date: Sun, 22 Sep 2019 10:46:11 -0500 Subject: [PATCH 01/19] Add initial GK-2A AMI L1B Reader --- satpy/etc/readers/ami_l1b.yaml | 293 +++++++++++++++++++++++++++++++++ satpy/readers/ami_l1b.py | 110 +++++++++++++ 2 files changed, 403 insertions(+) create mode 100644 satpy/etc/readers/ami_l1b.yaml create mode 100644 satpy/readers/ami_l1b.py diff --git a/satpy/etc/readers/ami_l1b.yaml b/satpy/etc/readers/ami_l1b.yaml new file mode 100644 index 0000000000..47cb47515e --- /dev/null +++ b/satpy/etc/readers/ami_l1b.yaml @@ -0,0 +1,293 @@ +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: + 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.576, 0.590, 0.614] + resolution: 1000 + calibration: + 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, .639, 0.679] + resolution: 500 + calibration: + 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: + 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: + 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: + 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: + 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: + 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: + 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: + 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: + 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: + 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: + 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: + 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: + 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: + 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..c3ad48775c --- /dev/null +++ b/satpy/readers/ami_l1b.py @@ -0,0 +1,110 @@ +#!/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 + +from pyresample import geometry +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): + """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_dims({'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' + + @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'] + 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'] + area_extent = ( + (0 - coff - 0.5) * cfac, + (0 - loff - 0.5) * lfac, + (cols - coff + 0.5) * cfac, + (rows - loff + 0.5) * 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_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] + return data From db0242c63244584bb93685e10a76488f0963bdc2 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Sun, 22 Sep 2019 10:50:08 -0500 Subject: [PATCH 02/19] Add AMI reader to docs reader table --- doc/source/index.rst | 3 +++ 1 file changed, 3 insertions(+) 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 From 2c8f6b315d623553c6add854bc662dcaab093c4c Mon Sep 17 00:00:00 2001 From: David Hoese Date: Sun, 22 Sep 2019 19:45:46 -0500 Subject: [PATCH 03/19] Add AMI calibration steps --- satpy/etc/composites/ami.yaml | 1 + satpy/etc/readers/ami_l1b.yaml | 2 +- satpy/readers/ami_l1b.py | 105 +++++++++++++++++++++++++++++---- 3 files changed, 96 insertions(+), 12 deletions(-) create mode 100644 satpy/etc/composites/ami.yaml diff --git a/satpy/etc/composites/ami.yaml b/satpy/etc/composites/ami.yaml new file mode 100644 index 0000000000..7f212ad759 --- /dev/null +++ b/satpy/etc/composites/ami.yaml @@ -0,0 +1 @@ +sensor_name: visir/ami diff --git a/satpy/etc/readers/ami_l1b.yaml b/satpy/etc/readers/ami_l1b.yaml index 47cb47515e..ef338c8a8d 100644 --- a/satpy/etc/readers/ami_l1b.yaml +++ b/satpy/etc/readers/ami_l1b.yaml @@ -97,7 +97,7 @@ datasets: C03: name: VI006 - wavelength: [0.599, .639, 0.679] + wavelength: [0.599, 0.639, 0.679] resolution: 500 calibration: radiance: diff --git a/satpy/readers/ami_l1b.py b/satpy/readers/ami_l1b.py index c3ad48775c..dd2c721b84 100644 --- a/satpy/readers/ami_l1b.py +++ b/satpy/readers/ami_l1b.py @@ -22,8 +22,10 @@ import numpy as np import xarray as xr +import dask.array as da 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 @@ -34,11 +36,36 @@ 'GK-2B': 'GEO-KOMPSAT-2B', } +# Copied from 20190415_GK-2A_AMI_Conversion_Table_v3.0.xlsx +# Sheet: coeff.& equation_WN +# Visible channels +# channel_name -> (DN2Rad_Gain, DN2Rad_Offset, Rad. to Albedo) +# IR channels +# channel_name -> (DN2Rad_Gain, DN2Rad_Offset, c0, c1, c2) +CALIBRATION_COEFFS = { + "VI004": (0.363545805215835, -7.27090454101562, 0.001558245), + "VI005": (0.343625485897064, -6.87249755859375, 0.0016595767), + "VI006": (0.154856294393539, -6.19424438476562, 0.001924484), + "VI008": (0.0457241721451282, -3.65792846679687, 0.0032723873), + "NR013": (0.0346878096461296, -1.38751220703125, 0.0087081313), + "NR016": (0.0498007982969284, -0.996017456054687, 0.0129512876), + "SW038": (-0.00108296517282724, 17.699987411499, -0.447843939824124, 1.00065568090389, -0.0000000633824089912448), + "WV063": (-0.0108914673328399, 44.1777038574218, -1.76279494011147, 1.00414910562278, -0.000000983310914319385), + "WV069": (-0.00818779878318309, 66.7480773925781, -0.334311414359106, 1.00097359874468, -0.000000494603070252304), + "WV073": (-0.0096982717514038, 79.0608520507812, -0.0613124859696595, 1.00019008722941, -0.000000105863656750499), + "IR087": (-0.0144806550815701, 118.050903320312, -0.141418528203155, 1.00052232906885, -0.00000036287276076109), + "IR096": (-0.0178435463458299, 145.464874267578, -0.114017728158198, 1.00047380585402, -0.000000374931509928403), + "IR105": (-0.0198196955025196, 161.580139160156, -0.142866448475177, 1.00064069572049, -0.000000550443294960498), + "IR112": (-0.0216744858771562, 176.713439941406, -0.249111718496148, 1.00121166873756, -0.00000113167964011665), + "IR123": (-0.023379972204566, 190.649627685546, -0.458113885722738, 1.00245520975535, -0.00000253064314720476), + "IR133": (-0.0243037566542625, 198.224365234375, -0.0938521568527657, 1.00053982112966, -0.000000594913715312849), +} + class AMIL1bNetCDF(BaseFileHandler): """Base reader for AMI L1B NetCDF4 files.""" - def __init__(self, filename, filename_info, filetype_info): + def __init__(self, filename, filename_info, filetype_info, 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, @@ -50,6 +77,7 @@ def __init__(self, filename, filename_info, filetype_info): 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 @property def start_time(self): @@ -78,19 +106,22 @@ def get_area_def(self, dsid): coff = self.nc.attrs['coff'] lfac = self.nc.attrs['lfac'] loff = self.nc.attrs['loff'] + bit_shift = 2**16 area_extent = ( - (0 - coff - 0.5) * cfac, - (0 - loff - 0.5) * lfac, - (cols - coff + 0.5) * cfac, - (rows - loff + 0.5) * lfac, + 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'} + 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()), @@ -107,4 +138,56 @@ 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 + # noticing better results for some bands when using: + # data &= 2**14 - 1 + # only take "no error" pixels as valid + data = data.where(qf == 0) + + coeffs = CALIBRATION_COEFFS.get(dataset_id.name) + if coeffs is None and dataset_id.calibration is not None: + raise ValueError("No coefficients configured for {}".format(dataset_id)) + if dataset_id.calibration in ('radiance', 'reflectance', 'brightness_temperature'): + gain = coeffs[0] + offset = coeffs[1] + data = gain * data + offset + if dataset_id.calibration == 'reflectance': + # depends on the radiance calibration above + rad_to_alb = coeffs[2] + if ds_info.get('units') == '%': + rad_to_alb *= 100 + data = data * rad_to_alb + # print(da.compute(np.nanmin(data.data), np.nanmax(data.data))) + elif dataset_id.calibration == 'brightness_temperature': + # 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 + + for attr_name in ('standard_name', 'units'): + attrs[attr_name] = ds_info[attr_name] + attrs.update(dataset_id.to_dict()) + data.attrs = attrs return data From d3bf76f6e7ff05f95f0e0da5af82d1a9996711f8 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Wed, 25 Sep 2019 13:46:16 -0500 Subject: [PATCH 04/19] Fix typo in VI005 AMI wavelength --- satpy/etc/readers/ami_l1b.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/satpy/etc/readers/ami_l1b.yaml b/satpy/etc/readers/ami_l1b.yaml index ef338c8a8d..d9858bd979 100644 --- a/satpy/etc/readers/ami_l1b.yaml +++ b/satpy/etc/readers/ami_l1b.yaml @@ -83,7 +83,7 @@ datasets: C02: name: VI005 - wavelength: [0.576, 0.590, 0.614] + wavelength: [0.495, 0.509, 0.523] resolution: 1000 calibration: radiance: From 9dca4072c515f4aef1522f59d147f64685164cd4 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Wed, 25 Sep 2019 13:47:17 -0500 Subject: [PATCH 05/19] Make AHI green corrector compositor take fraction parameters --- satpy/composites/ahi.py | 8 ++++++-- satpy/etc/composites/ami.yaml | 18 ++++++++++++++++++ 2 files changed, 24 insertions(+), 2 deletions(-) diff --git a/satpy/composites/ahi.py b/satpy/composites/ahi.py index ab069609b2..2f3d031401 100644 --- a/satpy/composites/ahi.py +++ b/satpy/composites/ahi.py @@ -29,6 +29,10 @@ class GreenCorrector(GenericCompositor): """Corrector of the AHI green band to compensate for the deficit of chlorophyl signal. """ + def __init__(self, *args, **kwargs): + # 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.""" @@ -36,7 +40,7 @@ def __call__(self, projectables, optional_datasets=None, **attrs): 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() + new_green.data = new_green.data.persist() return super(GreenCorrector, self).__call__((new_green,), **attrs) diff --git a/satpy/etc/composites/ami.yaml b/satpy/etc/composites/ami.yaml index 7f212ad759..60ed80bed4 100644 --- a/satpy/etc/composites/ami.yaml +++ b/satpy/etc/composites/ami.yaml @@ -1 +1,19 @@ sensor_name: visir/ami + +composites: + green: + compositor: !!python/name:satpy.composites.ahi.GreenCorrector + prerequisites: + - VI005 + - VI008 + standard_name: toa_bidirectional_reflectance + fractions: [0.85, 0.15] + + true_color_raw: + compositor: !!python/name:satpy.composites.GenericCompositor + prerequisites: + - VI006 +# - VI005 + - green + - VI004 + standard_name: true_color From 6c4aa7643044b7d94a159bfbece36a6ee4eef993 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Thu, 26 Sep 2019 12:02:23 -0500 Subject: [PATCH 06/19] Use channel_name attribute instead of dataset name in AMI reader --- satpy/composites/ahi.py | 8 +++----- satpy/readers/ami_l1b.py | 4 ++-- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/satpy/composites/ahi.py b/satpy/composites/ahi.py index 2f3d031401..13e34925ab 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,9 +25,8 @@ 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): # XXX: Should this be 0.93 and 0.07 self.fractions = kwargs.pop('fractions', [0.85, 0.15]) diff --git a/satpy/readers/ami_l1b.py b/satpy/readers/ami_l1b.py index dd2c721b84..be433d3f51 100644 --- a/satpy/readers/ami_l1b.py +++ b/satpy/readers/ami_l1b.py @@ -158,7 +158,8 @@ def get_dataset(self, dataset_id, ds_info): # only take "no error" pixels as valid data = data.where(qf == 0) - coeffs = CALIBRATION_COEFFS.get(dataset_id.name) + channel_name = attrs.get('channel_name', dataset_id.name) + coeffs = CALIBRATION_COEFFS.get(channel_name) if coeffs is None and dataset_id.calibration is not None: raise ValueError("No coefficients configured for {}".format(dataset_id)) if dataset_id.calibration in ('radiance', 'reflectance', 'brightness_temperature'): @@ -171,7 +172,6 @@ def get_dataset(self, dataset_id, ds_info): if ds_info.get('units') == '%': rad_to_alb *= 100 data = data * rad_to_alb - # print(da.compute(np.nanmin(data.data), np.nanmax(data.data))) elif dataset_id.calibration == 'brightness_temperature': # depends on the radiance calibration above # Convert um to m^-1 (SI units for pyspectral) From da31be8c99fb535858b39cc877f4e1d880aa7990 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Thu, 3 Oct 2019 13:32:35 -0500 Subject: [PATCH 07/19] Add AMI orbital_parameters and true_color recipes --- satpy/composites/ahi.py | 3 +-- satpy/etc/composites/ami.yaml | 36 +++++++++++++++++++++++++++-------- satpy/readers/ami_l1b.py | 31 +++++++++++++++++++++++++++++- 3 files changed, 59 insertions(+), 11 deletions(-) diff --git a/satpy/composites/ahi.py b/satpy/composites/ahi.py index 13e34925ab..c170eb0543 100644 --- a/satpy/composites/ahi.py +++ b/satpy/composites/ahi.py @@ -28,17 +28,16 @@ class GreenCorrector(GenericCompositor): """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') new_green = green * self.fractions[0] + nir * self.fractions[1] new_green.attrs = green.attrs.copy() - new_green.data = new_green.data.persist() return super(GreenCorrector, self).__call__((new_green,), **attrs) diff --git a/satpy/etc/composites/ami.yaml b/satpy/etc/composites/ami.yaml index 60ed80bed4..b8f88fb250 100644 --- a/satpy/etc/composites/ami.yaml +++ b/satpy/etc/composites/ami.yaml @@ -1,19 +1,39 @@ 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: - - VI005 - - VI008 + - 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.GenericCompositor + 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: - - VI006 -# - VI005 - - green - - VI004 + - name: VI006 + modifiers: [sunz_corrected, rayleigh_corrected] + - name: green + - name: VI004 + modifiers: [sunz_corrected, rayleigh_corrected] standard_name: true_color diff --git a/satpy/readers/ami_l1b.py b/satpy/readers/ami_l1b.py index be433d3f51..a45b1e3a91 100644 --- a/satpy/readers/ami_l1b.py +++ b/satpy/readers/ami_l1b.py @@ -23,6 +23,7 @@ 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 @@ -95,7 +96,7 @@ 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'] + 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'] @@ -134,6 +135,31 @@ def get_area_def(self, dsid): 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_nominal_longitude': sc_position[0], + 'satellite_nominal_latitude': sc_position[1], + 'satellite_nominal_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) @@ -189,5 +215,8 @@ def get_dataset(self, dataset_id, ds_info): 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 From b1dc6d3e67bb54c34e1195b29596807702f68de0 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Sun, 6 Oct 2019 13:57:16 -0500 Subject: [PATCH 08/19] Add AMI reader tests and counts calibration --- satpy/etc/readers/ami_l1b.yaml | 48 ++++ satpy/readers/ami_l1b.py | 14 +- satpy/tests/reader_tests/__init__.py | 4 +- satpy/tests/reader_tests/test_ami_l1b.py | 292 +++++++++++++++++++++++ 4 files changed, 350 insertions(+), 8 deletions(-) create mode 100644 satpy/tests/reader_tests/test_ami_l1b.py diff --git a/satpy/etc/readers/ami_l1b.yaml b/satpy/etc/readers/ami_l1b.yaml index d9858bd979..b0311e7b35 100644 --- a/satpy/etc/readers/ami_l1b.yaml +++ b/satpy/etc/readers/ami_l1b.yaml @@ -72,6 +72,9 @@ datasets: 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 @@ -86,6 +89,9 @@ datasets: 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 @@ -100,6 +106,9 @@ datasets: 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 @@ -114,6 +123,9 @@ datasets: 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 @@ -128,6 +140,9 @@ datasets: 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 @@ -142,6 +157,9 @@ datasets: 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 @@ -156,6 +174,9 @@ datasets: 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 @@ -170,6 +191,9 @@ datasets: 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 @@ -184,6 +208,9 @@ datasets: 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 @@ -198,6 +225,9 @@ datasets: 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 @@ -212,6 +242,9 @@ datasets: 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 @@ -226,6 +259,9 @@ datasets: 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 @@ -240,6 +276,9 @@ datasets: 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 @@ -254,6 +293,9 @@ datasets: 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 @@ -268,6 +310,9 @@ datasets: 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 @@ -282,6 +327,9 @@ datasets: 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 diff --git a/satpy/readers/ami_l1b.py b/satpy/readers/ami_l1b.py index a45b1e3a91..abc407d575 100644 --- a/satpy/readers/ami_l1b.py +++ b/satpy/readers/ami_l1b.py @@ -73,7 +73,7 @@ def __init__(self, filename, filename_info, filetype_info, allow_conditional_pix decode_cf=True, mask_and_scale=False, chunks={'dim_image_x': CHUNK_SIZE, 'dim_image_y': CHUNK_SIZE}) - self.nc = self.nc.rename_dims({'dim_image_x': 'x', 'dim_image_y': 'y'}) + 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) @@ -154,9 +154,9 @@ def get_orbital_parameters(self): 'projection_longitude': float(lon_0), 'projection_latitude': 0.0, 'projection_altitude': h, - 'satellite_nominal_longitude': sc_position[0], - 'satellite_nominal_latitude': sc_position[1], - 'satellite_nominal_altitude': sc_position[2] / 1000.0, # km + 'satellite_actual_longitude': sc_position[0], + 'satellite_actual_latitude': sc_position[1], + 'satellite_actual_altitude': sc_position[2] / 1000.0, # km } return orbital_parameters @@ -179,8 +179,6 @@ def get_dataset(self, dataset_id, ds_info): # mask DQF bits bits = attrs['number_of_valid_bits_per_pixel'] data &= 2**bits - 1 - # noticing better results for some bands when using: - # data &= 2**14 - 1 # only take "no error" pixels as valid data = data.where(qf == 0) @@ -201,7 +199,7 @@ def get_dataset(self, dataset_id, ds_info): elif dataset_id.calibration == 'brightness_temperature': # depends on the radiance calibration above # Convert um to m^-1 (SI units for pyspectral) - wn = 1 / dataset_id.wavelength[1] / 1e6 + 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) @@ -211,6 +209,8 @@ def get_dataset(self, dataset_id, ds_info): else: # new versions of pyspectral can do dask arrays data.data = bt_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] 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..5f1fdadf92 --- /dev/null +++ b/satpy/tests/reader_tests/test_ami_l1b.py @@ -0,0 +1,292 @@ +#!/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", + } + ) + + 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 + res = self.reader.get_dataset( + DatasetID(name='IR087', wavelength=[8.415, 8.59, 8.765], + calibration='brightness_temperature'), + { + 'file_key': 'image_pixel_values', + 'wavelength': [8.415, 8.59, 8.765], + 'standard_name': 'toa_brightness_temperature', + 'units': 'K', + }) + + expected = np.array([[238.34385135, 238.31443527, 238.28500087, 238.25554813, 238.22607701], + [238.1965875, 238.16707956, 238.13755317, 238.10800829, 238.07844489]]) + self.assertTrue(np.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') + + +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() From a2e4de0902be33a8aba8630b74552d542aab4dbb Mon Sep 17 00:00:00 2001 From: David Hoese Date: Tue, 8 Oct 2019 08:40:31 -0500 Subject: [PATCH 09/19] Fix AMI area definition being flipped vertically --- satpy/readers/ami_l1b.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/satpy/readers/ami_l1b.py b/satpy/readers/ami_l1b.py index abc407d575..6cf22c83ba 100644 --- a/satpy/readers/ami_l1b.py +++ b/satpy/readers/ami_l1b.py @@ -110,9 +110,9 @@ def get_area_def(self, dsid): 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((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), + -h * np.deg2rad((rows - loff + 0.5) * bit_shift / lfac), ) proj_dict = { From 2b30720f3d2cdc6f83ea2600f4a780a25eef39b7 Mon Sep 17 00:00:00 2001 From: simonrp84 Date: Tue, 8 Oct 2019 14:51:42 +0100 Subject: [PATCH 10/19] Use internal calibration coefficients if available --- satpy/readers/ami_l1b.py | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/satpy/readers/ami_l1b.py b/satpy/readers/ami_l1b.py index abc407d575..c01895e780 100644 --- a/satpy/readers/ami_l1b.py +++ b/satpy/readers/ami_l1b.py @@ -110,9 +110,9 @@ def get_area_def(self, dsid): 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((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), + -h * np.deg2rad((rows - loff + 0.5) * bit_shift / lfac), ) proj_dict = { @@ -183,12 +183,19 @@ def get_dataset(self, dataset_id, ds_info): data = data.where(qf == 0) channel_name = attrs.get('channel_name', dataset_id.name) - coeffs = CALIBRATION_COEFFS.get(channel_name) - if coeffs is None and dataset_id.calibration is not None: - raise ValueError("No coefficients configured for {}".format(dataset_id)) + + # 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'): - gain = coeffs[0] - offset = coeffs[1] + coeffs = CALIBRATION_COEFFS.get(channel_name) + if coeffs is None and dataset_id.calibration is not None: + raise ValueError("No coefficients configured for {}".format(dataset_id)) + if (gain < -900 or offset < -900): + print("WARNING: No calib") + gain = coeffs[0] + offset = coeffs[1] data = gain * data + offset if dataset_id.calibration == 'reflectance': # depends on the radiance calibration above From e714ef222b8a8e806bdde9b50351433ba6e4949f Mon Sep 17 00:00:00 2001 From: simonrp84 Date: Tue, 8 Oct 2019 16:45:05 +0100 Subject: [PATCH 11/19] Remove built-in slope/offset, use values in file. Add option to switch between pyspectral and AMI brightness temperature conversion --- satpy/readers/ami_l1b.py | 89 +++++++++++++++++++--------------------- 1 file changed, 43 insertions(+), 46 deletions(-) diff --git a/satpy/readers/ami_l1b.py b/satpy/readers/ami_l1b.py index c01895e780..04048c2d2a 100644 --- a/satpy/readers/ami_l1b.py +++ b/satpy/readers/ami_l1b.py @@ -37,36 +37,12 @@ 'GK-2B': 'GEO-KOMPSAT-2B', } -# Copied from 20190415_GK-2A_AMI_Conversion_Table_v3.0.xlsx -# Sheet: coeff.& equation_WN -# Visible channels -# channel_name -> (DN2Rad_Gain, DN2Rad_Offset, Rad. to Albedo) -# IR channels -# channel_name -> (DN2Rad_Gain, DN2Rad_Offset, c0, c1, c2) -CALIBRATION_COEFFS = { - "VI004": (0.363545805215835, -7.27090454101562, 0.001558245), - "VI005": (0.343625485897064, -6.87249755859375, 0.0016595767), - "VI006": (0.154856294393539, -6.19424438476562, 0.001924484), - "VI008": (0.0457241721451282, -3.65792846679687, 0.0032723873), - "NR013": (0.0346878096461296, -1.38751220703125, 0.0087081313), - "NR016": (0.0498007982969284, -0.996017456054687, 0.0129512876), - "SW038": (-0.00108296517282724, 17.699987411499, -0.447843939824124, 1.00065568090389, -0.0000000633824089912448), - "WV063": (-0.0108914673328399, 44.1777038574218, -1.76279494011147, 1.00414910562278, -0.000000983310914319385), - "WV069": (-0.00818779878318309, 66.7480773925781, -0.334311414359106, 1.00097359874468, -0.000000494603070252304), - "WV073": (-0.0096982717514038, 79.0608520507812, -0.0613124859696595, 1.00019008722941, -0.000000105863656750499), - "IR087": (-0.0144806550815701, 118.050903320312, -0.141418528203155, 1.00052232906885, -0.00000036287276076109), - "IR096": (-0.0178435463458299, 145.464874267578, -0.114017728158198, 1.00047380585402, -0.000000374931509928403), - "IR105": (-0.0198196955025196, 161.580139160156, -0.142866448475177, 1.00064069572049, -0.000000550443294960498), - "IR112": (-0.0216744858771562, 176.713439941406, -0.249111718496148, 1.00121166873756, -0.00000113167964011665), - "IR123": (-0.023379972204566, 190.649627685546, -0.458113885722738, 1.00245520975535, -0.00000253064314720476), - "IR133": (-0.0243037566542625, 198.224365234375, -0.0938521568527657, 1.00053982112966, -0.000000594913715312849), -} - class AMIL1bNetCDF(BaseFileHandler): """Base reader for AMI L1B NetCDF4 files.""" - def __init__(self, filename, filename_info, filetype_info, allow_conditional_pixels=False): + 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, @@ -79,6 +55,11 @@ def __init__(self, filename, filename_info, filetype_info, allow_conditional_pix 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): @@ -182,39 +163,55 @@ def get_dataset(self, dataset_id, ds_info): # only take "no error" pixels as valid data = data.where(qf == 0) - channel_name = attrs.get('channel_name', dataset_id.name) - # 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'): - coeffs = CALIBRATION_COEFFS.get(channel_name) - if coeffs is None and dataset_id.calibration is not None: - raise ValueError("No coefficients configured for {}".format(dataset_id)) - if (gain < -900 or offset < -900): - print("WARNING: No calib") - gain = coeffs[0] - offset = coeffs[1] data = gain * data + offset if dataset_id.calibration == 'reflectance': # depends on the radiance calibration above - rad_to_alb = coeffs[2] + 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': - # 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) + print(self.calib_mode) + 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: - # new versions of pyspectral can do dask arrays + # 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 elif dataset_id.calibration not in ('counts', 'radiance'): raise ValueError("Unknown calibration: '{}'".format(dataset_id.calibration)) From e6872b8da5c756c816fb6ac20a9983f60c096e48 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Tue, 8 Oct 2019 19:35:27 -0500 Subject: [PATCH 12/19] Fix AMI area extent test to match code --- satpy/tests/reader_tests/test_ami_l1b.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/satpy/tests/reader_tests/test_ami_l1b.py b/satpy/tests/reader_tests/test_ami_l1b.py index 5f1fdadf92..59d0e8c608 100644 --- a/satpy/tests/reader_tests/test_ami_l1b.py +++ b/satpy/tests/reader_tests/test_ami_l1b.py @@ -190,7 +190,7 @@ def test_get_area_def(self, adef): 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]) + [-5511523.904082, -5511523.904082, 5511022.902, 5511022.902]) def test_get_dataset_vis(self): """Test get visible calibrated data.""" From 46cd7d5afe3f09d23a88471c0d65b1127e3942fe Mon Sep 17 00:00:00 2001 From: simonrp84 Date: Tue, 8 Oct 2019 14:51:42 +0100 Subject: [PATCH 13/19] Use internal calibration coefficients if available --- satpy/readers/ami_l1b.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/satpy/readers/ami_l1b.py b/satpy/readers/ami_l1b.py index 6cf22c83ba..c01895e780 100644 --- a/satpy/readers/ami_l1b.py +++ b/satpy/readers/ami_l1b.py @@ -183,12 +183,19 @@ def get_dataset(self, dataset_id, ds_info): data = data.where(qf == 0) channel_name = attrs.get('channel_name', dataset_id.name) - coeffs = CALIBRATION_COEFFS.get(channel_name) - if coeffs is None and dataset_id.calibration is not None: - raise ValueError("No coefficients configured for {}".format(dataset_id)) + + # 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'): - gain = coeffs[0] - offset = coeffs[1] + coeffs = CALIBRATION_COEFFS.get(channel_name) + if coeffs is None and dataset_id.calibration is not None: + raise ValueError("No coefficients configured for {}".format(dataset_id)) + if (gain < -900 or offset < -900): + print("WARNING: No calib") + gain = coeffs[0] + offset = coeffs[1] data = gain * data + offset if dataset_id.calibration == 'reflectance': # depends on the radiance calibration above From 35de05887fafc37e68b430c546acb78984cdd797 Mon Sep 17 00:00:00 2001 From: simonrp84 Date: Tue, 8 Oct 2019 16:45:05 +0100 Subject: [PATCH 14/19] Remove built-in slope/offset, use values in file. Add option to switch between pyspectral and AMI brightness temperature conversion --- satpy/readers/ami_l1b.py | 89 +++++++++++++++++++--------------------- 1 file changed, 43 insertions(+), 46 deletions(-) diff --git a/satpy/readers/ami_l1b.py b/satpy/readers/ami_l1b.py index c01895e780..04048c2d2a 100644 --- a/satpy/readers/ami_l1b.py +++ b/satpy/readers/ami_l1b.py @@ -37,36 +37,12 @@ 'GK-2B': 'GEO-KOMPSAT-2B', } -# Copied from 20190415_GK-2A_AMI_Conversion_Table_v3.0.xlsx -# Sheet: coeff.& equation_WN -# Visible channels -# channel_name -> (DN2Rad_Gain, DN2Rad_Offset, Rad. to Albedo) -# IR channels -# channel_name -> (DN2Rad_Gain, DN2Rad_Offset, c0, c1, c2) -CALIBRATION_COEFFS = { - "VI004": (0.363545805215835, -7.27090454101562, 0.001558245), - "VI005": (0.343625485897064, -6.87249755859375, 0.0016595767), - "VI006": (0.154856294393539, -6.19424438476562, 0.001924484), - "VI008": (0.0457241721451282, -3.65792846679687, 0.0032723873), - "NR013": (0.0346878096461296, -1.38751220703125, 0.0087081313), - "NR016": (0.0498007982969284, -0.996017456054687, 0.0129512876), - "SW038": (-0.00108296517282724, 17.699987411499, -0.447843939824124, 1.00065568090389, -0.0000000633824089912448), - "WV063": (-0.0108914673328399, 44.1777038574218, -1.76279494011147, 1.00414910562278, -0.000000983310914319385), - "WV069": (-0.00818779878318309, 66.7480773925781, -0.334311414359106, 1.00097359874468, -0.000000494603070252304), - "WV073": (-0.0096982717514038, 79.0608520507812, -0.0613124859696595, 1.00019008722941, -0.000000105863656750499), - "IR087": (-0.0144806550815701, 118.050903320312, -0.141418528203155, 1.00052232906885, -0.00000036287276076109), - "IR096": (-0.0178435463458299, 145.464874267578, -0.114017728158198, 1.00047380585402, -0.000000374931509928403), - "IR105": (-0.0198196955025196, 161.580139160156, -0.142866448475177, 1.00064069572049, -0.000000550443294960498), - "IR112": (-0.0216744858771562, 176.713439941406, -0.249111718496148, 1.00121166873756, -0.00000113167964011665), - "IR123": (-0.023379972204566, 190.649627685546, -0.458113885722738, 1.00245520975535, -0.00000253064314720476), - "IR133": (-0.0243037566542625, 198.224365234375, -0.0938521568527657, 1.00053982112966, -0.000000594913715312849), -} - class AMIL1bNetCDF(BaseFileHandler): """Base reader for AMI L1B NetCDF4 files.""" - def __init__(self, filename, filename_info, filetype_info, allow_conditional_pixels=False): + 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, @@ -79,6 +55,11 @@ def __init__(self, filename, filename_info, filetype_info, allow_conditional_pix 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): @@ -182,39 +163,55 @@ def get_dataset(self, dataset_id, ds_info): # only take "no error" pixels as valid data = data.where(qf == 0) - channel_name = attrs.get('channel_name', dataset_id.name) - # 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'): - coeffs = CALIBRATION_COEFFS.get(channel_name) - if coeffs is None and dataset_id.calibration is not None: - raise ValueError("No coefficients configured for {}".format(dataset_id)) - if (gain < -900 or offset < -900): - print("WARNING: No calib") - gain = coeffs[0] - offset = coeffs[1] data = gain * data + offset if dataset_id.calibration == 'reflectance': # depends on the radiance calibration above - rad_to_alb = coeffs[2] + 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': - # 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) + print(self.calib_mode) + 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: - # new versions of pyspectral can do dask arrays + # 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 elif dataset_id.calibration not in ('counts', 'radiance'): raise ValueError("Unknown calibration: '{}'".format(dataset_id.calibration)) From 552b2e7402d64630d823829903e4261e729ecb4b Mon Sep 17 00:00:00 2001 From: simonrp84 Date: Fri, 11 Oct 2019 17:01:07 +0100 Subject: [PATCH 15/19] Add tests --- satpy/tests/reader_tests/test_ami_l1b.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/satpy/tests/reader_tests/test_ami_l1b.py b/satpy/tests/reader_tests/test_ami_l1b.py index 59d0e8c608..dbb81ebf13 100644 --- a/satpy/tests/reader_tests/test_ami_l1b.py +++ b/satpy/tests/reader_tests/test_ami_l1b.py @@ -115,6 +115,15 @@ def setUp(self, xr_, counts=None): "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, } ) From 6088e5f76c0dacd05e913eca185093add2536a26 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Fri, 11 Oct 2019 12:11:15 -0500 Subject: [PATCH 16/19] Move AMI IR calibration to separate method and add test for file calib --- satpy/readers/ami_l1b.py | 78 +++++++++++++----------- satpy/tests/reader_tests/test_ami_l1b.py | 29 +++++---- 2 files changed, 59 insertions(+), 48 deletions(-) diff --git a/satpy/readers/ami_l1b.py b/satpy/readers/ami_l1b.py index 04048c2d2a..7867d679d2 100644 --- a/satpy/readers/ami_l1b.py +++ b/satpy/readers/ami_l1b.py @@ -176,43 +176,7 @@ def get_dataset(self, dataset_id, ds_info): rad_to_alb *= 100 data = data * rad_to_alb elif dataset_id.calibration == 'brightness_temperature': - print(self.calib_mode) - 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 + data = self._calibrate_ir(dataset_id, data) elif dataset_id.calibration not in ('counts', 'radiance'): raise ValueError("Unknown calibration: '{}'".format(dataset_id.calibration)) @@ -224,3 +188,43 @@ def get_dataset(self, dataset_id, ds_info): 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/test_ami_l1b.py b/satpy/tests/reader_tests/test_ami_l1b.py index dbb81ebf13..df0e878850 100644 --- a/satpy/tests/reader_tests/test_ami_l1b.py +++ b/satpy/tests/reader_tests/test_ami_l1b.py @@ -271,19 +271,26 @@ def setUp(self): def test_ir_calibrate(self): """Test IR calibration.""" from satpy import DatasetID - res = self.reader.get_dataset( - DatasetID(name='IR087', wavelength=[8.415, 8.59, 8.765], - calibration='brightness_temperature'), - { - 'file_key': 'image_pixel_values', - 'wavelength': [8.415, 8.59, 8.765], - 'standard_name': 'toa_brightness_temperature', - 'units': 'K', - }) - + 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', + } + res = self.reader.get_dataset(ds_id, ds_info) expected = np.array([[238.34385135, 238.31443527, 238.28500087, 238.25554813, 238.22607701], [238.1965875, 238.16707956, 238.13755317, 238.10800829, 238.07844489]]) - self.assertTrue(np.allclose(res.data.compute(), expected, equal_nan=True)) + 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' + res = self.reader.get_dataset(ds_id, ds_info) + # 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') From 610436661112f3de45ef0d94c48f22d416be92bd Mon Sep 17 00:00:00 2001 From: David Hoese Date: Fri, 11 Oct 2019 15:30:51 -0500 Subject: [PATCH 17/19] Copy AHI composites to AMI composite configs --- satpy/composites/__init__.py | 11 +++-- satpy/etc/composites/ami.yaml | 86 +++++++++++++++++++++++++++++++++++ 2 files changed, 93 insertions(+), 4 deletions(-) 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/etc/composites/ami.yaml b/satpy/etc/composites/ami.yaml index b8f88fb250..e90c90beb1 100644 --- a/satpy/etc/composites/ami.yaml +++ b/satpy/etc/composites/ami.yaml @@ -10,6 +10,7 @@ composites: modifiers: [sunz_corrected] standard_name: toa_bidirectional_reflectance fractions: [0.85, 0.15] + green: compositor: !!python/name:satpy.composites.ahi.GreenCorrector prerequisites: @@ -19,6 +20,7 @@ composites: modifiers: [sunz_corrected] standard_name: toa_bidirectional_reflectance fractions: [0.85, 0.15] + true_color_raw: compositor: !!python/name:satpy.composites.SelfSharpenedRGB prerequisites: @@ -28,6 +30,7 @@ composites: - name: VI004 modifiers: [sunz_corrected] standard_name: true_color + true_color: compositor: !!python/name:satpy.composites.SelfSharpenedRGB prerequisites: @@ -37,3 +40,86 @@ composites: - 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 From 7aeb0d5559f72a7378069eb61066f5cf163f0795 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Sun, 13 Oct 2019 09:05:49 -0500 Subject: [PATCH 18/19] Add AMI airmass and ash recipes --- satpy/etc/composites/ami.yaml | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/satpy/etc/composites/ami.yaml b/satpy/etc/composites/ami.yaml index e90c90beb1..ae670df403 100644 --- a/satpy/etc/composites/ami.yaml +++ b/satpy/etc/composites/ami.yaml @@ -123,3 +123,31 @@ composites: 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 From 4cfec8220f594c9d3c02ba45971fd4243c4e9c2c Mon Sep 17 00:00:00 2001 From: David Hoese Date: Mon, 21 Oct 2019 12:08:12 -0500 Subject: [PATCH 19/19] Add assertion to make sure AMI calibration comes from file --- satpy/tests/reader_tests/test_ami_l1b.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/satpy/tests/reader_tests/test_ami_l1b.py b/satpy/tests/reader_tests/test_ami_l1b.py index df0e878850..ff2667fb90 100644 --- a/satpy/tests/reader_tests/test_ami_l1b.py +++ b/satpy/tests/reader_tests/test_ami_l1b.py @@ -271,6 +271,7 @@ def setUp(self): 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 = { @@ -279,7 +280,9 @@ def test_ir_calibrate(self): 'standard_name': 'toa_brightness_temperature', 'units': 'K', } - res = self.reader.get_dataset(ds_id, ds_info) + 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) @@ -288,7 +291,9 @@ def test_ir_calibrate(self): # test builtin coefficients self.reader.calib_mode = 'FILE' - res = self.reader.get_dataset(ds_id, ds_info) + 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