diff --git a/changelog/7118.feature.rst b/changelog/7118.feature.rst new file mode 100644 index 00000000000..40d41a45aef --- /dev/null +++ b/changelog/7118.feature.rst @@ -0,0 +1 @@ +Added a :meth:`sunpy.coordinates.Helioprojective.is_visible` method to return whether the coordinate is visible (i.e., not obscured from the observer assuming that the Sun is an opaque sphere). diff --git a/sunpy/coordinates/frames.py b/sunpy/coordinates/frames.py index 4606188946f..04143b0ffbf 100644 --- a/sunpy/coordinates/frames.py +++ b/sunpy/coordinates/frames.py @@ -609,6 +609,85 @@ def make_3d(self): lat=lat, distance=d)) + @u.quantity_input + def is_visible(self, *, tolerance: u.m = 1*u.m): + """ + Returns whether the coordinate is on the visible side of the Sun. + + A coordinate is visible if it can been seen from the observer (the ``observer`` + frame attribute) assuming that the Sun is an opaque sphere with a fixed radius + (the ``rsun`` frame attribute). The visible side of the Sun is always smaller + than a full hemisphere. + + Parameters + ---------- + tolerance : `~astropy.units.Quantity` + The depth below the surface of the Sun that should be treated as + transparent. + + Notes + ----- + If the coordinate is 2D, it is automatically deemed visible. A 2D coordinate + describes a look direction from the observer, who would simply see whatever is + in "front", and thus cannot correspond to a point hidden from the observer. + + The ``tolerance`` parameter accommodates situations where the limitations of + numerical precision would falsely conclude that a coordinate is not visible. + For example, a coordinate that is expressly created to be on the solar surface + may be calculated to be slightly below the surface, and hence not visible if + there is no tolerance. However, a consequence of the ``tolerance`` parameter + is that a coordinate that is formally on the far side of the Sun but is + extremely close to the solar limb can be evaluated as visible. With the + default ``tolerance`` value of 1 meter, a coordinate on the surface of the Sun + can be up to 11 arcseconds of heliographic longitude past the solar limb and + still be evaluated as visible. + + Examples + -------- + >>> import numpy as np + >>> import astropy.units as u + >>> from astropy.coordinates import SkyCoord + >>> from sunpy.coordinates import HeliographicStonyhurst, Helioprojective + + >>> hpc_frame = Helioprojective(observer='earth', obstime='2023-08-03') + + >>> in_front = SkyCoord(0*u.arcsec, 0*u.arcsec, 0.5*u.AU, frame=hpc_frame) + >>> print(in_front.is_visible()) + True + + >>> behind = SkyCoord(0*u.arcsec, 0*u.arcsec, 1.5*u.AU, frame=hpc_frame) + >>> print(behind.is_visible()) + False + + >>> hgs_array = SkyCoord(np.arange(-180, 180, 60)*u.deg, [0]*6*u.deg, + ... frame='heliographic_stonyhurst', obstime='2023-08-03') + >>> print(hgs_array) + + >>> print(hgs_array.transform_to(hpc_frame).is_visible()) + [False False True True True False] + """ + # If the coordinate is 2D, it must be visible + if self._is_2d: + return np.ones_like(self.data, dtype=bool) + + # Use a slightly smaller solar radius to accommodate numerical precision + solar_radius = self.rsun - tolerance + + data = self.cartesian + data_to_sun = self.observer.radius * CartesianRepresentation(1, 0, 0) - data + + # When representing the helioprojective point as true Cartesian, the X value is the + # distance from the observer to the point in the sunward direction + is_behind_observer = data.x < 0 + # When comparing heliocentric angles, we compare the sine values and hence avoid calling arcsin() + is_beyond_limb = np.sqrt(data.y **2 + data.z **2) / data.norm() > solar_radius / self.observer.radius + is_above_surface = data_to_sun.norm() >= solar_radius + is_on_near_side = data.dot(data_to_sun) >= 0 + + return is_behind_observer | is_beyond_limb | (is_on_near_side & is_above_surface) + _spherical_screen = None @classmethod diff --git a/sunpy/coordinates/tests/test_frames.py b/sunpy/coordinates/tests/test_frames.py index b8db7236356..fe90fe708f5 100644 --- a/sunpy/coordinates/tests/test_frames.py +++ b/sunpy/coordinates/tests/test_frames.py @@ -242,6 +242,29 @@ def test_hpc_obstime_from_observer(): assert hpc.obstime is None +def test_hpc_is_visible_2d(): + hpc = Helioprojective(2000*u.arcsec, 2000*u.arcsec, + observer='earth', obstime='2023-08-03') + assert hpc.is_visible() + + +def test_hpc_is_visible(): + hpc = Helioprojective([0]*2*u.arcsec, [0]*2*u.arcsec, [0.5, 1.5]*u.AU, + observer='earth', obstime='2023-08-03') + assert (hpc.is_visible() == [True, False]).all() + + +def test_hpc_is_visible_tolerance(): + hpc = Helioprojective(200*u.arcsec, 0*u.arcsec, + observer='earth', obstime='2023-08-03').make_3d() + + # Due to the limitations of numerical precision, the coordinate will be computed to be slightly + # below the solar surface, and thus invisible when the tolerance is set to zero + assert not hpc.is_visible(tolerance=0*u.m) + + assert hpc.is_visible(tolerance=1*u.m) + + # ============================================================================== # ## Heliographic Tests # ==============================================================================