diff --git a/examples/multi-instrument.py b/examples/multi-instrument.py index e9d2f8b..a231de2 100644 --- a/examples/multi-instrument.py +++ b/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() diff --git a/synthesizAR/instruments/base.py b/synthesizAR/instruments/base.py index 63ad4ee..fa209de 100644 --- a/synthesizAR/instruments/base.py +++ b/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 diff --git a/synthesizAR/util/util.py b/synthesizAR/util/util.py index 9be99b5..81f56c7 100644 --- a/synthesizAR/util/util.py +++ b/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 `_ diff --git a/synthesizAR/visualize/fieldlines.py b/synthesizAR/visualize/fieldlines.py index 59c4370..9a128f2 100644 --- a/synthesizAR/visualize/fieldlines.py +++ b/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