#!/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/>.
"""Geostationary Projection / Area computations.
This module computes properties and area definitions for geostationary
satellites. It is designed to be a common module that can be called by
all geostationary satellite readers and uses commonly-included parameters
such as the CFAC/LFAC values, satellite position, etc, to compute the
correct area definition.
"""
import numpy as np
from pyresample import geometry
[docs]
def get_xy_from_linecol(line, col, offsets, factors):
"""Get the intermediate coordinates from line & col.
Intermediate coordinates are actually the instruments scanning angles.
"""
loff, coff = offsets
lfac, cfac = factors
x__ = float(col - coff) / (float(cfac) / 2 ** 16)
y__ = float(line - loff) / (float(lfac) / 2 ** 16)
return x__, y__
[docs]
def make_ext(ll_x, ur_x, ll_y, ur_y, h):
"""Create the area extent from computed ll and ur.
Args:
ll_x: The lower left x coordinate (m)
ur_x: The upper right x coordinate (m)
ll_y: The lower left y coordinate (m)
ur_y: The upper right y coordinate (m)
h: The satellite altitude above the Earth's surface
Returns:
aex: An area extent for the scene
"""
aex = (np.deg2rad(ll_x) * h, np.deg2rad(ll_y) * h,
np.deg2rad(ur_x) * h, np.deg2rad(ur_y) * h)
return aex
[docs]
def get_area_extent(pdict):
"""Get the area extent seen by a geostationary satellite.
Args:
pdict: A dictionary containing common parameters:
nlines: Number of lines in image
ncols: Number of columns in image
cfac: Column scaling factor
lfac: Line scaling factor
coff: Column offset factor
loff: Line offset factor
scandir: 'N2S' for standard (N->S), 'S2N' for inverse (S->N)
h: Altitude of satellite above the Earth's surface (m)
Returns:
aex: An area extent for the scene
"""
# count starts at 1
cols = 1 - 0.5
if pdict["scandir"] == "S2N":
lines = 0.5 - 1
scanmult = -1
else:
lines = 1 - 0.5
scanmult = 1
# Lower left x, y scanning angles in degrees
ll_x, ll_y = get_xy_from_linecol(lines * scanmult,
cols,
(pdict["loff"], pdict["coff"]),
(pdict["lfac"], pdict["cfac"]))
cols += pdict["ncols"]
lines += pdict["nlines"]
# Upper right x, y scanning angles in degrees
ur_x, ur_y = get_xy_from_linecol(lines * scanmult,
cols,
(pdict["loff"], pdict["coff"]),
(pdict["lfac"], pdict["cfac"]))
if pdict["scandir"] == "S2N":
ll_y *= -1
ur_y *= -1
# Convert degrees to radians and create area extent
aex = make_ext(ll_x=ll_x, ur_x=ur_x, ll_y=ll_y, ur_y=ur_y, h=pdict["h"])
return aex
[docs]
def get_area_definition(pdict, a_ext):
"""Get the area definition for a geo-sat.
Args:
pdict: A dictionary containing common parameters:
nlines: Number of lines in image
ncols: Number of columns in image
ssp_lon: Subsatellite point longitude (deg)
a: Earth equatorial radius (m)
b: Earth polar radius (m)
h: Platform height (m)
a_name: Area name
a_desc: Area description
p_id: Projection id
a_ext: A four element tuple containing the area extent (scan angle)
for the scene in radians
Returns:
a_def: An area definition for the scene
.. note::
The AreaDefinition `proj_id` attribute is being deprecated.
"""
proj_dict = {"a": float(pdict["a"]),
"b": float(pdict["b"]),
"lon_0": float(pdict["ssp_lon"]),
"h": float(pdict["h"]),
"proj": "geos",
"units": "m"}
a_def = geometry.AreaDefinition(
pdict["a_name"],
pdict["a_desc"],
pdict["p_id"],
proj_dict,
int(pdict["ncols"]),
int(pdict["nlines"]),
a_ext)
return a_def
[docs]
def sampling_to_lfac_cfac(sampling):
"""Convert angular sampling to line/column scaling factor (aka LFAC/CFAC).
Reference: `MSG Ground Segment LRIT HRIT Mission Specific Implementation`_,
Appendix E.2.
.. _MSG Ground Segment LRIT HRIT Mission Specific Implementation:
https://www-cdn.eumetsat.int/files/2020-04/pdf_ten_05057_spe_msg_lrit_hri.pdf
Args:
sampling: float
Angular sampling (rad)
Returns:
Line/column scaling factor (deg-1)
"""
return 2.0 ** 16 / np.rad2deg(sampling)
[docs]
def get_geos_area_naming(input_dict):
"""Get a dictionary containing formatted AreaDefinition naming.
Args:
input_dict: dict
Dictionary with keys `platform_name`, `instrument_name`, `service_name`, `service_desc`, `resolution` .
The resolution is expected in meters.
Returns:
area_naming_dict with `area_id`, `description` keys, values are strings.
.. note::
The AreaDefinition `proj_id` attribute is being deprecated and is therefore not formatted here.
An empty string is to be used until the attribute is fully removed.
"""
area_naming_dict = {}
resolution_strings = get_resolution_and_unit_strings(input_dict["resolution"])
area_naming_dict["area_id"] = "{}_{}_{}_{}{}".format(input_dict["platform_name"].lower(),
input_dict["instrument_name"].lower(),
input_dict["service_name"].lower(),
resolution_strings["value"],
resolution_strings["unit"]
)
area_naming_dict["description"] = "{} {} {} area definition " \
"with {} {} resolution".format(input_dict["platform_name"].upper(),
input_dict["instrument_name"].upper(),
input_dict["service_desc"],
resolution_strings["value"],
resolution_strings["unit"]
)
return area_naming_dict
[docs]
def get_resolution_and_unit_strings(resolution):
"""Get the resolution value and unit as strings.
If the resolution is larger than 1000 m, use kilometer as unit. If lower, use meter.
Args:
resolution: scalar
Resolution in meters.
Returns:
Dictionary with `value` and `unit` keys, values are strings.
"""
if resolution >= 1000:
return {"value": "{:.0f}".format(resolution*1e-3),
"unit": "km"}
return {"value": "{:.0f}".format(resolution),
"unit": "m"}