New helper function for standardized GEO area definitions

[strandgren]

Feature Request

Is your feature request related to a problem? Please describe. The list of standardized area definitions in areas.yaml is growing, with e.g. different resolutions, operating modes and sub-satellite points. Currently 36 standardized area definitions are hard coded in areas.yaml, with quite a bit of duplicate information.

This also allows for limited flexibility. A recent example is the change in sub-satellite point of the MSG SEVIRI IODC service from 41.5E to 45.5E. An easy solution would be to just change the lon_0 in the IODC-related area definitions. This would however lead to errors if users use this for older data for which the lon_0 is still 41.5E, for which the sub-satellite point was still 41.5E. Of course these is also the possibility to add a new set of (five) area definitions for the IODC service at this new longitude. But that's what I would like to avoid.

Describe the solution you'd like An extended or new helper function that compute/construct the area definitions dynamically when called, instead of just looking for hard coded entries in areas.yaml.

Example Code Below is a small prototype of the proposed helper function together with a dictionary containing all information currently used for the 36 standardized SEVIRI, FCI and ABI area definitions.

from pyresample.area_config import create_area_def

AREA_INFO = {}

# MSG SEVIRI
AREA_INFO['seviri'] = {}
AREA_INFO['seviri']['fes'] = {'name': "Full Earth Scanning", 'projection': {}, 'shape': {}, 'area_extent': {}}
AREA_INFO['seviri']['fes']['projection']['proj'] = 'geos'
AREA_INFO['seviri']['fes']['projection']['a'] = 6378169.0
AREA_INFO['seviri']['fes']['projection']['b'] = 6356583.8
AREA_INFO['seviri']['fes']['projection']['h'] = 35785831.0
AREA_INFO['seviri']['fes']['shape']['1km'] = (11136, 11136)
AREA_INFO['seviri']['fes']['shape']['3km'] = (3712, 3712)
AREA_INFO['seviri']['fes']['shape']['9km'] = (1237, 1237)
AREA_INFO['seviri']['fes']['shape']['48km'] = (232, 232)
AREA_INFO['seviri']['fes']['area_extent']['1km'] = (-5571248.412732527, -5566247.740968115, 5566247.740968115, 5571248.412732527)
AREA_INFO['seviri']['fes']['area_extent']['3km'] = (-5570248.686685662, -5567248.28340708, 5567248.28340708, 5570248.686685662)
AREA_INFO['seviri']['fes']['area_extent']['9km'] = (-5567248.28351984, -5567248.28340708, 5567248.28340708, 5567248.28351984)
AREA_INFO['seviri']['fes']['area_extent']['48km'] = AREA_INFO['seviri']['fes']['area_extent']['3km']

AREA_INFO['seviri']['rss'] = {'name': "Rapid Scanning Service", 'projection': {}, 'shape': {}, 'area_extent': {}}
AREA_INFO['seviri']['rss']['projection'] = AREA_INFO['seviri']['fes']['projection']
AREA_INFO['seviri']['rss']['shape'] = AREA_INFO['seviri']['fes']['shape']
AREA_INFO['seviri']['rss']['area_extent'] = AREA_INFO['seviri']['fes']['area_extent']

AREA_INFO['seviri']['iodc'] = {'name': "Indian Ocean Data Coverage", 'projection': {}, 'shape': {}, 'area_extent': {}}
AREA_INFO['seviri']['iodc']['projection'] = AREA_INFO['seviri']['fes']['projection']
AREA_INFO['seviri']['iodc']['shape'] = AREA_INFO['seviri']['fes']['shape']
AREA_INFO['seviri']['iodc']['area_extent'] = AREA_INFO['seviri']['fes']['area_extent']

# MTG FCI
AREA_INFO['fci'] = {}
AREA_INFO['fci']['fdss'] = {'name': "Full Disk Scanning Service", 'projection': {}, 'shape': {}, 'area_extent': {}}
AREA_INFO['fci']['fdss']['projection']['proj'] = 'geos'
AREA_INFO['fci']['fdss']['projection']['h'] = 35786400
AREA_INFO['fci']['fdss']['projection']['x_0'] = 0
AREA_INFO['fci']['fdss']['projection']['y_0'] = 0
AREA_INFO['fci']['fdss']['projection']['ellps'] = 'WGS84'
AREA_INFO['fci']['fdss']['projection']['no_defs'] = 'null'
AREA_INFO['fci']['fdss']['shape']['1km'] = (11136, 11136)
AREA_INFO['fci']['fdss']['shape']['2km'] = (5568, 5568)
AREA_INFO['fci']['fdss']['shape']['6km'] = (1856, 1856)
AREA_INFO['fci']['fdss']['shape']['32km'] = (348, 348)
AREA_INFO['fci']['fdss']['area_extent']['1km'] = (-5567999.998577303, -5567999.998577303, 5567999.998527619,  5567999.998527619)
AREA_INFO['fci']['fdss']['area_extent']['2km'] = (-5567999.994200589, -5567999.994200589, 5567999.994206558,  5567999.994206558)
AREA_INFO['fci']['fdss']['area_extent']['6km'] = AREA_INFO['fci']['fdss']['area_extent']['2km']
AREA_INFO['fci']['fdss']['area_extent']['32km'] = AREA_INFO['fci']['fdss']['area_extent']['2km']

