outer_boundary_corners and __contains__ raise ValueError on some areas
Checking the outer boundary corners of some areas fails with ValueError.  In this case, checking if a certain point is contained also fails.

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

```python
from pyresample.area_config import parse_area_file
area = parse_area_file(
    "/data/gholl/checkouts/satpy/satpy/etc/areas.yaml",
    "seviri_0deg")[0]
print((0, 0) in area)
print(area.outer_boundary_corners)
```
#### Problem description

raises ValueError, probably because the Earth is round and has no corners.

#### Expected Output

For the `outer_boundary_corners`, maybe NaN would be appropriate, or maybe a `ValueError` is actually correct.  However, I don't think `(lon, lat) in area` should raise a `ValueError` for any valid lon/lat (see also #268).  For valid input, this should always return either True or False.

#### Actual Result, Traceback if applicable

```
Traceback (most recent call last):
  File "/data/gholl/checkouts/protocode/mwe/outer-corners.py", line 5, in <module>
    print(area.outer_boundary_corners)
  File "/data/gholl/checkouts/pyresample/pyresample/geometry.py", line 2134, in outer_boundary_corners
    return [Coordinate(corner_lons[0], corner_lats[0]),
  File "/data/gholl/checkouts/pyresample/pyresample/spherical_geometry.py", line 50, 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

pyresample v1.22.0-8-gee7c1b9


diff --git a/.github/workflows/deploy.yaml b/.github/workflows/deploy.yaml
index 1bd1f7f..cf90d7b 100644
--- a/.github/workflows/deploy.yaml
+++ b/.github/workflows/deploy.yaml
@@ -46,10 +46,6 @@ jobs:
             os: ubuntu-latest
             python-version: '3.8'
             docker-image: manylinux2014_x86_64
-          - name: manylinux 32-bit
-            os: ubuntu-latest
-            python-version: '3.8'
-            docker-image: manylinux2014_i686
 
     steps:
       - uses: actions/checkout@v3
@@ -142,4 +138,3 @@ jobs:
           user: __token__
           password: ${{ secrets.pypi_password }}
           skip_existing: true
-
diff --git a/pyresample/geometry.py b/pyresample/geometry.py
index 9e64242..de07947 100644
--- a/pyresample/geometry.py
+++ b/pyresample/geometry.py
@@ -36,10 +36,9 @@ from pyresample import CHUNK_SIZE
 from pyresample._spatial_mp import Cartesian, Cartesian_MP, Proj_MP
 from pyresample.area_config import create_area_def
 from pyresample.boundary import AreaDefBoundary, Boundary, SimpleBoundary
-from pyresample.utils import (
-    check_slice_orientation,
+from pyresample.utils import check_slice_orientation, load_cf_area
+from pyresample.utils.proj4 import (
     get_geostationary_height,
-    load_cf_area,
     proj4_dict_to_str,
     proj4_radius_parameters,
 )
@@ -1556,7 +1555,7 @@ class AreaDefinition(_ProjectionDefinition):
         # Ensure an even number of vertices for side creation
         if (frequency % 2) != 0:
             frequency = frequency + 1
-        lons, lats = get_geostationary_bounding_box(self, nb_points=frequency)
+        lons, lats = get_geostationary_bounding_box_in_lonlats(self, nb_points=frequency)
         # Retrieve dummy sides for GEO (side1 and side3 always of length 2)
         side02_step = int(frequency / 2) - 1
         lon_sides = [lons[slice(0, side02_step + 1)],
@@ -2772,40 +2771,49 @@ def get_geostationary_angle_extent(geos_area):
     aeq = 1 - req ** 2 / (h ** 2)
     ap_ = 1 - rp ** 2 / (h ** 2)
 
-    # generate points around the north hemisphere in satellite projection
-    # make it a bit smaller so that we stay inside the valid area
-    xmax = np.arccos(np.sqrt(aeq))
-    ymax = np.arccos(np.sqrt(ap_))
-    return xmax, ymax
+    x_angle = np.arccos(np.sqrt(aeq))
+    y_angle = np.arccos(np.sqrt(ap_))
+    return x_angle, y_angle
 
 
 def get_geostationary_bounding_box_in_proj_coords(geos_area, nb_points=50):
     """Get the bbox in geos projection coordinates of the valid pixels inside `geos_area`.
 
-    Notes:
-    - The first and last element of the output vectors are equal.
-    - If nb_points is even, it will return x and y vectors of length nb_points + 1.
-
     Parameters
     ----------
     nb_points : Number of points on the polygon.
 
     """
-    xmax, ymax = get_geostationary_angle_extent(geos_area)
+    x, y = get_full_geostationary_bounding_box_in_proj_coords(geos_area, nb_points)
+    ll_x, ll_y, ur_x, ur_y = geos_area.area_extent
+
+    from shapely.geometry import Polygon
+    geo_bbox = Polygon(np.vstack((x, y)).T)
+    area_bbox = Polygon(((ll_x, ll_y), (ll_x, ur_y), (ur_x, ur_y), (ur_x, ll_y)))
+    intersection = area_bbox.intersection(geo_bbox)
+    try:
+        x, y = intersection.boundary.xy
+    except NotImplementedError:
+        return [], []
+    return np.asanyarray(x[:-1]), np.asanyarray(y[:-1])
+
+
+def get_full_geostationary_bounding_box_in_proj_coords(geos_area, nb_points=50):
+    """Get the bbox in geos projection coordinates of the full disk in `geos_area` projection.
+
+    Args:
+      nb_points: Number of points on the polygon
+    """
+    x_max_angle, y_max_angle = get_geostationary_angle_extent(geos_area)
     h = get_geostationary_height(geos_area.crs)
 
     # generate points around the north hemisphere in satellite projection
     # make it a bit smaller so that we stay inside the valid area
-    x = np.cos(np.linspace(-np.pi, 0, int(nb_points / 2.0) + 1)) * (xmax - 0.0001)
-    y = -np.sin(np.linspace(-np.pi, 0, int(nb_points / 2.0) + 1)) * (ymax - 0.0001)
-
-    ll_x, ll_y, ur_x, ur_y = geos_area.area_extent
-
+    points_around = np.linspace(-np.pi, np.pi, nb_points, endpoint=False)
+    x = np.cos(points_around) * (x_max_angle - 0.0001)
+    y = -np.sin(points_around) * (y_max_angle - 0.0001)
     x *= h
     y *= h
-    # We remove one element with [:-1] to avoid duplicate values at the equator
-    x = np.clip(np.concatenate([x[:-1], x[::-1]]), min(ll_x, ur_x), max(ll_x, ur_x))
-    y = np.clip(np.concatenate([y[:-1], -y]), min(ll_y, ur_y), max(ll_y, ur_y))
 
     return x, y
 
diff --git a/pyresample/gradient/__init__.py b/pyresample/gradient/__init__.py
index 68d5bf2..57c67c2 100644
--- a/pyresample/gradient/__init__.py
+++ b/pyresample/gradient/__init__.py
@@ -39,7 +39,7 @@ from pyresample import CHUNK_SIZE
 from pyresample.geometry import (
     AreaDefinition,
     SwathDefinition,
-    get_geostationary_bounding_box,
+    get_geostationary_bounding_box_in_lonlats,
 )
 from pyresample.gradient._gradient_search import (
     one_step_gradient_indices,
@@ -375,7 +375,7 @@ def _check_input_coordinates(dst_x, dst_y,
 def get_border_lonlats(geo_def):
     """Get the border x- and y-coordinates."""
     if geo_def.proj_dict['proj'] == 'geos':
-        lon_b, lat_b = get_geostationary_bounding_box(geo_def, 3600)
+        lon_b, lat_b = get_geostationary_bounding_box_in_lonlats(geo_def, 3600)
     else:
         lons, lats = geo_def.get_boundary_lonlats()
         lon_b = np.concatenate((lons.side1, lons.side2, lons.side3, lons.side4))
diff --git a/pyresample/test/test_geometry.py b/pyresample/test/test_geometry.py
index 791cd48..883b64e 100644
--- a/pyresample/test/test_geometry.py
+++ b/pyresample/test/test_geometry.py
@@ -2591,8 +2591,100 @@ def _get_fake_antimeridian_lonlats(use_dask: bool) -> tuple:
     return lons, lats
 
 
-class TestCrop(unittest.TestCase):
-    """Test the area helpers."""
+@pytest.fixture
+def truncated_geos_area():
+    """Create a truncated geostationary area."""
+    projection = {'a': '6378169', 'h': '35785831', 'lon_0': '9.5', 'no_defs': 'None', 'proj': 'geos',
+                  'rf': '295.488065897014', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
+    area_extent = (5567248.0742, 5570248.4773, -5570248.4773, 1393687.2705)
+    width = 3712
+    height = 1392
+    geos_area = geometry.AreaDefinition('msg_rss', "msg_rss", "msg_rss", projection, width, height, area_extent)
+    return geos_area
+
+
+@pytest.fixture
+def truncated_geos_area_in_space():
+    """Create a truncated geostationary area."""
+    projection = {'a': '6378169', 'h': '35785831', 'lon_0': '9.5', 'no_defs': 'None', 'proj': 'geos',
+                  'rf': '295.488065897014', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
+    area_extent = (5575000, 5575000, 5570000, 5570000)
+    width = 10
+    height = 10
+    geos_area = geometry.AreaDefinition('msg_rss', "msg_rss", "msg_rss", projection, width, height, area_extent)
+    return geos_area
+
+
+class TestGeostationaryTools:
+    """Test the geostationary bbox tools."""
+
+    def test_get_full_geostationary_bbox(self, truncated_geos_area):
+        nb_points = 20
+        x, y = geometry.get_full_geostationary_bounding_box_in_proj_coords(truncated_geos_area, nb_points)
+        assert len(x) == nb_points
+        assert len(y) == nb_points
+
+        assert x[0] != x[-1]
+        assert y[0] != y[-1]
+
+        expected_x = np.array([-5.43062255e+06, -5.16482897e+06, -4.39346593e+06, -3.19203985e+06,
+                               -1.67815466e+06, 3.32529726e-10, 1.67815466e+06, 3.19203985e+06,
+                               4.39346593e+06, 5.16482897e+06, 5.43062255e+06, 5.16482897e+06,
+                               4.39346593e+06, 3.19203985e+06, 1.67815466e+06, 3.32529726e-10,
+                               -1.67815466e+06, -3.19203985e+06, -4.39346593e+06, -5.16482897e+06])
+
+        expected_y = np.array([6.62789871e-10, 1.67242779e+06, 3.18114670e+06, 4.37847280e+06,
+                               5.14720348e+06, 5.41209002e+06, 5.14720348e+06, 4.37847280e+06,
+                               3.18114670e+06, 1.67242779e+06, -0.00000000e+00, -1.67242779e+06,
+                               -3.18114670e+06, -4.37847280e+06, -5.14720348e+06, -5.41209002e+06,
+                               -5.14720348e+06, -4.37847280e+06, -3.18114670e+06, -1.67242779e+06])
+
+        np.testing.assert_allclose(x, expected_x)
+        np.testing.assert_allclose(y, expected_y)
+
+    def test_get_geostationary_bbox_works_with_truncated_area(self, truncated_geos_area):
+        """Ensure the geostationary bbox works when truncated."""
+        lon, lat = geometry.get_geostationary_bounding_box_in_lonlats(truncated_geos_area, 20)
+
+        expected_lon = np.array(
+            [-64.24072434653284, -68.69662326361153, -65.92516214783112, -60.726360278290336,
+             -47.39851775032484, 9.500000000000018, 66.39851775032487, 79.72636027829033,
+             84.92516214783113, 87.69662326361151, 83.24072434653286])
+        expected_lat = np.array(
+            [14.554922655532085, 17.768795771961937, 35.34328897185421, 52.597860701318254, 69.00533141646078,
+             79.1481121862375, 69.00533141646076, 52.597860701318254, 35.34328897185421, 17.768795771961933,
+             14.554922655532085])
+        np.testing.assert_allclose(lon, expected_lon)
+        np.testing.assert_allclose(lat, expected_lat)
+
+    def test_get_geostationary_bbox_works_with_truncated_area_proj_coords(self, truncated_geos_area):
+        """Ensure the geostationary bbox works when truncated."""
+        x, y = geometry.get_geostationary_bounding_box_in_proj_coords(truncated_geos_area, 20)
+
+        expected_x = np.array(
+            [-5209128.302753595, -5164828.965702432, -4393465.934674804, -3192039.8468840676, -1678154.6586309497,
+             3.325297262895822e-10, 1678154.6586309501, 3192039.846884068, 4393465.934674805, 5164828.965702432,
+             5209128.302753594])
+        expected_y = np.array(
+            [1393687.2705, 1672427.7900638399, 3181146.6955466354, 4378472.798117005, 5147203.47659387,
+             5412090.016106332, 5147203.476593869, 4378472.798117005, 3181146.695546635, 1672427.7900638392,
+             1393687.2705])
+
+        np.testing.assert_allclose(x, expected_x)
+        np.testing.assert_allclose(y, expected_y)
+
+    def test_get_geostationary_bbox_does_not_contain_inf(self, truncated_geos_area):
+        """Ensure the geostationary bbox does not contain np.inf."""
+        lon, lat = geometry.get_geostationary_bounding_box_in_lonlats(truncated_geos_area, 20)
+        assert not any(np.isinf(lon))
+        assert not any(np.isinf(lat))
+
+    def test_get_geostationary_bbox_returns_empty_lonlats_in_space(self, truncated_geos_area_in_space):
+        """Ensure the geostationary bbox is empty when in space."""
+        lon, lat = geometry.get_geostationary_bounding_box_in_lonlats(truncated_geos_area_in_space, 20)
+
+        assert len(lon) == 0
+        assert len(lat) == 0
 
     def test_get_geostationary_bbox(self):
         """Get the geostationary bbox."""
@@ -2606,23 +2698,20 @@ class TestCrop(unittest.TestCase):
         geos_area.crs = CRS(proj_dict)
         geos_area.area_extent = [-5500000., -5500000., 5500000., 5500000.]
 
-        lon, lat = geometry.get_geostationary_bounding_box(geos_area, 20)
-        # This musk be equal to lon.
-        elon = np.array([-79.23372832, -78.19662326, -75.42516215, -70.22636028,
-                         -56.89851775, 0., 56.89851775, 70.22636028,
-                         75.42516215, 78.19662326, 79.23372832, 78.19662326,
-                         75.42516215, 70.22636028, 56.89851775, 0.,
-                        -56.89851775, -70.22636028, -75.42516215, -78.19662326,
-                        -79.23372832])
-        elat = np.array([0., 17.76879577, 35.34328897, 52.5978607,
-                         69.00533142, 79.14811219, 69.00533142, 52.5978607,
-                         35.34328897, 17.76879577, -0., -17.76879577,
-                         -35.34328897, -52.5978607, -69.00533142, -79.14811219,
-                         -69.00533142, -52.5978607, -35.34328897, -17.76879577,
-                         0.])
-
-        np.testing.assert_allclose(lon, elon, atol=1e-07)
-        np.testing.assert_allclose(lat, elat, atol=1e-07)
+        lon, lat = geometry.get_geostationary_bounding_box_in_lonlats(geos_area, 20)
+        expected_lon = np.array([-78.19662326, -75.42516215, -70.22636028,
+                                 -56.89851775, 0., 56.89851775, 70.22636028,
+                                 75.42516215, 78.19662326, 79.23372832, 78.19662326,
+                                 75.42516215, 70.22636028, 56.89851775, 0.,
+                                 -56.89851775, -70.22636028, -75.42516215, -78.19662326, -79.23372832, ])
+        expected_lat = np.array([17.76879577, 35.34328897, 52.5978607,
+                                 69.00533142, 79.14811219, 69.00533142, 52.5978607,
+                                 35.34328897, 17.76879577, -0., -17.76879577,
+                                 -35.34328897, -52.5978607, -69.00533142, -79.14811219,
+                                 -69.00533142, -52.5978607, -35.34328897, -17.76879577, 0.])
+
+        np.testing.assert_allclose(lon, expected_lon, atol=1e-07)
+        np.testing.assert_allclose(lat, expected_lat, atol=1e-07)
 
         geos_area = MagicMock()
         lon_0 = 10
@@ -2634,8 +2723,8 @@ class TestCrop(unittest.TestCase):
         geos_area.crs = CRS(proj_dict)
         geos_area.area_extent = [-5500000., -5500000., 5500000., 5500000.]
 
-        lon, lat = geometry.get_geostationary_bounding_box(geos_area, 20)
-        np.testing.assert_allclose(lon, elon + lon_0)
+        lon, lat = geometry.get_geostationary_bounding_box_in_lonlats(geos_area, 20)
+        np.testing.assert_allclose(lon, expected_lon + lon_0)
 
     def test_get_geostationary_angle_extent(self):
         """Get max geostationary angles."""
@@ -2675,6 +2764,10 @@ class TestCrop(unittest.TestCase):
         np.testing.assert_allclose(expected,
                                    geometry.get_geostationary_angle_extent(geos_area))
 
+
+class TestCrop(unittest.TestCase):
+    """Test the area helpers."""
+
     def test_sub_area(self):
         """Sub area slicing."""
         area = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD',
@@ -3156,20 +3249,19 @@ class TestBoundary(unittest.TestCase):
         # Check boundary vertices
         n_vertices = 10
         boundary = areadef.boundary(frequency=n_vertices, force_clockwise=False)
-        boundary.vertices.shape
 
         # Check boundary vertices is in correct order
-        expected_vertices = np.array([[-7.92337283e+01, 6.94302533e-15],
-                                      [-7.54251621e+01, 3.53432890e+01],
+        expected_vertices = np.array([[-7.54251621e+01, 3.53432890e+01],
                                       [-5.68985178e+01, 6.90053314e+01],
                                       [5.68985178e+01, 6.90053314e+01],
                                       [7.54251621e+01, 3.53432890e+01],
-                                      [7.92337283e+01, -6.94302533e-15],
+                                      [7.92337283e+01, -0.00000000e+00],
                                       [7.54251621e+01, -3.53432890e+01],
                                       [5.68985178e+01, -6.90053314e+01],
                                       [-5.68985178e+01, -6.90053314e+01],
-                                      [-7.54251621e+01, -3.53432890e+01]])
-        assert np.allclose(expected_vertices, boundary.vertices)
+                                      [-7.54251621e+01, -3.53432890e+01],
+                                      [-7.92337283e+01, 6.94302533e-15]])
+        np.testing.assert_allclose(expected_vertices, boundary.vertices)
 
     def test_global_platee_caree_projection(self):
         """Test boundary for global platee caree projection."""
diff --git a/pyresample/test/test_gradient.py b/pyresample/test/test_gradient.py
index c28a190..59d9b0a 100644
--- a/pyresample/test/test_gradient.py
+++ b/pyresample/test/test_gradient.py
@@ -599,7 +599,7 @@ def test_check_overlap():
     assert check_overlap(poly1, poly2) is False
 
 
-@mock.patch('pyresample.gradient.get_geostationary_bounding_box')
+@mock.patch('pyresample.gradient.get_geostationary_bounding_box_in_lonlats')
 def test_get_border_lonlats(get_geostationary_bounding_box):
     """Test that correct methods are called in get_border_lonlats()."""
     from pyresample.gradient import get_border_lonlats
diff --git a/setup.py b/setup.py
index 4dd9ba4..7a0bae9 100644
--- a/setup.py
+++ b/setup.py
@@ -29,13 +29,14 @@ import versioneer
 
 requirements = ['setuptools>=3.2', 'pyproj>=3.0', 'configobj',
                 'pykdtree>=1.3.1', 'pyyaml', 'numpy>=1.10.0',
+                "shapely",
                 ]
 
 if sys.version_info < (3, 10):
     requirements.append('importlib_metadata')
 
-test_requires = ['rasterio', 'dask', 'xarray', 'cartopy', 'pillow', 'matplotlib', 'scipy', 'zarr',
-                 'pytest-lazy-fixtures']
+test_requires = ['rasterio', 'dask', 'xarray', 'cartopy>=0.20.0', 'pillow', 'matplotlib', 'scipy', 'zarr',
+                 'pytest-lazy-fixtures', 'shapely']
 extras_require = {'numexpr': ['numexpr'],
                   'quicklook': ['matplotlib', 'cartopy>=0.20.0', 'pillow'],
                   'rasterio': ['rasterio'],
@@ -46,8 +47,6 @@ extras_require = {'numexpr': ['numexpr'],
                   'tests': test_requires}
 
 setup_requires = ['numpy>=1.10.0', 'cython']
-test_requires = ['rasterio', 'dask', 'xarray', 'cartopy>=0.20.0', 'pillow', 'matplotlib', 'scipy', 'zarr',
-                 'pytest-lazy-fixture']
 
 if sys.platform.startswith("win"):
     extra_compile_args = []
['pyresample/gradient/__init__.py', 'setup.py', 'pyresample/test/test_gradient.py', 'pyresample/test/test_geometry.py', 'pyresample/geometry.py', '.github/workflows/deploy.yaml']