Use solar_angle_equivalency and Helioprojective.is_visible() (#100)

* use solar angle equivalency

* use sunpy is_visible method

* small example fix

Author: wtbarnes

examples/multi-instrument.py

@@ -136,8 +136,7 @@ def get_instrument_name(self, *args):
 ###############################################################################
 # We can then use our custom instrument class in the exact same way as our
 # predefined classes to model the emission from EUVI. Note that we'll only do
-# this for the 171 $\AA$ channel.
+# this for the 171 Å channel.
 euvi = InstrumentSTEREOEUVI([0, 1]*u.s, stereo_a, fov_center=center, fov_width=(250, 250)*u.arcsec)
 euvi_images = euvi.observe(skeleton, channels=euvi.channels[2:3])
-for k in euvi_images:
-    euvi_images[k][0].peek()
+euvi_images['171'][0].peek()

synthesizAR/instruments/base.py

@@ -12,10 +12,11 @@
 from dataclasses import dataclass
 from scipy.interpolate import interp1d
 from scipy.ndimage import gaussian_filter
-from sunpy.coordinates.frames import HeliographicStonyhurst, Helioprojective
+from sunpy.coordinates import HeliographicStonyhurst, Helioprojective
+from sunpy.coordinates.utils import solar_angle_equivalency
 from sunpy.map import make_fitswcs_header, Map
 
-from synthesizAR.util import find_minimum_fov, is_visible
+from synthesizAR.util import find_minimum_fov
 from synthesizAR.util.decorators import return_quantity_as_tuple
 
 __all__ = ['ChannelBase', 'InstrumentBase']
@@ -129,8 +130,9 @@ def pixel_area(self) -> u.cm**2:
         """
         Pixel area
         """
-        w_x, w_y = (1*u.pix * self.resolution).to(u.radian).value * self.observer.radius
-        return w_x * w_y
+        sa_equiv = solar_angle_equivalency(self.observer)
+        res = (1*u.pix * self.resolution).to('cm', equivalencies=sa_equiv)
+        return res[0] * res[1]
 
     def convolve_with_psf(self, smap, channel):
         """
@@ -345,7 +347,7 @@ def integrate_los(self, time, channel, skeleton, coordinates_centers, bins, bin_
         if not self.average_over_los:
             kernels *= (skeleton.all_cross_sectional_areas / self.pixel_area).decompose() * skeleton.all_widths
         if check_visible:
-            visible = is_visible(coordinates_centers, self.observer)
+            visible = coordinates_centers.is_visible()
         else:
             visible = np.ones(kernels.shape)
         # Bin

synthesizAR/util/util.py

@@ -6,7 +6,6 @@
 import astropy.units as u
 import numpy as np
 import sunpy.coordinates
-import sunpy.sun.constants as sun_const
 import warnings
 
 from astropy.coordinates import SkyCoord
@@ -20,7 +19,6 @@
     'los_velocity',
     'coord_in_fov',
     'find_minimum_fov',
-    'is_visible',
     'from_pfsspack',
     'from_pfsspy',
     'change_obstime',
@@ -98,28 +96,6 @@ def find_minimum_fov(coordinates, padding=None):
     return bottom_left_corner, top_right_corner
 
 
-def is_visible(coords, observer):
-    """
-    Create mask of coordinates not blocked by the solar disk.
-
-    Parameters
-    ----------
-    coords : `~astropy.coordinates.SkyCoord`
-        Helioprojective oordinates of the object(s) of interest
-    observer : `~astropy.coordinates.SkyCoord`
-        Heliographic-Stonyhurst Location of the observer
-    """
-    theta_x = coords.Tx
-    theta_y = coords.Ty
-    distance = coords.distance
-    rsun_obs = ((sun_const.radius / (observer.radius - sun_const.radius)).decompose()
-                * u.radian).to(u.arcsec)
-    off_disk = np.sqrt(theta_x**2 + theta_y**2) > rsun_obs
-    in_front_of_disk = distance - observer.radius < 0.
-
-    return np.any(np.stack([off_disk, in_front_of_disk], axis=1), axis=1)
-
-
 def from_pfsspack(pfss_fieldlines):
     """
     Convert fieldline coordinates output from the SSW package `pfss <http://www.lmsal.com/~derosa/pfsspack/>`_

synthesizAR/visualize/fieldlines.py

@@ -1,18 +1,18 @@
 """
 Visualizaition functions related to 1D fieldlines
 """
-import numpy as np
-import matplotlib.pyplot as plt
-from matplotlib.colors import Normalize
 import astropy.units as u
-from astropy.time import Time
+import matplotlib.pyplot as plt
+import numpy as np
+
 from astropy.coordinates import SkyCoord
+from astropy.time import Time
 from astropy.visualization import ImageNormalize
-from sunpy.map import GenericMap, make_fitswcs_header
 from sunpy.coordinates import Helioprojective
 from sunpy.coordinates.ephemeris import get_earth
+from sunpy.map import GenericMap, make_fitswcs_header
 
-from synthesizAR.util import is_visible, find_minimum_fov
+from synthesizAR.util import find_minimum_fov
 
 __all__ = ['set_ax_lims', 'plot_fieldlines']
 
@@ -55,20 +55,13 @@ def plot_fieldlines(*coords,
         If True, draw the HGS grid
     axes_limits : `tuple`, optional
         Tuple of world coordinates (axis1, axis2)
-
-    Other Parameters
-    ----------------
-    plot_kwargs : `dict`
+    plot_kwargs : `dict`, optional
         Additional parameters to pass to `~matplotlib.pyplot.plot` when
         drawing field lines.
-    grid_kwargs : `dict`
+    grid_kwargs : `dict`, optional
         Additional parameters to pass to `~sunpy.map.Map.draw_grid`
-    imshow_kwargs : `dict`
+    imshow_kwargs : `dict`, optional
         Additional parameters to pass to `~sunpy.map.Map.plot`
-
-    See Also
-    --------
-    synthesizAR.util.is_visible
     """
     plot_kwargs = {'color': 'k', 'lw': 1}
     grid_kwargs = {'grid_spacing': 10*u.deg, 'color': 'k', 'alpha': 0.75}
@@ -100,7 +93,7 @@ def plot_fieldlines(*coords,
     for coord in coords:
         c = coord.transform_to(image_map.coordinate_frame)
         if check_visible:
-            c = c[is_visible(c, image_map.observer_coordinate)]
+            c = c[c.is_visible()]
         transformed_coords.append(c)
         if len(c) == 0:
             continue  # Matplotlib throws exception when no points are visible