pytroll / pyresample

Geospatial image resampling in Python
http://pyresample.readthedocs.org
GNU Lesser General Public License v3.0
349 stars 94 forks source link

get_geostationary_bounding_box* contains duplicated vertices at the equator #474

Closed ghiggi closed 1 year ago

ghiggi commented 1 year ago

The get_geostationary_bounding_box returns duplicated vertices at the equator. The code snippet below shows the "bad" behaviour.

Code Sample, a minimal, complete, and verifiable piece of code

from pyresample import AreaDefinition
from pyresample.geometry import (
    get_geostationary_bounding_box_in_proj_coords,
    get_geostationary_bounding_box_in_lonlats
)
# Define geostationary projection
proj_dict_geostationary = {
        'proj_id': "dummy_geo_projection",
        "area_id" : 'dummy_geo_projection',
        "description" : 'geostationary projection',
        "projection" : {'a': 6378169.00,  'b': 6356583.80, 'h': 35785831.00, 'lon_0': 0, 'proj': 'geos'},
        "area_extent" : (-5500000., -5500000., 5500000., 5500000.),
        "width" : 100,
        "height" : 100,
}
areadef = AreaDefinition(**proj_dict_geostationary)
geos_area = areadef

#### If nb_points is 4, it requires only 2 points at the equator
nb_points = 4
x, y = get_geostationary_bounding_box_in_proj_coords(geos_area=geos_area,
                                                     nb_points=nb_points) 
print(x)
print(y) # all 0 

lons, lats = get_geostationary_bounding_box_in_lonlats(geos_area=geos_area,
                                                       nb_points=nb_points)
print(lons)
print(lats) # all 0 

#### get_geostationary_bounding_box_in_proj_coords  returns duplicated values at the equator 
lons, lats = get_geostationary_bounding_box_in_lonlats(geos_area=geos_area,
                                                       nb_points=6)
print(lons)
print(lats)  
lons[2:4] # [79.23372832, 79.23372832]
lats[2:4]  # [-0.00000000e+00, 0.00000000e+00]

Expected Output

The get_geostationary_bounding_box functions in case nbpoints=4 should return the top, left, bottom and right "extent" of the geo disk, instead of the current left and right duplicates.