#!/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 <http://www.gnu.org/licenses/>.
"""Interface to VIRR (Visible and Infra-Red Radiometer) level 1b format.
The file format is HDF5. Important attributes:
- Latitude
- Longitude
- SolarZenith
- EV_Emissive
- EV_RefSB
- Emissive_Radiance_Offsets
- Emissive_Radiance_Scales
- RefSB_Cal_Coefficients
- RefSB_Effective_Wavelength
- Emmisive_Centroid_Wave_Number
Supported satellites:
- FY-3B and FY-3C.
For more information:
- https://www.wmo-sat.info/oscar/instruments/view/607.
"""
import logging
from datetime import datetime
import dask.array as da
import numpy as np
from pyspectral.blackbody import blackbody_wn_rad2temp as rad2temp
from satpy.readers.hdf5_utils import HDF5FileHandler
LOG = logging.getLogger(__name__)
# PROVIDED BY NIGEL ATKINSON - 2013
# FY3B_REF_COEFFS = [
# 0.12640, -1.43200, #channel1#
# 0.13530, -1.62360, #channel2#
# 0.09193, -2.48207, #channel6#
# 0.07480, -0.90980, #channel7#
# 0.07590, -0.91080, #channel8#
# 0.07460, -0.89520, #channel9#
# 0.06300, -0.76280] #channel10#
# CMA - 2015 - http://www.nsmc.org.cn/en/NSMC/Contents/100089.html
FY3B_REF_COEFFS = [
0.1264, -1.4320,
0.1353, -1.6236,
0.0919, -2.4821,
0.0938, -1.1494,
0.0857, -1.0280,
0.0803, -0.9636,
0.0630, -0.7628]
[docs]class VIRR_L1B(HDF5FileHandler):
"""VIRR Level 1b reader."""
def __init__(self, filename, filename_info, filetype_info):
"""Open file and perform initial setup."""
super(VIRR_L1B, self).__init__(filename, filename_info, filetype_info)
LOG.debug('day/night flag for {0}: {1}'.format(filename, self['/attr/Day Or Night Flag']))
self.geolocation_prefix = filetype_info['geolocation_prefix']
self.platform_id = filename_info['platform_id']
self.l1b_prefix = 'Data/'
self.wave_number = 'Emissive_Centroid_Wave_Number'
# Else filename_info['platform_id'] == FY3C.
if filename_info['platform_id'] == 'FY3B':
self.l1b_prefix = ''
self.wave_number = 'Emmisive_Centroid_Wave_Number'
[docs] def get_dataset(self, dataset_id, ds_info):
"""Create DataArray from file content for `dataset_id`."""
file_key = self.geolocation_prefix + ds_info.get('file_key', dataset_id['name'])
if self.platform_id == 'FY3B':
file_key = file_key.replace('Data/', '')
data = self[file_key]
band_index = ds_info.get('band_index')
valid_range = data.attrs.pop('valid_range', None)
if isinstance(valid_range, np.ndarray):
valid_range = valid_range.tolist()
if band_index is not None:
data = data[band_index]
if valid_range:
data = data.where((data >= valid_range[0]) &
(data <= valid_range[1]))
if 'Emissive' in file_key:
self._calibrate_emissive(data, band_index)
elif 'RefSB' in file_key:
data = self._calibrate_reflective(data, band_index)
else:
slope = self._correct_slope(self[file_key + '/attr/Slope'])
intercept = self[file_key + '/attr/Intercept']
if valid_range:
data = data.where((data >= valid_range[0]) &
(data <= valid_range[1]))
data = data * slope + intercept
new_dims = {old: new for old, new in zip(data.dims, ('y', 'x'))}
data = data.rename(new_dims)
# use lowercase sensor name to be consistent with the rest of satpy
data.attrs.update({'platform_name': self['/attr/Satellite Name'],
'sensor': self['/attr/Sensor Identification Code'].lower()})
data.attrs.update(ds_info)
units = self.get(file_key + '/attr/units')
if units is not None and str(units).lower() != 'none':
data.attrs.update({'units': self.get(file_key + '/attr/units')})
elif data.attrs.get('calibration') == 'reflectance':
data.attrs.update({'units': '%'})
else:
data.attrs.update({'units': '1'})
return data
def _calibrate_reflective(self, data, band_index):
if self.platform_id == 'FY3B':
coeffs = da.from_array(FY3B_REF_COEFFS, chunks=-1)
else:
coeffs = self['/attr/RefSB_Cal_Coefficients']
slope = self._correct_slope(coeffs[0::2])
intercept = coeffs[1::2]
data = data * slope[band_index] + intercept[band_index]
return data
def _calibrate_emissive(self, data, band_index):
slope = self._correct_slope(self[self.l1b_prefix + 'Emissive_Radiance_Scales'].
data[:, band_index][:, np.newaxis])
intercept = self[self.l1b_prefix + 'Emissive_Radiance_Offsets'].data[:, band_index][:, np.newaxis]
# Converts cm^-1 (wavenumbers) and (mW/m^2)/(str/cm^-1) (radiance data)
# to SI units m^-1, mW*m^-3*str^-1.
wave_number = self['/attr/' + self.wave_number][band_index] * 100
bt_data = rad2temp(wave_number, (data.data * slope + intercept) * 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
def _correct_slope(self, slope):
# 0 slope is invalid. Note: slope can be a scalar or array.
return da.where(slope == 0, 1, slope)
@property
def start_time(self):
"""Get starting observation time."""
start_time = self['/attr/Observing Beginning Date'] + 'T' + self['/attr/Observing Beginning Time'] + 'Z'
return datetime.strptime(start_time, '%Y-%m-%dT%H:%M:%S.%fZ')
@property
def end_time(self):
"""Get ending observation time."""
end_time = self['/attr/Observing Ending Date'] + 'T' + self['/attr/Observing Ending Time'] + 'Z'
return datetime.strptime(end_time, '%Y-%m-%dT%H:%M:%S.%fZ')