From dad40548da5bd427895dac16a80560f1239a38bb Mon Sep 17 00:00:00 2001 From: BENR0 Date: Wed, 17 Jul 2024 12:39:41 +0200 Subject: [PATCH] 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):