Fix spill over of ocean/and features in cartopy plots in case of geostationary full disc plot.

pyresample/utils/cartopy.py

@@ -73,65 +73,71 @@ def __init__(self,
 
         # For geostationary full disc the boundary needs to be the actuall boundary of the earth disc otherwise
         # when ocean or land features are added to the cartopy plot these spill over.
-        if "Geostationary Satellite" in crs.to_wkt():
-            # This code is copied over from the 'Geostationary' class in  'cartopy/lib/cartopy/crs.py'.
-            satellite_height = 35785831
-            false_easting = 0
-            false_northing = 0
-            a = float(ccrs.WGS84_SEMIMAJOR_AXIS)
-            b = float(ccrs.WGS84_SEMIMINOR_AXIS)
-            h = float(satellite_height)
-
-            # To find the bound we trace around where the line from the satellite
-            # is tangent to the surface. This involves trigonometry on a sphere
-            # centered at the satellite. The two scanning angles form two legs of
-            # triangle on this sphere--the hypotenuse "c" (angle arc) is controlled
-            # by distance from center to the edge of the ellipse being seen.
-
-            # This is one of the angles in the spherical triangle and used to
-            # rotate around and "scan" the boundary
-            angleA = np.linspace(0, -2 * np.pi, 91)  # Clockwise boundary.
-
-            # Convert the angle around center to the proper value to use in the
-            # parametric form of an ellipse
-            th = np.arctan(a / b * np.tan(angleA))
-
-            # Given the position on the ellipse, what is the distance from center
-            # to the ellipse--and thus the tangent point
-            r = np.hypot(a * np.cos(th), b * np.sin(th))
-            sat_dist = a + h
-
-            # Using this distance, solve for sin and tan of c in the triangle that
-            # includes the satellite, Earth center, and tangent point--we need to
-            # figure out the location of this tangent point on the elliptical
-            # cross-section through the Earth towards the satellite, where the
-            # major axis is a and the minor is r. With the ellipse centered on the
-            # Earth and the satellite on the y-axis (at y = a + h = sat_dist), the
-            # equation for an ellipse and some calculus gives us the tangent point
-            # (x0, y0) as:
-            # y0 = a**2 / sat_dist
-            # x0 = r * np.sqrt(1 - a**2 / sat_dist**2)
-            # which gives:
-            # sin_c = x0 / np.hypot(x0, sat_dist - y0)
-            # tan_c = x0 / (sat_dist - y0)
-            # A bit of algebra combines these to give directly:
-            sin_c = r / np.sqrt(sat_dist ** 2 - a ** 2 + r ** 2)
-            tan_c = r / np.sqrt(sat_dist ** 2 - a ** 2)
-
-            # Using Napier's rules for right spherical triangles R2 and R6,
-            # (See https://en.wikipedia.org/wiki/Spherical_trigonometry), we can
-            # solve for arc angles b and a, which are our x and y scanning angles,
-            # respectively.
-            coords = np.vstack([np.arctan(np.cos(angleA) * tan_c),  # R6
-                                np.arcsin(np.sin(angleA) * sin_c)])  # R2
-
-            # Need to multiply scanning angles by satellite height to get to the
-            # actual native coordinates for the projection.
-            coords *= h
-            coords += np.array([[false_easting], [false_northing]])
-            self._boundary = sgeom.LinearRing(coords.T)
-        else:
+        if "Geostationary Satellite" not in crs.to_wkt():
             self._boundary = super().boundary
+        else:
+            self._boundary = self._geos_boundary()
+
+    def _geos_boundary(self):
+        """Calculate full disk boundary.
+
+        This code is copied over from the 'Geostationary' class in  'cartopy/lib/cartopy/crs.py'.
+        """
+        satellite_height = 35785831
+        false_easting = 0
+        false_northing = 0
+        a = float(ccrs.WGS84_SEMIMAJOR_AXIS)
+        b = float(ccrs.WGS84_SEMIMINOR_AXIS)
+        h = float(satellite_height)
+
+        # To find the bound we trace around where the line from the satellite
+        # is tangent to the surface. This involves trigonometry on a sphere
+        # centered at the satellite. The two scanning angles form two legs of
+        # triangle on this sphere--the hypotenuse "c" (angle arc) is controlled
+        # by distance from center to the edge of the ellipse being seen.
+
+        # This is one of the angles in the spherical triangle and used to
+        # rotate around and "scan" the boundary
+        angleA = np.linspace(0, -2 * np.pi, 91)  # Clockwise boundary.
+
+        # Convert the angle around center to the proper value to use in the
+        # parametric form of an ellipse
+        th = np.arctan(a / b * np.tan(angleA))
+
+        # Given the position on the ellipse, what is the distance from center
+        # to the ellipse--and thus the tangent point
+        r = np.hypot(a * np.cos(th), b * np.sin(th))
+        sat_dist = a + h
+
+        # Using this distance, solve for sin and tan of c in the triangle that
+        # includes the satellite, Earth center, and tangent point--we need to
+        # figure out the location of this tangent point on the elliptical
+        # cross-section through the Earth towards the satellite, where the
+        # major axis is a and the minor is r. With the ellipse centered on the
+        # Earth and the satellite on the y-axis (at y = a + h = sat_dist), the
+        # equation for an ellipse and some calculus gives us the tangent point
+        # (x0, y0) as:
+        # y0 = a**2 / sat_dist
+        # x0 = r * np.sqrt(1 - a**2 / sat_dist**2)
+        # which gives:
+        # sin_c = x0 / np.hypot(x0, sat_dist - y0)
+        # tan_c = x0 / (sat_dist - y0)
+        # A bit of algebra combines these to give directly:
+        sin_c = r / np.sqrt(sat_dist ** 2 - a ** 2 + r ** 2)
+        tan_c = r / np.sqrt(sat_dist ** 2 - a ** 2)
+
+        # Using Napier's rules for right spherical triangles R2 and R6,
+        # (See https://en.wikipedia.org/wiki/Spherical_trigonometry), we can
+        # solve for arc angles b and a, which are our x and y scanning angles,
+        # respectively.
+        coords = np.vstack([np.arctan(np.cos(angleA) * tan_c),  # R6
+                            np.arcsin(np.sin(angleA) * sin_c)])  # R2
+
+        # Need to multiply scanning angles by satellite height to get to the
+        # actual native coordinates for the projection.
+        coords *= h
+        coords += np.array([[false_easting], [false_northing]])
+        return sgeom.LinearRing(coords.T)
 
     @property
     def boundary(self):