#!/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 old NWCSAF/Geo (v2013 and earlier) cloud product format.
References:
- The NWCSAF GEO 2013 products documentation:
http://www.nwcsaf.org/web/guest/archive - Search for Code "ICD/3"; Type
"MSG" and the box to the right should say 'Status' (which means any
status). Version 7.0 seems to be for v2013
http://www.nwcsaf.org/aemetRest/downloadAttachment/2623
"""
import logging
from datetime import datetime
import h5py
import numpy as np
from pyresample.geometry import AreaDefinition
from satpy.readers.hdf5_utils import HDF5FileHandler
logger = logging.getLogger(__name__)
PLATFORM_NAMES = {"MSG1": "Meteosat-8",
"MSG2": "Meteosat-9",
"MSG3": "Meteosat-10",
"MSG4": "Meteosat-11", }
[docs]
class Hdf5NWCSAF(HDF5FileHandler):
"""NWCSAF MSG hdf5 reader."""
def __init__(self, filename, filename_info, filetype_info):
"""Init method."""
super(Hdf5NWCSAF, self).__init__(filename, filename_info,
filetype_info)
self.cache = {}
[docs]
def get_dataset(self, dataset_id, ds_info):
"""Load a dataset."""
file_key = ds_info.get("file_key", dataset_id["name"])
data = self[file_key]
nodata = None
if "SCALING_FACTOR" in data.attrs and "OFFSET" in data.attrs:
dtype = np.dtype(data.data)
if dataset_id["name"] in ["ctth_alti"]:
data.attrs["valid_range"] = (0, 27000)
data.attrs["_FillValue"] = np.nan
if dataset_id["name"] in ["ctth_alti", "ctth_pres", "ctth_tempe", "ctth_effective_cloudiness"]:
dtype = np.dtype("float32")
nodata = 255
if dataset_id["name"] in ["ct"]:
data.attrs["valid_range"] = (0, 20)
data.attrs["_FillValue"] = 255
# data.attrs['palette_meanings'] = list(range(21))
attrs = data.attrs
scaled_data = (data * data.attrs["SCALING_FACTOR"] + data.attrs["OFFSET"]).astype(dtype)
if nodata:
scaled_data = scaled_data.where(data != nodata)
scaled_data = scaled_data.where(scaled_data >= 0)
data = scaled_data
data.attrs = attrs
for key in list(data.attrs.keys()):
val = data.attrs[key]
if isinstance(val, h5py.h5r.Reference):
del data.attrs[key]
return data
[docs]
def get_area_def(self, dsid):
"""Get the area definition of the datasets in the file."""
if dsid["name"].endswith("_pal"):
raise NotImplementedError
cfac = self.file_content["/attr/CFAC"]
lfac = self.file_content["/attr/LFAC"]
coff = self.file_content["/attr/COFF"]
loff = self.file_content["/attr/LOFF"]
numcols = int(self.file_content["/attr/NC"])
numlines = int(self.file_content["/attr/NL"])
aex = get_area_extent(cfac, lfac, coff, loff, numcols, numlines)
pname = self.file_content["/attr/PROJECTION_NAME"]
proj = {}
if pname.startswith("GEOS"):
proj["proj"] = "geos"
proj["a"] = "6378169.0"
proj["b"] = "6356583.8"
proj["h"] = "35785831.0"
proj["lon_0"] = str(float(pname.split("<")[1][:-1]))
else:
raise NotImplementedError("Only geos projection supported yet.")
area_def = AreaDefinition(self.file_content["/attr/REGION_NAME"],
self.file_content["/attr/REGION_NAME"],
pname,
proj,
numcols,
numlines,
aex)
return area_def
@property
def start_time(self):
"""Return the start time of the object."""
return datetime.strptime(self.file_content["/attr/IMAGE_ACQUISITION_TIME"], "%Y%m%d%H%M")
[docs]
def get_area_extent(cfac, lfac, coff, loff, numcols, numlines):
"""Get the area extent from msg parameters."""
xur = (numcols - coff) * 2 ** 16 / (cfac * 1.0)
xur = np.deg2rad(xur) * 35785831.0
xll = (-1 - coff) * 2 ** 16 / (cfac * 1.0)
xll = np.deg2rad(xll) * 35785831.0
xres = (xur - xll) / numcols
xur, xll = xur - xres / 2, xll + xres / 2
yll = (numlines - loff) * 2 ** 16 / (-lfac * 1.0)
yll = np.deg2rad(yll) * 35785831.0
yur = (-1 - loff) * 2 ** 16 / (-lfac * 1.0)
yur = np.deg2rad(yur) * 35785831.0
yres = (yur - yll) / numlines
yll, yur = yll + yres / 2, yur - yres / 2
return xll, yll, xur, yur