Adds a common function for geostationary projection / area definition calculations

satpy/readers/_geos_area.py

@@ -0,0 +1,146 @@
+#!/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
+
+
+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__
+
+
+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
+
+
+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)
+    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
+
+
+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
+    """
+    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

satpy/readers/agri_l1.py

@@ -28,7 +28,7 @@
 import xarray as xr
 import dask.array as da
 from datetime import datetime
-from pyresample import geometry
+from satpy.readers._geos_area import get_area_extent, get_area_definition
 from satpy.readers.hdf5_utils import HDF5FileHandler
 
 logger = logging.getLogger(__name__)
@@ -124,46 +124,50 @@ def get_area_def(self, key):
         # Coordination Group for Meteorological Satellites LRIT/HRIT Global Specification
         # https://www.cgms-info.org/documents/cgms-lrit-hrit-global-specification-(v2-8-of-30-oct-2013).pdf
         res = key.resolution
-        coff = _COFF_list[_resolution_list.index(res)]
-        loff = _LOFF_list[_resolution_list.index(res)]
-        cfac = _CFAC_list[_resolution_list.index(res)]
-        lfac = _LFAC_list[_resolution_list.index(res)]
-        a = self.file_content['/attr/dEA'] * 1E3  # equator radius (m)
-        b = a * (1 - 1 / self.file_content['/attr/dObRecFlat'])  # polar radius (m)
-        h = self.file_content['/attr/NOMSatHeight']  # the altitude of satellite (m)
-
-        lon_0 = self.file_content['/attr/NOMCenterLon']
-        nlines = self.file_content['/attr/RegLength']
-        ncols = self.file_content['/attr/RegWidth']
-
-        cols = 0
-        left_x = (cols - coff) * (2.**16 / cfac)
-        cols += ncols - 1
-        right_x = (cols - coff) * (2.**16 / cfac)
-
-        lines = self.file_content['/attr/Begin Line Number']
-        upper_y = -(lines - loff) * (2.**16 / lfac)
-        lines = self.file_content['/attr/End Line Number']
-        lower_y = -(lines - loff) * (2.**16 / lfac)
-        area_extent = (np.deg2rad(left_x) * h, np.deg2rad(lower_y) * h,
-                       np.deg2rad(right_x) * h, np.deg2rad(upper_y) * h)
-
-        proj_dict = {'a': float(a),
-                     'b': float(b),
-                     'lon_0': float(lon_0),
-                     'h': float(h),
-                     'proj': 'geos',
-                     'units': 'm',
-                     'sweep': 'y'}
-
-        area = geometry.AreaDefinition(
-            self.filename_info['observation_type'],
-            "AGRI {} area".format(self.filename_info['observation_type']),
-            'FY-4A',
-            proj_dict,
-            ncols,
-            nlines,
-            area_extent)
+        pdict = {}
+        pdict['coff'] = _COFF_list[_resolution_list.index(res)]
+        pdict['loff'] = _LOFF_list[_resolution_list.index(res)]
+        pdict['cfac'] = _CFAC_list[_resolution_list.index(res)]
+        pdict['lfac'] = _LFAC_list[_resolution_list.index(res)]
+        pdict['a'] = self.file_content['/attr/dEA'] * 1E3  # equator radius (m)
+        pdict['b'] = pdict['a'] * (1 - 1 / self.file_content['/attr/dObRecFlat'])  # polar radius (m)
+        pdict['h'] = self.file_content['/attr/NOMSatHeight']  # the altitude of satellite (m)
+
+        pdict['ssp_lon'] = self.file_content['/attr/NOMCenterLon']
+        pdict['nlines'] = self.file_content['/attr/RegLength']
+        pdict['ncols'] = self.file_content['/attr/RegWidth']
+
+        pdict['scandir'] = 'S2N'
+
+        b500 = ['C02']
+        b1000 = ['C01', 'C03']
+        b2000 = ['C04', 'C05', 'C06', 'C07']
+
+        pdict['a_desc'] = "AGRI {} area".format(self.filename_info['observation_type'])
+
+        if (key.name in b500):
+            pdict['a_name'] = self.filename_info['observation_type']+'_500m'
+            pdict['p_id'] = 'FY-4A, 500m'
+        elif (key.name in b1000):
+            pdict['a_name'] = self.filename_info['observation_type']+'_1000m'
+            pdict['p_id'] = 'FY-4A, 1000m'
+        elif (key.name in b2000):
+            pdict['a_name'] = self.filename_info['observation_type']+'_2000m'
+            pdict['p_id'] = 'FY-4A, 2000m'
+        else:
+            pdict['a_name'] = self.filename_info['observation_type']+'_2000m'
+            pdict['p_id'] = 'FY-4A, 4000m'
+
+        pdict['coff'] = pdict['coff'] + 0.5
+        pdict['nlines'] = pdict['nlines'] - 1
+        pdict['ncols'] = pdict['ncols'] - 1
+        pdict['loff'] = (pdict['loff'] - self.file_content['/attr/End Line Number'] + 0.5)
+        area_extent = get_area_extent(pdict)
+        area_extent = (area_extent[0] + 2000, area_extent[1], area_extent[2] + 2000, area_extent[3])
+
+        pdict['nlines'] = pdict['nlines'] + 1
+        pdict['ncols'] = pdict['ncols'] + 1
+        area = get_area_definition(pdict, area_extent)
 
         return area
 

satpy/readers/ahi_hsd.py

@@ -40,10 +40,10 @@
 import xarray as xr
 import warnings
 
-from pyresample import geometry
 from satpy import CHUNK_SIZE
 from satpy.readers.file_handlers import BaseFileHandler
 from satpy.readers.utils import get_geostationary_mask, np2str, get_earth_radius
+from satpy.readers._geos_area import get_area_extent, get_area_definition
 
 AHI_CHANNEL_NAMES = ("1", "2", "3", "4", "5",
                      "6", "7", "8", "9", "10",
@@ -316,45 +316,29 @@ def get_dataset(self, key, info):
 
     def get_area_def(self, dsid):
         del dsid
-        cfac = np.uint32(self.proj_info['CFAC'])
-        lfac = np.uint32(self.proj_info['LFAC'])
-        coff = np.float32(self.proj_info['COFF'])
-        loff = np.float32(self.proj_info['LOFF'])
-        a = float(self.proj_info['earth_equatorial_radius'] * 1000)
-        h = float(self.proj_info['distance_from_earth_center'] * 1000 - a)
-        b = float(self.proj_info['earth_polar_radius'] * 1000)
-        lon_0 = float(self.proj_info['sub_lon'])
-        nlines = int(self.data_info['number_of_lines'])
-        ncols = int(self.data_info['number_of_columns'])
-
-        # count starts at 1
-        cols = 1 - 0.5
-        left_x = (cols - coff) * (2.**16 / cfac)
-        cols += ncols
-        right_x = (cols - coff) * (2.**16 / cfac)
-
-        lines = (self.segment_number - 1) * nlines + 1 - 0.5
-        upper_y = -(lines - loff) * (2.**16 / lfac)
-        lines += nlines
-        lower_y = -(lines - loff) * (2.**16 / lfac)
-        area_extent = (np.deg2rad(left_x) * h, np.deg2rad(lower_y) * h,
-                       np.deg2rad(right_x) * h, np.deg2rad(upper_y) * h)
-
-        proj_dict = {'a': float(a),
-                     'b': float(b),
-                     'lon_0': float(lon_0),
-                     'h': float(h),
-                     'proj': 'geos',
-                     'units': 'm'}
-
-        area = geometry.AreaDefinition(
-            self.observation_area,
-            "AHI {} area".format(self.observation_area),
-            'geosh8',
-            proj_dict,
-            ncols,
-            nlines,
-            area_extent)
+
+        pdict = {}
+        pdict['cfac'] = np.uint32(self.proj_info['CFAC'])
+        pdict['lfac'] = np.uint32(self.proj_info['LFAC'])
+        pdict['coff'] = np.float32(self.proj_info['COFF'])
+        pdict['loff'] = -np.float32(self.proj_info['LOFF']) + 1
+        pdict['a'] = float(self.proj_info['earth_equatorial_radius'] * 1000)
+        pdict['h'] = float(self.proj_info['distance_from_earth_center'] * 1000 - pdict['a'])
+        pdict['b'] = float(self.proj_info['earth_polar_radius'] * 1000)
+        pdict['ssp_lon'] = float(self.proj_info['sub_lon'])
+        pdict['nlines'] = int(self.data_info['number_of_lines'])
+        pdict['ncols'] = int(self.data_info['number_of_columns'])
+        pdict['scandir'] = 'N2S'
+
+        pdict['loff'] = pdict['loff'] + (self.segment_number * pdict['nlines'])
+
+        aex = get_area_extent(pdict)
+
+        pdict['a_name'] = self.observation_area
+        pdict['a_desc'] = "AHI {} area".format(self.observation_area)
+        pdict['p_id'] = 'geosh8'
+
+        area = get_area_definition(pdict, aex)
 
         self.area = area
         return area