Merge pull request #601 from BENR0/feat_html_repr

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

pyresample/_formatting_html.py

@@ -25,6 +25,7 @@
 import numpy as np
 
 import pyresample.geometry as geom
+from pyresample.utils.proj4 import ignore_pyproj_proj_warnings
 
 try:
     import cartopy
@@ -65,15 +66,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 +99,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)
 
@@ -208,7 +194,8 @@ def proj_area_attrs_section(area: 'geom.AreaDefinition') -> str: # noqa F821
 
     """
     resolution_str = "/".join([str(round(x, 1)) for x in area.resolution])
-    proj_dict = area.proj_dict
+    with ignore_pyproj_proj_warnings():
+        proj_dict = area.proj_dict
     proj_str = "{{{}}}".format(", ".join(["'%s': '%s'" % (str(k), str(proj_dict[k])) for k in
                                           sorted(proj_dict.keys())]))
     area_units = proj_dict.get("units", "")
@@ -222,7 +209,7 @@ def proj_area_attrs_section(area: 'geom.AreaDefinition') -> str: # noqa F821
                   f"<dt>Width/Height</dt><dd>{area.width}/{area.height} Pixel</dd>"
                   f"<dt>Resolution x/y (SSP)</dt><dd>{resolution_str} {area_units}</dd>"
                   f"<dt>Extent (ll_x, ll_y, ur_x, ur_y)</dt>"
-                  f"<dd>{tuple(round(x, 4) for x in area.area_extent)}</dd>"
+                  f"<dd>{tuple(round(float(x), 4) for x in area.area_extent)}</dd>"
                   "</dl>"
                   )
 
@@ -319,14 +306,15 @@ def swath_area_attrs_section(area: 'geom.SwathDefinition') -> str: # noqa F821
     return f"{coll}"
 
 
-def area_repr(area: Union['geom.AreaDefinition', 'geom.SwathDefinition'], include_header: Optional[bool] = True, # noqa F821
-              include_static_files: Optional[bool] = True):
+def area_repr(area: Union['geom.AreaDefinition', 'geom.SwathDefinition'],
+              include_header: bool = True,
+              include_static_files: bool = True):
     """Return html repr of an AreaDefinition.
 
     Args:
         area : Area definition.
         include_header : If true a header with object type will be included in
-            the html. This is mainly intented for display in Jupyter Notebooks. For the
+            the html. This is mainly intended for display in Jupyter Notebooks. For the
             display in the overview of area definitions for the Satpy documentation this
             should be set to false.
         include_static_files : Load and include css and html needed for representation.

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,77 @@ 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" not in crs.to_wkt():
+            self._boundary = super().boundary
+        else:
+            self._boundary = self._geos_boundary()
+
+    @staticmethod
+    def _geos_boundary():
+        """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):
+        """Return boundary."""
+        return self._boundary