#!/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/>.
"""Reader for the FY-3D MERSI-2 L1B file format.
The files for this reader are HDF5 and come in four varieties; band data
and geolocation data, both at 250m and 1000m resolution.
This reader was tested on FY-3D MERSI-2 data, but should work on future
platforms as well assuming no file format changes.
"""
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
N_TOT_IR_CHANS_LL = 6
[docs]
class MERSIL1B(HDF5FileHandler):
"""MERSI-2/MERSI-LL/MERSI-RM L1B file reader."""
[docs]
def _strptime(self, date_attr, time_attr):
"""Parse date/time strings."""
date = self[date_attr]
time = self[time_attr] # "18:27:39.720"
# cuts off microseconds because of unknown meaning
# is .720 == 720 microseconds or 720000 microseconds
return datetime.strptime(date + " " + time.split(".")[0], "%Y-%m-%d %H:%M:%S")
@property
def start_time(self):
"""Time for first observation."""
return self._strptime("/attr/Observing Beginning Date", "/attr/Observing Beginning Time")
@property
def end_time(self):
"""Time for final observation."""
return self._strptime("/attr/Observing Ending Date", "/attr/Observing Ending Time")
@property
def sensor_name(self):
"""Map sensor name to Satpy 'standard' sensor names."""
file_sensor = self["/attr/Sensor Identification Code"]
sensor = {
"MERSI": "mersi-2",
"MERSI LL": "mersi-ll",
"MERSI RM": "mersi-rm",
}.get(file_sensor, file_sensor)
return sensor
[docs]
def get_refl_mult(self):
"""Get reflectance multiplier."""
if self.sensor_name == "mersi-rm":
# MERSI-RM has reflectance in the range 0-1, so we need to convert
return 100.
else:
return 1.
[docs]
def _get_single_slope_intercept(self, slope, intercept, cal_index):
try:
# convert scalar arrays to scalar
return slope.item(), intercept.item()
except ValueError:
# numpy array but has more than one element
return slope[cal_index], intercept[cal_index]
[docs]
def _get_coefficients(self, cal_key, cal_index):
coeffs = self[cal_key][cal_index]
slope = coeffs.attrs.pop("Slope", None)
intercept = coeffs.attrs.pop("Intercept", None)
if slope is not None:
slope, intercept = self._get_single_slope_intercept(
slope, intercept, cal_index)
coeffs = coeffs * slope + intercept
return coeffs
[docs]
def get_dataset(self, dataset_id, ds_info):
"""Load data variable and metadata and calibrate if needed."""
file_key = ds_info.get("file_key", dataset_id["name"])
band_index = ds_info.get("band_index")
data = self[file_key]
if band_index is not None:
data = data[band_index]
if data.ndim >= 2:
data = data.rename({data.dims[-2]: "y", data.dims[-1]: "x"})
attrs = data.attrs.copy() # avoid contaminating other band loading
attrs.update(ds_info)
if "rows_per_scan" in self.filetype_info:
attrs.setdefault("rows_per_scan", self.filetype_info["rows_per_scan"])
data = self._mask_data(data, dataset_id, attrs)
slope = attrs.pop("Slope", None)
intercept = attrs.pop("Intercept", None)
if slope is not None and dataset_id.get("calibration") != "counts":
if band_index is not None and slope.size > 1:
slope = slope[band_index]
intercept = intercept[band_index]
data = data * slope + intercept
if dataset_id.get("calibration") == "reflectance":
coeffs = self._get_coefficients(ds_info["calibration_key"],
ds_info["calibration_index"])
data = coeffs[0] + coeffs[1] * data + coeffs[2] * data ** 2
data = data * self.get_refl_mult()
elif dataset_id.get("calibration") == "brightness_temperature":
calibration_index = ds_info["calibration_index"]
# Converts um^-1 (wavenumbers) and (mW/m^2)/(str/cm^-1) (radiance data)
# to SI units m^-1, mW*m^-3*str^-1.
wave_number = 1. / (dataset_id["wavelength"][1] / 1e6)
data = self._get_bt_dataset(data, calibration_index, wave_number)
data.attrs = attrs
# convert bytes to str
for key, val in attrs.items():
# python 3 only
if bytes is not str and isinstance(val, bytes):
data.attrs[key] = val.decode("utf8")
data.attrs.update({
"platform_name": self["/attr/Satellite Name"],
"sensor": self.sensor_name,
})
return data
[docs]
def _mask_data(self, data, dataset_id, attrs):
"""Mask the data using fill_value and valid_range attributes."""
fill_value = attrs.pop("FillValue", np.nan) # covered by valid_range
valid_range = attrs.pop("valid_range", None)
if dataset_id.get("calibration") == "counts":
# preserve integer type of counts if possible
attrs["_FillValue"] = fill_value
new_fill = data.dtype.type(fill_value)
else:
new_fill = np.nan
if valid_range is not None:
# Due to a bug in the valid_range upper limit in the 10.8(24) and 12.0(25)
# in the HDF data, this is hardcoded here.
if dataset_id["name"] in ["24", "25"] and valid_range[1] == 4095:
valid_range[1] = 25000
# typically bad_values == 65535, saturated == 65534
# dead detector == 65533
data = data.where((data >= valid_range[0]) &
(data <= valid_range[1]), new_fill)
return data
[docs]
def _get_bt_dataset(self, data, calibration_index, wave_number):
"""Get the dataset as brightness temperature.
Apparently we don't use these calibration factors for Rad -> BT::
coeffs = self._get_coefficients(ds_info['calibration_key'], calibration_index)
# coefficients are per-scan, we need to repeat the values for a
# clean alignment
coeffs = np.repeat(coeffs, data.shape[0] // coeffs.shape[1], axis=1)
coeffs = coeffs.rename({
coeffs.dims[0]: 'coefficients', coeffs.dims[1]: 'y'
}) # match data dims
data = coeffs[0] + coeffs[1] * data + coeffs[2] * data**2 + coeffs[3] * data**3
"""
# pass the dask array
bt_data = rad2temp(wave_number, data.data * 1e-5) # brightness temperature
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
# Some BT bands seem to have 0 in the first 10 columns
# and it is an invalid measurement, so let's mask
data = data.where(data != 0)
# additional corrections from the file
if self.sensor_name == "mersi-2":
corr_coeff_a = float(self["/attr/TBB_Trans_Coefficient_A"][calibration_index])
corr_coeff_b = float(self["/attr/TBB_Trans_Coefficient_B"][calibration_index])
elif self.sensor_name == "mersi-ll":
# MERSI-LL stores these coefficients differently
try:
coeffs = self["/attr/TBB_Trans_Coefficient"]
corr_coeff_a = coeffs[calibration_index]
corr_coeff_b = coeffs[calibration_index + N_TOT_IR_CHANS_LL]
except KeyError:
return data
else:
# MERSI-RM has no correction coefficients
corr_coeff_a = 0
if corr_coeff_a != 0:
data = (data - corr_coeff_b) / corr_coeff_a
# some bands have 0 counts for the first N columns and
# seem to be invalid data points
data = data.where(data != 0)
return data