# GEO ABI
AREA_INFO['abi'] = {}
AREA_INFO['abi']['f'] = {'name': "Full Disk", 'projection': {}, 'shape': {}, 'area_extent': {}}
AREA_INFO['abi']['f']['projection']['proj'] = 'geos'
AREA_INFO['abi']['f']['projection']['sweep'] = 'x'
AREA_INFO['abi']['f']['projection']['h'] = 35786023
AREA_INFO['abi']['f']['projection']['x_0'] = 0
AREA_INFO['abi']['f']['projection']['y_0'] = 0
AREA_INFO['abi']['f']['projection']['ellps'] = 'GRS80'
AREA_INFO['abi']['f']['projection']['no_defs'] = 'null'
AREA_INFO['abi']['f']['projection']['type'] = 'crs'
AREA_INFO['abi']['f']['shape']['500m'] = (21696, 21696)
AREA_INFO['abi']['f']['shape']['1km'] = (10848, 10848)
AREA_INFO['abi']['f']['shape']['2km'] = (5424, 5424)
AREA_INFO['abi']['f']['area_extent']['500m'] = (-5434894.885056, -5434894.885056, 5434894.885056, 5434894.885056)
AREA_INFO['abi']['f']['area_extent']['1km'] = AREA_INFO['abi']['f']['area_extent']['500m']
AREA_INFO['abi']['f']['area_extent']['2km'] = AREA_INFO['abi']['f']['area_extent']['500m']

AREA_INFO['abi']['c'] = {'name': "CONUS", 'projection': {}, 'shape': {}, 'area_extent': {}}
AREA_INFO['abi']['c']['projection'] = AREA_INFO['abi']['f']['projection']
AREA_INFO['abi']['c']['shape']['500m'] = (6000, 10000)
AREA_INFO['abi']['c']['shape']['1km'] = (3000, 5000)
AREA_INFO['abi']['c']['shape']['2km'] = (1500, 2500)
AREA_INFO['abi']['c']['area_extent']['500m'] = (-3627271.29128, 1583173.65752, 1382771.92872, 4589199.58952)
AREA_INFO['abi']['c']['area_extent']['1km'] = AREA_INFO['abi']['c']['area_extent']['500m']
AREA_INFO['abi']['c']['area_extent']['2km'] = AREA_INFO['abi']['c']['area_extent']['500m']

AREA_INFO['abi']['p'] = {'name': "PACUS", 'projection': {}, 'shape': {}, 'area_extent': {}}
AREA_INFO['abi']['p']['projection'] = AREA_INFO['abi']['f']['projection']
AREA_INFO['abi']['p']['shape']['500m'] = (6000, 10000)
AREA_INFO['abi']['p']['shape']['1km'] = (3000, 5000)
AREA_INFO['abi']['p']['shape']['2km'] = (1500, 2500)
AREA_INFO['abi']['p']['area_extent']['500m'] = (-2505021.61, 1583173.65752, 2505021.61, 4589199.58952)
AREA_INFO['abi']['p']['area_extent']['1km'] = AREA_INFO['abi']['p']['area_extent']['500m']
AREA_INFO['abi']['p']['area_extent']['2km'] = AREA_INFO['abi']['p']['area_extent']['500m']

def new_helper_function(platform, intrument, service, resolution, lon_0):
    if lon_0 >= 0.0:
        lon0_str = f'E{int(lon_0 * 10):04d}'
    else:
        lon0_str = f'W{int(abs(lon_0) * 10):04d}'

    area_id = f'{platform}_{intrument}_{service}_{resolution}_{lon0_str}'.lower()
    description = f'Area definition for {platform.upper()} {intrument.upper()} {AREA_INFO[intrument][service]["name"]} mode at {resolution} resolution with sub-satellite longitude at {lon_0} degrees East.'
    proj_dict = AREA_INFO[intrument][service]['projection']
    shape = AREA_INFO[intrument][service]['shape'][resolution]
    area_extent = AREA_INFO[intrument][service]['area_extent'][resolution]

    proj_dict['lon_0'] = lon_0

    return create_area_def(area_id, proj_dict, description=description, shape=shape, area_extent=area_extent)

# Two example calls
adef = new_helper_function('msg', 'seviri', 'fes', '3km', 45.5)
adef = new_helper_function('goes_west', 'abi', 'c', '500m', -137)

Describe any changes to existing user workflow Either this helper function could be used stand-alone. Or, if possible, it could be integrated in satpy.resample.get_area_def with the following workflow:

satpy.resample import get_area_def
adef = get_area_def(area_id)
if adef is None:  # No area definition with that name found in areas.yaml
    check if area_id matches an expected pre-defined format
    chop up area_id into platform, instrument, service, resolution and sub-satellite point longitude assuming a pre-defined format
    adef = new_helper_function(platform, instrument, service, resolution, lon_0)

[djhoese]

As discussed in the meeting today, here are some options that might be good ideas:

1. Auto-generate the YAML file with some utility script. We get the benefit of python for-loops and shared information (ex. extents that are the same between multiple areas) while still having YAML that can be saved in git, linked to on github (for user support/questions), and other tools can use the YAML directly instead of having to use dynamically created areas (ex. a page on the sphinx docs which lists the areas and shows the user a boundary).
2. A utility function that takes the information the user knows (satellite, instrument, resolution, etc) and finds the best area of the satpy builtin "standard" areas that matches that description. Or perhaps this function returns a list of possible ones.
3. A utility function that lists all the names of the areas in the builtin Satpy areas file. Could be the same function as described above, but with no additional parameters specified.
4. Refactor area loading in Satpy and/or pyresample to load from multiple files where YAML files are separated based on some property (satellite or instrument or something else).

Also discussed was the possibility of using YAML aliases to reuse information across multiple definitions. This makes things harder to read but means that things that are easy to mistype aren't duplicated (ex. area extents).

I also brought up the possibility of adding something to satpy's pre-commit config to detect changes in the auto-generate python script and then run it to generate the YAML on commit. This way the YAML should never be out of sync with the script.