AreaDefinition.__contains__ raises ValueError

[gerritholl]

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

# Your code here
from satpy import Scene
from satpy.utils import debug_on
from glob import glob
debug_on()
sc = Scene(filenames=glob("/data/gholl/cache/fogtools/abi/2017/03/14/00/06/*/*.nc"), reader="abi_l1b")
sc.load(["C10"])
print((0, -40) in sc["C10"].attrs["area"])

Problem description

The code fails with a ValueError: Illegal (lon, lat) coordinates: (inf, inf), thrown by pyresample's spherical_geometry module. This happens both when the coordinates are inside or outside the area.

Expected Output

I would expect the output to be either True or False depending on to whether the coordinates (0, 40) are within the area, in other words, visible by ABI. I can still achieve this by attempting area.get_xy_from_lonlat(...) and catching the possibly resulting ValueError, as follows:

try:
    sc["C10"].attrs["area"].get_xy_from_lonlat(45, 0)
except ValueError:
    return False
else:
    return True

I would expect an OverflowError or ValueError if lat/lon passed are out of valid range, a ValueError if they are nan, and a TypeError if they are not numbers.

Actual Result, Traceback if applicable

[DEBUG: 2020-04-15 17:19:48 : satpy.scene] Setting 'PPP_CONFIG_DIR' to '/home/gholl/checkouts/satpy/satpy/etc/'
[DEBUG: 2020-04-15 17:19:48 : satpy.readers] Reading ['/data/gholl/miniconda3/envs/py38/lib/python3.8/site-packages/satpy/etc/readers/abi_l1b.yaml', '/home/gholl/checkouts/satpy/satpy/etc/readers/abi_l1b.yaml']                                                                       
[DEBUG: 2020-04-15 17:19:48 : satpy.readers.yaml_reader] Assigning to abi_l1b: ['/data/gholl/cache/fogtools/abi/2017/03/14/00/06/2/OR_ABI-L1b-RadF-M3C02_G16_s20170730006100_e20170730016467_c20170730016509.nc', '/data/gholl/cache/fogtools/abi/2017/03/14/00/06/3/OR_ABI-L1b-RadF-M3C03_G16_s20170730006100_e20170730016467_c20170730016517.nc', '/data/gholl/cache/fogtools/abi/2017/03/14/00/06/4/OR_ABI-L1b-RadF-M3C04_G16_s20170730006100_e20170730016467_c20170730016511.nc', '/data/gholl/cache/fogtools/abi/2017/03/14/00/06/5/OR_ABI-L1b-RadF-M3C05_G16_s20170730006100_e20170730016467_c20170730016515.nc', '/data/gholl/cache/fogtools/abi/2017/03/14/00/06/6/OR_ABI-L1b-RadF-M3C06_G16_s20170730006100_e20170730016472_c20170730016515.nc', '/data/gholl/cache/fogtools/abi/2017/03/14/00/06/7/OR_ABI-L1b-RadF-M3C07_G16_s20170730006100_e20170730016478_c20170730016513.nc', '/data/gholl/cache/fogtools/abi/2017/03/14/00/06/8/OR_ABI-L1b-RadF-M3C08_G16_s20170730006100_e20170730016467_c20170730016516.nc', '/data/gholl/cache/fogtools/abi/2017/03/14/00/06/9/OR_ABI-L1b-RadF-M3C09_G16_s20170730006100_e20170730016472_c20170730016531.nc', '/data/gholl/cache/fogtools/abi/2017/03/14/00/06/10/OR_ABI-L1b-RadF-M3C10_G16_s20170730006100_e20170730016478_c20170730016523.nc', '/data/gholl/cache/fogtools/abi/2017/03/14/00/06/11/OR_ABI-L1b-RadF-M3C11_G16_s20170730006100_e20170730016467_c20170730016525.nc', '/data/gholl/cache/fogtools/abi/2017/03/14/00/06/13/OR_ABI-L1b-RadF-M3C13_G16_s20170730006100_e20170730016478_c20170730016532.nc', '/data/gholl/cache/fogtools/abi/2017/03/14/00/06/14/OR_ABI-L1b-RadF-M3C14_G16_s20170730006100_e20170730016467_c20170730016530.nc', '/data/gholl/cache/fogtools/abi/2017/03/14/00/06/15/OR_ABI-L1b-RadF-M3C15_G16_s20170730006100_e20170730016473_c20170730016532.nc', '/data/gholl/cache/fogtools/abi/2017/03/14/00/06/16/OR_ABI-L1b-RadF-M3C16_G16_s20170730006100_e20170730016478_c20170730016532.nc']                                                                    
[DEBUG: 2020-04-15 17:19:48 : satpy.composites] Looking for composites config file abi.yaml
[DEBUG: 2020-04-15 17:19:48 : satpy.composites] Looking for composites config file visir.yaml
/data/gholl/miniconda3/envs/py38/lib/python3.8/site-packages/pyproj/crs/crs.py:543: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems
  proj_string = self.to_proj4()
