Source code for sunpy.visualization.drawing

"""
This module provides functions that draw on Astropy's WCSAxes.
"""
import numpy as np
from matplotlib import patches
from matplotlib.path import Path

import astropy.units as u
from astropy.constants import R_sun
from astropy.coordinates import SkyCoord
from astropy.visualization.wcsaxes.wcsapi import wcsapi_to_celestial_frame

from sunpy.coordinates import HeliographicCarrington, HeliographicStonyhurst
from sunpy.coordinates.frames import HeliocentricInertial, Helioprojective
from sunpy.coordinates.sun import _angular_radius
from sunpy.coordinates.utils import get_limb_coordinates
from sunpy.visualization import wcsaxes_compat

__all__ = ["limb", "equator", "prime_meridian"]


[docs]@u.quantity_input def limb(axes, observer, *, rsun: u.m = R_sun, resolution=1000, **kwargs): """ Draws the solar limb as seen by the specified observer. The limb is a circle for only the simplest plots. If the specified observer of the limb is different from the observer of the coordinate frame of the plot axes, not only may the limb not be a true circle, a portion of the limb may be hidden from the observer. In that case, the circle is divided into visible and hidden segments, represented by solid and dotted lines, respectively. Parameters ---------- axes : `~matplotlib.axes` or ``None`` Axes to plot limb on. observer : `astropy.coordinates.SkyCoord` Observer coordinate for which the limb is drawn. rsun : `~astropy.units.Quantity` Solar radius (in physical length units) at which to draw the limb. Defaults to the standard photospheric radius. resolution : `int` The number of points to use to represent the limb. Returns ------- visible : `~matplotlib.patches.Polygon` or `~matplotlib.patches.Circle` or None The patch added to the axes for the visible part of the limb (i.e., the "near" side of the Sun). hidden : `~matplotlib.patches.Polygon` or None The patch added to the axes for the hidden part of the limb (i.e., the "far" side of the Sun). Notes ----- Keyword arguments are passed onto the patches. If the limb is a true circle, ``visible`` will instead be `~matplotlib.patches.Circle` and ``hidden`` will be ``None``. If there are no hidden points (e.g., on a synoptic map any limb is fully visible) ``hidden`` will be ``None``. If there are no visible points (e.g., for an observer on the opposite side of the Sun to the map observer) ``visible`` will be ``None``. To avoid triggering Matplotlib auto-scaling, these patches are added as artists instead of patches. One consequence is that the plot legend is not populated automatically when the limb is specified with a text label. See :ref:`sphx_glr_gallery_text_labels_and_annotations_custom_legends.py` in the Matplotlib documentation for examples of creating a custom legend. """ if not wcsaxes_compat.is_wcsaxes(axes): raise ValueError('axes must be a WCSAxes') c_kw = {'fill': False, 'color': 'white', 'zorder': 100} c_kw.update(kwargs) transform = axes.get_transform('world') # transform is always passed on as a keyword argument c_kw.setdefault('transform', transform) # If the observer matches the axes's observer frame and is Helioprojective, use a Circle axes_frame = wcsapi_to_celestial_frame(axes.wcs) if isinstance(axes_frame, Helioprojective): axes_observer = SkyCoord(axes_frame.observer) if axes_observer.separation_3d(observer) < 1 * u.m: distance = observer.transform_to(HeliocentricInertial).spherical.distance # Obtain the solar radius and the world->pixel transform angular_radius = _angular_radius(rsun, distance) circ = patches.Circle([0, 0], radius=angular_radius.to_value(u.deg), **c_kw) axes.add_artist(circ) return circ, None # Otherwise, we use Polygon to be able to distort the limb # Create the limb coordinate array using Heliocentric Radial limb = get_limb_coordinates(observer, rsun, resolution) # Transform the limb to the axes frame and get the 2D vertices visible, hidden = _plot_vertices(limb, axes, axes_frame, rsun, **kwargs) return visible, hidden
[docs]@u.quantity_input def equator(axes, *, rsun: u.m = R_sun, resolution=500, **kwargs): """ Draws the solar equator as seen by the axes observer. Hidden parts are drawn as a dotted line. Parameters ---------- axes : `matplotlib.axes.Axes` The axes to plot the equator on. rsun : `~astropy.units.Quantity` Solar radius (in physical length units) at which to draw the solar equator. Defaults to the standard photospheric radius. resolution : `int` The number of points used to represent the equator. Returns ------- visible : `~matplotlib.patches.Polygon` The patch added to the axes for the visible part of the solar equator. hidden : `~matplotlib.patches.Polygon` The patch added to the axes for the hidden part of the solar equator. Examples -------- >>> import sunpy.map >>> from sunpy.visualization import drawing >>> import sunpy.data.sample # doctest: +REMOTE_DATA >>> aia = sunpy.map.Map(sunpy.data.sample.AIA_171_IMAGE) # doctest: +REMOTE_DATA >>> fig = plt.figure() # doctest: +SKIP >>> ax = fig.add_subplot(projection=aia) # doctest: +SKIP >>> aia.plot() # doctest: +SKIP >>> drawing.equator(ax) # doctest: +SKIP """ if not wcsaxes_compat.is_wcsaxes(axes): raise ValueError('axes must be a WCSAxes') axes_frame = wcsapi_to_celestial_frame(axes.wcs) lat = 0*u.deg lat0 = SkyCoord(np.linspace(-180, 179, resolution) * u.deg, np.ones(resolution) * lat, radius=rsun, frame=HeliographicStonyhurst(obstime=axes_frame.obstime)) visible, hidden = _plot_vertices(lat0, axes, axes_frame, rsun, **kwargs) return visible, hidden
[docs]@u.quantity_input def prime_meridian(axes, *, rsun: u.m = R_sun, resolution=500, **kwargs): """ Draws the solar prime meridian (zero Carrington longitude) as seen by the axes observer. Hidden parts are drawn as a dotted line. Parameters ---------- axes : `matplotlib.axes.Axes` The axes to plot the prime meridian on, or "None" to use current axes. rsun : `~astropy.units.Quantity` Solar radius (in physical length units) at which to draw the solar prime meridian. Defaults to the standard photospheric radius. resolution : `int` The number of points used to represent the prime meridian. Returns ------- visible : `~matplotlib.patches.Polygon` The patch added to the axes for the visible part of the solar equator. hidden : `~matplotlib.patches.Polygon` The patch added to the axes for the hidden part of the solar equator. """ if not wcsaxes_compat.is_wcsaxes(axes): raise ValueError('axes must be a WCSAxes') axes_frame = wcsapi_to_celestial_frame(axes.wcs) if not hasattr(axes_frame, 'observer'): raise ValueError('the coordinate frame of the WCSAxes does not have an observer, ' 'so zero Carrington longitude cannot be determined.') observer = axes_frame.observer lon = 0*u.deg lon0 = SkyCoord(np.ones(resolution) * lon, np.linspace(-90, 90, resolution) * u.deg, radius=rsun, frame=HeliographicCarrington(observer=observer, obstime=axes_frame.obstime)) visible, hidden = _plot_vertices(lon0, axes, axes_frame, rsun, close_path=False, **kwargs) return visible, hidden
def _plot_vertices(coord, axes, frame, rsun, close_path=True, **kwargs): """ Draws the provided SkyCoord on the WCSAxes as `~matplotlib.patches.Polygon` objects depending on visibility. """ c_kw = {'fill': False, 'color': 'white', 'zorder': 100} c_kw.update(kwargs) transform = axes.get_transform("world") c_kw.setdefault('transform', transform) # Get the 2D vertices of the coordinates coord = coord.transform_to(frame) Tx = coord.spherical.lon.to_value(u.deg) Ty = coord.spherical.lat.to_value(u.deg) vertices = np.array([Tx, Ty]).T # Determine which points are visible if hasattr(frame, 'observer'): # The reference distance is the distance to the limb for the axes # observer rsun = getattr(frame, 'rsun', rsun) reference_distance = np.sqrt(frame.observer.radius**2 - rsun**2) is_visible = coord.spherical.distance <= reference_distance else: # If the axes has no observer, the entire limb is considered visible is_visible = np.ones_like(coord.spherical.distance, bool, subok=False) # Identify discontinuities. Uses the same approach as # astropy.visualization.wcsaxes.grid_paths.get_lon_lat_path() step = np.sqrt((vertices[1:, 0] - vertices[:-1, 0]) ** 2 + (vertices[1:, 1] - vertices[:-1, 1]) ** 2) continuous = np.concatenate([[True, True], step[1:] < 100 * step[:-1]]) if not close_path is True: continuous[0] = False continuous[-1] = False visible, hidden = None, None if np.sum(is_visible) > 0: if 'linestyle' not in kwargs: c_kw['linestyle'] = '-' visible = patches.Polygon(vertices, **c_kw) _modify_polygon_visibility(visible, is_visible & continuous) axes.add_artist(visible) if np.sum(~is_visible) > 0: if 'linestyle' not in kwargs: c_kw['linestyle'] = ':' hidden = patches.Polygon(vertices, **c_kw) _modify_polygon_visibility(hidden, ~is_visible & continuous) axes.add_artist(hidden) return visible, hidden def _modify_polygon_visibility(polygon, keep): polygon_codes = polygon.get_path().codes polygon_codes[:-1][~keep] = Path.MOVETO polygon_codes[-1] = Path.MOVETO if not keep[0] else Path.LINETO