From f8d483a39ba4e99086a5f833e4b68bf9482e96c1 Mon Sep 17 00:00:00 2001 From: BENR0 Date: Fri, 14 Jun 2024 10:57:44 +0200 Subject: [PATCH 1/6] remove custom feature setup in favour of automatic cartopy features --- pyresample/_formatting_html.py | 23 ++++------------------- 1 file changed, 4 insertions(+), 19 deletions(-) diff --git a/pyresample/_formatting_html.py b/pyresample/_formatting_html.py index d1a087fe..c1fefcc9 100644 --- a/pyresample/_formatting_html.py +++ b/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) From d13528d4d69ac130eb55313118c08dfa4a96a5e4 Mon Sep 17 00:00:00 2001 From: BENR0 Date: Fri, 14 Jun 2024 11:01:20 +0200 Subject: [PATCH 2/6] fix: spill over of cartopy features in case of full disc --- pyresample/utils/cartopy.py | 71 +++++++++++++++++++++++++++++++++++++ 1 file changed, 71 insertions(+) diff --git a/pyresample/utils/cartopy.py b/pyresample/utils/cartopy.py index d298241e..128a8c6c 100644 --- a/pyresample/utils/cartopy.py +++ b/pyresample/utils/cartopy.py @@ -31,6 +31,10 @@ DataArray = np.ndarray import cartopy.crs as ccrs +import shapely.geometry as sgeom + +WGS84_SEMIMAJOR_AXIS = 6378137.0 +WGS84_SEMIMINOR_AXIS = 6356752.3142 class Projection(ccrs.Projection): @@ -69,3 +73,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(WGS84_SEMIMAJOR_AXIS) + b = float(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 From e16203a4e5877d3ad196d5e23bb15e4f382d08c6 Mon Sep 17 00:00:00 2001 From: BENR0 Date: Fri, 14 Jun 2024 11:16:56 +0200 Subject: [PATCH 3/6] refactor: import WGS parameters from cartopy --- pyresample/utils/cartopy.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/pyresample/utils/cartopy.py b/pyresample/utils/cartopy.py index 128a8c6c..29dea4e4 100644 --- a/pyresample/utils/cartopy.py +++ b/pyresample/utils/cartopy.py @@ -33,9 +33,6 @@ import cartopy.crs as ccrs import shapely.geometry as sgeom -WGS84_SEMIMAJOR_AXIS = 6378137.0 -WGS84_SEMIMINOR_AXIS = 6356752.3142 - class Projection(ccrs.Projection): """Flexible generic Projection with optionally specified bounds.""" @@ -81,8 +78,8 @@ def __init__(self, satellite_height = 35785831 false_easting = 0 false_northing = 0 - a = float(WGS84_SEMIMAJOR_AXIS) - b = float(WGS84_SEMIMINOR_AXIS) + 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 From dad40548da5bd427895dac16a80560f1239a38bb Mon Sep 17 00:00:00 2001 From: BENR0 Date: Wed, 17 Jul 2024 12:39:41 +0200 Subject: [PATCH 4/6] refactor projection class --- pyresample/utils/cartopy.py | 122 +++++++++++++++++++----------------- 1 file changed, 64 insertions(+), 58 deletions(-) diff --git a/pyresample/utils/cartopy.py b/pyresample/utils/cartopy.py index 29dea4e4..77ebe9d5 100644 --- a/pyresample/utils/cartopy.py +++ b/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): From f31cc528fd43cb4dba8d5aad5341e04853ba0ba9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Benjamin=20R=C3=B6sner?= Date: Thu, 25 Jul 2024 13:12:48 +0200 Subject: [PATCH 5/6] refactor: make geos boundary static method --- pyresample/utils/cartopy.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyresample/utils/cartopy.py b/pyresample/utils/cartopy.py index 77ebe9d5..5e29d779 100644 --- a/pyresample/utils/cartopy.py +++ b/pyresample/utils/cartopy.py @@ -78,7 +78,8 @@ def __init__(self, else: self._boundary = self._geos_boundary() - def _geos_boundary(self): + @staticmethod + def _geos_boundary(): """Calculate full disk boundary. This code is copied over from the 'Geostationary' class in 'cartopy/lib/cartopy/crs.py'. From 9bf142ca211997ee63815fd3dce61766a9ee4de8 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Tue, 30 Jul 2024 14:29:14 -0500 Subject: [PATCH 6/6] Cleanup HTML repr proj warnings and extent floats --- pyresample/_formatting_html.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/pyresample/_formatting_html.py b/pyresample/_formatting_html.py index c1fefcc9..6af39adf 100644 --- a/pyresample/_formatting_html.py +++ b/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 @@ -193,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", "") @@ -207,7 +209,7 @@ def proj_area_attrs_section(area: 'geom.AreaDefinition') -> str: # noqa F821 f"
Width/Height
{area.width}/{area.height} Pixel
" f"
Resolution x/y (SSP)
{resolution_str} {area_units}
" f"
Extent (ll_x, ll_y, ur_x, ur_y)
" - f"
{tuple(round(x, 4) for x in area.area_extent)}
" + f"
{tuple(round(float(x), 4) for x in area.area_extent)}
" "" ) @@ -299,14 +301,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.