[DEBUG: 2020-04-15 17:19:48 : satpy.readers.abi_l1b] Reading in get_dataset C10.
[DEBUG: 2020-04-15 17:19:48 : satpy.readers.abi_l1b] Calibrating to brightness temperatures
Traceback (most recent call last):
  File "mwe35.py", line 7, in <module>
    print((0, -40) in sc["C10"].attrs["area"])
  File "/data/gholl/miniconda3/envs/py38/lib/python3.8/site-packages/pyresample/geometry.py", line 334, in __contains__
    corners = self.corners
  File "/data/gholl/miniconda3/envs/py38/lib/python3.8/site-packages/pyresample/geometry.py", line 322, in corners
    return [Coordinate(*self.get_lonlat(0, 0)),
  File "/data/gholl/miniconda3/envs/py38/lib/python3.8/site-packages/pyresample/spherical_geometry.py", line 55, in __init__
    raise ValueError('Illegal (lon, lat) coordinates: (%s, %s)'
ValueError: Illegal (lon, lat) coordinates: (inf, inf)

Versions of Python, package at hand and relevant dependencies

Python 3.8.2, pyresample 1.15.0.

[djhoese]

Are the input coordinates (the ones passed to contains) meant to be X/Y or lon/lat?

Side note: I didn't know this existed.

[gerritholl]

It's supposed to be (lon, lat). From https://pyresample.readthedocs.io/en/latest/geo_def.html:

It can be tested if a (lon, lat) point is inside a GeometryDefinition

>>> print((0, -90) in area_def)
True

It's missing from the API documentation though, perhaps due to special methods such as __contains__ being excluded in the sphinx configuration?

[gerritholl]

Although that doc applies to GeometryDefinition, I assume that AreaDefinition uses the same units (the docstring for AreaDefinition.__contains__ doesn't say).

[djhoese]

Is this full disk ABI data?

I have a feeling a check needs to be added to see if the corners are valid lon/lats:

pyresample/pyresample/geometry.py

Lines 318 to 339 in 3d6d38d

	@property
	def corners(self):
	"""Return the corners of the current area."""
	from pyresample.spherical_geometry import Coordinate
	return [Coordinate(*self.get_lonlat(0, 0)),
	Coordinate(*self.get_lonlat(0, -1)),
	Coordinate(*self.get_lonlat(-1, -1)),
	Coordinate(*self.get_lonlat(-1, 0))]
	def __contains__(self, point):
	"""Check if a point is inside the 4 corners of the current area.
	This uses great circle arcs as area boundaries.
	"""
	from pyresample.spherical_geometry import point_inside, Coordinate
	corners = self.corners
	if isinstance(point, tuple):
	return point_inside(Coordinate(*point), corners)
	else:
	return point_inside(point, corners)

The AreaDefinition is a subclass of BaseDefinition where these methods are defined.

[gerritholl]

Yes, it's full disk ABI, but I also encountered the problem with less-than-full-disk areas, probably whenever not all corners are valid lat/lons.

[djhoese]

probably whenever not all corners are valid lat/lons.

Exactly