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

pyresample/_formatting_html.py

@@ -65,15 +65,11 @@ def _icon(icon_name):
 
 
 def plot_area_def(area: Union['geom.AreaDefinition', 'geom.SwathDefinition'], # noqa F821
-                  feature_res: Optional[str] = "110m",
                   fmt: Optional[Literal["svg", "png", None]] = None) -> Union[str, None]:
     """Plot area.
 
     Args:
         area : Area/Swath to plot.
-        feature_res :
-            Resolution of the features added to the map. Argument is handed over
-            to `scale` parameter in cartopy.feature.
         fmt : Output format of the plot. The output is the string representation of
             the respective format xml for svg and base64 for png. Either svg or png.
             If None (default) plot is just shown.
@@ -102,21 +98,10 @@ def plot_area_def(area: Union['geom.AreaDefinition', 'geom.SwathDefinition'], #
         bounds = poly.buffer(5).bounds
         ax.set_extent([bounds[0], bounds[2], bounds[1], bounds[3]], crs=cartopy.crs.CRS(area.crs))
 
-    coastlines = cartopy.feature.NaturalEarthFeature(category="physical",
-                                                     name="coastline",
-                                                     scale=feature_res,
-                                                     linewidth=1,
-                                                     facecolor="never")
-    borders = cartopy.feature.NaturalEarthFeature(category="cultural",
-                                                  name="admin_0_boundary_lines_land", # noqa E114
-                                                  scale=feature_res,
-                                                  edgecolor="black",
-                                                  facecolor="never") # noqa E1>
-    ocean = cartopy.feature.OCEAN
-
-    ax.add_feature(borders)
-    ax.add_feature(coastlines)
-    ax.add_feature(ocean, color="lightgrey")
+    ax.add_feature(cartopy.feature.OCEAN)
+    ax.add_feature(cartopy.feature.LAND)
+    ax.add_feature(cartopy.feature.COASTLINE)
+    ax.add_feature(cartopy.feature.BORDERS)
 
     plt.tight_layout(pad=0)
 

pyresample/utils/cartopy.py

@@ -31,6 +31,7 @@
     DataArray = np.ndarray
 
 import cartopy.crs as ccrs
+import shapely.geometry as sgeom
 
 
 class Projection(ccrs.Projection):
@@ -69,3 +70,70 @@ def __init__(self,
         if self.bounds is not None:
             x0, x1, y0, y1 = self.bounds
             self.threshold = min(abs(x1 - x0), abs(y1 - y0)) / 100.
+
+        # 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:
+            self._boundary = super().boundary
+
+    @property
+    def boundary(self):
+        """Return boundary."""
+        return self._boundary