#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright (c) 2017 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/>.
"""HRIT format reader.
References:
ELECTRO-L GROUND SEGMENT MSU-GS INSTRUMENT,
LRIT/HRIT Mission Specific Implementation, February 2012
"""
import logging
from datetime import datetime
import numpy as np
import xarray as xr
from satpy.readers._geos_area import get_area_definition, get_area_extent
from satpy.readers.hrit_base import (
HRITFileHandler,
ancillary_text,
annotation_header,
base_hdr_map,
image_data_function,
time_cds_short,
)
logger = logging.getLogger("hrit_electrol")
# goms implementation:
key_header = np.dtype([("key_number", "u1"),
("seed", ">f8")])
segment_identification = np.dtype([("GP_SC_ID", ">i2"),
("spectral_channel_id", ">i1"),
("segment_sequence_number", ">u2"),
("planned_start_segment_number", ">u2"),
("planned_end_segment_number", ">u2"),
("data_field_representation", ">i1")])
image_segment_line_quality = np.dtype([("line_number_in_grid", ">i4"),
("line_mean_acquisition",
[("days", ">u2"),
("milliseconds", ">u4")]),
("line_validity", "u1"),
("line_radiometric_quality", "u1"),
("line_geometric_quality", "u1")])
goms_variable_length_headers = {
image_segment_line_quality: "image_segment_line_quality"}
goms_text_headers = {image_data_function: "image_data_function",
annotation_header: "annotation_header",
ancillary_text: "ancillary_text"}
goms_hdr_map = base_hdr_map.copy()
goms_hdr_map.update({7: key_header,
128: segment_identification,
129: image_segment_line_quality
})
orbit_coef = np.dtype([("StartTime", time_cds_short),
("EndTime", time_cds_short),
("X", ">f8", (8, )),
("Y", ">f8", (8, )),
("Z", ">f8", (8, )),
("VX", ">f8", (8, )),
("VY", ">f8", (8, )),
("VZ", ">f8", (8, ))])
attitude_coef = np.dtype([("StartTime", time_cds_short),
("EndTime", time_cds_short),
("XofSpinAxis", ">f8", (8, )),
("YofSpinAxis", ">f8", (8, )),
("ZofSpinAxis", ">f8", (8, ))])
cuc_time = np.dtype([("coarse", "u1", (4, )),
("fine", "u1", (3, ))])
time_cds_expanded = np.dtype([("days", ">u2"),
("milliseconds", ">u4"),
("microseconds", ">u2"),
("nanoseconds", ">u2")])
satellite_status = np.dtype([("TagType", "<u4"),
("TagLength", "<u4"),
("SatelliteID", "<u8"),
("SatelliteName", "S256"),
("NominalLongitude", "<f8"),
("SatelliteCondition", "<u4"),
("TimeOffset", "<f8")])
image_acquisition = np.dtype([("TagType", "<u4"),
("TagLength", "<u4"),
("Status", "<u4"),
("StartDelay", "<i4"),
("Cel", "<f8")])
prologue = np.dtype([("SatelliteStatus", satellite_status),
("ImageAcquisition", image_acquisition, (10, )),
("ImageCalibration", "<i4", (10, 1024))])
[docs]
def recarray2dict(arr):
"""Change record array to a dictionary."""
res = {}
for dtuple in arr.dtype.descr:
key = dtuple[0]
ntype = dtuple[1]
data = arr[key]
if isinstance(ntype, list):
res[key] = recarray2dict(data)
else:
res[key] = data
return res
[docs]
class HRITGOMSPrologueFileHandler(HRITFileHandler):
"""GOMS HRIT format reader."""
def __init__(self, filename, filename_info, filetype_info):
"""Initialize the reader."""
super(HRITGOMSPrologueFileHandler, self).__init__(filename, filename_info,
filetype_info,
(goms_hdr_map,
goms_variable_length_headers,
goms_text_headers))
self.prologue = {}
self.read_prologue()
[docs]
def read_prologue(self):
"""Read the prologue metadata."""
with open(self.filename, "rb") as fp_:
fp_.seek(self.mda["total_header_length"])
data = np.fromfile(fp_, dtype=prologue, count=1)[0]
self.prologue.update(recarray2dict(data))
self.process_prologue()
[docs]
def process_prologue(self):
"""Reprocess prologue to correct types."""
radiometric_processing = np.dtype([("TagType", "<u4"),
("TagLength", "<u4"),
("RPSummary",
[("Impulse", "<u4"),
("IsStrNoiseCorrection", "<u4"),
("IsOptic", "<u4"),
("IsBrightnessAligment", "<u4")]),
("OpticCorrection",
[("Degree", "<i4"),
("A", "<f8", (16, ))]),
("RPQuality",
[("EffDinRange", "<f8"),
("EathDarkening", "<f8"),
("Zone", "<f8"),
("Impulse", "<f8"),
("Group", "<f8"),
("DefectCount", "<u4"),
("DefectProcent", "<f8"),
("S_Noise_DT_Preflight", "<f8"),
("S_Noise_DT_Bort", "<f8"),
("S_Noise_DT_Video", "<f8"),
("S_Noise_DT_1_5", "<f8"),
("CalibrStability", "<f8"),
("TemnSKO", "<f8", (2, )),
("StructSKO", "<f8", (2, )),
("Struct_1_5", "<f8"),
("Zone_1_ 5", "<f8"),
("RadDif", "<f8")])])
geometric_processing = np.dtype([("TagType", "<u4"),
("TagLength", "<u4"),
("TGeomNormInfo",
[("IsExist", "<u4"),
("IsNorm", "<u4"),
("SubLon", "<f8"),
("TypeProjection", "<u4"),
("PixInfo", "<f8", (4, ))]),
("SatInfo",
[("TISO",
[("T0", "<f8"),
("dT", "<f8"),
("ASb", "<f8"),
("Evsk", "<f8", (3, 3, 4)),
("ARx", "<f8", (4, )),
("ARy", "<f8", (4, )),
("ARz", "<f8", (4, )),
("AVx", "<f8", (4, )),
("AVy", "<f8", (4, )),
("AVz", "<f8", (4, ))]),
("Type", "<i4")]),
("TimeProcessing", "<f8"),
("ApriorAccuracy", "<f8"),
("RelativeAccuracy", "<f8", (2, ))])
epilogue = np.dtype([("RadiometricProcessing", radiometric_processing, (10, )),
("GeometricProcessing", geometric_processing, (10, ))])
# FIXME: Add rest of the epilogue
[docs]
class HRITGOMSEpilogueFileHandler(HRITFileHandler):
"""GOMS HRIT format reader."""
def __init__(self, filename, filename_info, filetype_info):
"""Initialize the reader."""
super(HRITGOMSEpilogueFileHandler, self).__init__(filename, filename_info,
filetype_info,
(goms_hdr_map,
goms_variable_length_headers,
goms_text_headers))
self.epilogue = {}
self.read_epilogue()
[docs]
def read_epilogue(self):
"""Read the prologue metadata."""
with open(self.filename, "rb") as fp_:
fp_.seek(self.mda["total_header_length"])
data = np.fromfile(fp_, dtype=epilogue, count=1)[0]
self.epilogue.update(recarray2dict(data))
C1 = 1.19104273e-5
C2 = 1.43877523
# Defined in MSG Level 1.5 Image Data Format Description
# https://www-cdn.eumetsat.int/files/2020-05/pdf_ten_05105_msg_img_data.pdf
SPACECRAFTS = {19001: "Electro-L N1",
19002: "Electro-L N2",
19003: "Electro-L N3"}
[docs]
class HRITGOMSFileHandler(HRITFileHandler):
"""GOMS HRIT format reader."""
def __init__(self, filename, filename_info, filetype_info,
prologue, epilogue):
"""Initialize the reader."""
super(HRITGOMSFileHandler, self).__init__(filename, filename_info,
filetype_info,
(goms_hdr_map,
goms_variable_length_headers,
goms_text_headers))
self.prologue = prologue.prologue
self.epilogue = epilogue.epilogue
self.chid = self.mda["spectral_channel_id"]
sublon = self.epilogue["GeometricProcessing"]["TGeomNormInfo"]["SubLon"]
sublon = sublon[self.chid]
self.mda["projection_parameters"]["SSP_longitude"] = np.rad2deg(sublon)
self.mda["orbital_parameters"]["satellite_nominal_longitude"] = np.rad2deg(
self.prologue["SatelliteStatus"]["NominalLongitude"])
satellite_id = self.prologue["SatelliteStatus"]["SatelliteID"]
self.platform_name = SPACECRAFTS[satellite_id]
[docs]
def get_dataset(self, key, info):
"""Get the data from the files."""
res = super(HRITGOMSFileHandler, self).get_dataset(key, info)
res = self.calibrate(res, key["calibration"])
res.attrs["units"] = info["units"]
res.attrs["standard_name"] = info["standard_name"]
res.attrs["wavelength"] = info["wavelength"]
res.attrs["platform_name"] = self.platform_name
res.attrs["sensor"] = "msu-gs"
res.attrs["orbital_parameters"] = {
"satellite_nominal_longitude": self.mda["orbital_parameters"]["satellite_nominal_longitude"],
"satellite_nominal_latitude": 0.,
"projection_longitude": self.mda["projection_parameters"]["SSP_longitude"],
"projection_latitude": 0.,
"projection_altitude": 35785831.00
}
return res
[docs]
def calibrate(self, data, calibration):
"""Calibrate the data."""
tic = datetime.now()
if calibration == "counts":
res = data
elif calibration in ["radiance", "brightness_temperature"]:
res = self._calibrate(data)
else:
raise NotImplementedError("Don't know how to calibrate to " +
str(calibration))
res.attrs["standard_name"] = calibration
res.attrs["calibration"] = calibration
logger.debug("Calibration time " + str(datetime.now() - tic))
return res
[docs]
@staticmethod
def _getitem(block, lut):
return lut[block]
[docs]
def _calibrate(self, data):
"""Visible/IR channel calibration."""
lut = self.prologue["ImageCalibration"][self.chid]
if abs(lut).max() > 16777216:
lut = lut.astype(np.float64)
else:
lut = lut.astype(np.float32)
lut /= 1000
lut[0] = np.nan
# Dask/XArray don't support indexing in 2D (yet).
res = data.data.map_blocks(self._getitem, lut, dtype=lut.dtype)
res = xr.DataArray(res, dims=data.dims,
attrs=data.attrs, coords=data.coords)
res = res.where(data > 0)
return res
[docs]
def get_area_def(self, dsid):
"""Get the area definition of the band."""
pdict = {}
pdict["cfac"] = np.int32(self.mda["cfac"])
pdict["lfac"] = np.int32(self.mda["lfac"])
pdict["coff"] = np.float32(self.mda["coff"])
pdict["loff"] = np.float32(self.mda["loff"])
pdict["a"] = 6378169.00
pdict["b"] = 6356583.80
pdict["h"] = 35785831.00
pdict["scandir"] = "N2S"
pdict["ssp_lon"] = self.mda["projection_parameters"]["SSP_longitude"]
pdict["nlines"] = int(self.mda["number_of_lines"])
pdict["ncols"] = int(self.mda["number_of_columns"])
pdict["loff"] = pdict["nlines"] - pdict["loff"]
pdict["a_name"] = "geosgoms"
pdict["a_desc"] = "Electro-L/GOMS channel area"
pdict["p_id"] = "goms"
area_extent = get_area_extent(pdict)
area = get_area_definition(pdict, area_extent)
self.area = area
return area