Source code for sunpy.map.mapbase

"""
Map is a generic Map class from which all other Map classes inherit from.
"""
import re
import copy
import html
import inspect
import numbers
import textwrap
import itertools
import webbrowser
from tempfile import NamedTemporaryFile
from collections import namedtuple

import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.backend_bases import FigureCanvasBase
from matplotlib.figure import Figure

try:
    from dask.array import Array as DaskArray
    DASK_INSTALLED = True
except ImportError:
    DASK_INSTALLED = False

import astropy.units as u
import astropy.wcs
from astropy.coordinates import BaseCoordinateFrame, Longitude, SkyCoord, UnitSphericalRepresentation
from astropy.nddata import NDData
from astropy.utils.metadata import MetaData
from astropy.visualization import AsymmetricPercentileInterval, HistEqStretch, ImageNormalize
from astropy.visualization.wcsaxes import Quadrangle, WCSAxes

# The next two are not used but are called to register functions with external modules
import sunpy.coordinates
import sunpy.io as io
import sunpy.io._fits
import sunpy.visualization.colormaps
from sunpy import config, log
from sunpy.coordinates import HeliographicCarrington, get_earth, sun
from sunpy.coordinates.utils import get_rectangle_coordinates
from sunpy.image.resample import resample as sunpy_image_resample
from sunpy.image.resample import reshape_image_to_4d_superpixel
from sunpy.image.transform import _get_transform_method, _rotation_function_names, affine_transform
from sunpy.sun import constants
from sunpy.time import is_time, parse_time
from sunpy.util import MetaDict, expand_list
from sunpy.util.decorators import (
    add_common_docstring,
    cached_property_based_on,
    check_arithmetic_compatibility,
)
from sunpy.util.exceptions import warn_metadata, warn_user
from sunpy.util.functools import seconddispatch
from sunpy.util.util import _figure_to_base64, fix_duplicate_notes
from sunpy.visualization import axis_labels_from_ctype, peek_show, wcsaxes_compat
from sunpy.visualization.colormaps import cm as sunpy_cm

TIME_FORMAT = config.get("general", "time_format")
PixelPair = namedtuple('PixelPair', 'x y')
SpatialPair = namedtuple('SpatialPair', 'axis1 axis2')

_META_FIX_URL = 'https://docs.sunpy.org/en/stable/how_to/fix_map_metadata.html'

# Manually specify the ``.meta`` docstring. This is assigned to the .meta
# class attribute in GenericMap.__init__()
_meta_doc = """
The map metadata.

This is used to interpret the map data. It may
have been modified from the original metadata by sunpy. See the
`~sunpy.util.MetaDict.added_items`, `~sunpy.util.MetaDict.removed_items`
and `~sunpy.util.MetaDict.modified_items` properties of MetaDict
to query how the metadata has been modified.
"""

# The notes live here so we can reuse it in the source maps
_notes_doc = """

Notes
-----

A number of the properties of this class are returned as two-value named
tuples that can either be indexed by position ([0] or [1]) or be accessed
by the names (.x and .y) or (.axis1 and .axis2). Things that refer to pixel
axes use the ``.x``, ``.y`` convention, where x and y refer to the FITS
axes (x for columns y for rows). Spatial axes use ``.axis1`` and ``.axis2``
which correspond to the first and second axes in the header. ``axis1``
corresponds to the coordinate axis for ``x`` and ``axis2`` corresponds to
``y``.

This class assumes that the metadata adheres to the FITS 4 standard.
Where the CROTA2 metadata is provided (without PC_ij) it assumes a conversion
to the standard PC_ij described in section 6.1 of .
`Calabretta & Greisen (2002) <https://doi.org/10.1051/0004-6361:20021327>`_

.. warning::
    If a header has CD_ij values but no PC_ij values, CDELT values are required
    for this class to construct the WCS.
    If a file with more than two dimensions is feed into the class,
    only the first two dimensions (NAXIS1, NAXIS2) will be loaded and the
    rest will be discarded.
"""

__all__ = ['GenericMap', 'MapMetaValidationError', 'PixelPair']


[docs] class MapMetaValidationError(AttributeError): pass
[docs] class GenericMap(NDData): """ A Generic spatially-aware 2D data array Parameters ---------- data : `numpy.ndarray`, list A 2d list or ndarray containing the map data. header : dict A dictionary of the original image header tags. plot_settings : dict, optional Plot settings. Other Parameters ---------------- **kwargs : Additional keyword arguments are passed to `~astropy.nddata.NDData` init. Examples -------- >>> import sunpy.map >>> import sunpy.data.sample # doctest: +REMOTE_DATA >>> aia = sunpy.map.Map(sunpy.data.sample.AIA_171_IMAGE) # doctest: +REMOTE_DATA >>> aia # doctest: +REMOTE_DATA <sunpy.map.sources.sdo.AIAMap object at ...> SunPy Map --------- Observatory: SDO Instrument: AIA 3 Detector: AIA Measurement: 171.0 Angstrom Wavelength: 171.0 Angstrom Observation Date: 2011-06-07 06:33:02 Exposure Time: 0.234256 s Dimension: [1024. 1024.] pix Coordinate System: helioprojective Scale: [2.402792 2.402792] arcsec / pix Reference Pixel: [511.5 511.5] pix Reference Coord: [3.22309951 1.38578135] arcsec array([[ -95.92475 , 7.076416 , -1.9656711, ..., -127.96519 , -127.96519 , -127.96519 ], [ -96.97533 , -5.1167884, 0. , ..., -98.924576 , -104.04137 , -127.919716 ], [ -93.99607 , 1.0189276, -4.0757103, ..., -5.094638 , -37.95505 , -127.87541 ], ..., [-128.01454 , -128.01454 , -128.01454 , ..., -128.01454 , -128.01454 , -128.01454 ], [-127.899666 , -127.899666 , -127.899666 , ..., -127.899666 , -127.899666 , -127.899666 ], [-128.03072 , -128.03072 , -128.03072 , ..., -128.03072 , -128.03072 , -128.03072 ]], dtype=float32) >>> aia.spatial_units # doctest: +REMOTE_DATA SpatialPair(axis1=Unit("arcsec"), axis2=Unit("arcsec")) >>> aia.peek() # doctest: +SKIP """ _registry = dict() # This overrides the default doc for the meta attribute meta = MetaData(doc=_meta_doc, copy=False) # Enabling the GenericMap reflected operators is a bit subtle. The GenericMap # reflected operator will be used only if the Quantity non-reflected operator # returns NotImplemented. The Quantity operator strips the unit from the # Quantity and tries to combine the value with the GenericMap using NumPy's # __array_ufunc__(). If NumPy believes that it can proceed, this will result # in an error. We explicitly set __array_ufunc__ = None so that the NumPy # call, and consequently the Quantity operator, will return NotImplemented. __array_ufunc__ = None def __init_subclass__(cls, **kwargs): """ An __init_subclass__ hook initializes all of the subclasses of a given class. So for each subclass, it will call this block of code on import. This replicates some metaclass magic without the need to be aware of metaclasses. Here we use this to register each subclass in a dict that has the ``is_datasource_for`` attribute. This is then passed into the Map Factory so we can register them. """ super().__init_subclass__(**kwargs) if cls.__doc__ is None: # Set an empty string, to prevent an error adding None to str in the next line cls.__doc__ = '' cls.__doc__ = fix_duplicate_notes(_notes_doc, cls.__doc__) if hasattr(cls, 'is_datasource_for'): # NOTE: This conditional is due to overlapping map sources in sunpy and pfsspy that # lead to a MultipleMatchError if sunpy.map and pfsspy.map are imported. # See https://github.com/sunpy/sunpy/issues/7294 for more information. if f'{cls.__module__}.{cls.__name__}' != "pfsspy.map.GongSynopticMap": cls._registry[cls] = cls.is_datasource_for def __init__(self, data, header, plot_settings=None, **kwargs): # If the data has more than two dimensions, the first dimensions # (NAXIS1, NAXIS2) are used and the rest are discarded. ndim = data.ndim if ndim > 2: # We create a slice that removes all but the 'last' two # dimensions. (Note dimensions in ndarray are in reverse order) new_2d_slice = [0]*(ndim-2) new_2d_slice.extend([slice(None), slice(None)]) data = data[tuple(new_2d_slice)] # Warn the user that the data has been truncated warn_user("This file contains more than 2 dimensions. " "Data will be truncated to the first two dimensions.") params = list(inspect.signature(NDData).parameters) nddata_kwargs = {x: kwargs.pop(x) for x in params & kwargs.keys()} super().__init__(data, meta=MetaDict(header), **nddata_kwargs) # Setup some attributes self._nickname = None # These are placeholders for default attributes, which are only set # once if their data isn't present in the map metadata. self._default_time = None self._default_dsun = None self._default_carrington_longitude = None self._default_heliographic_latitude = None self._default_heliographic_longitude = None # Validate header # TODO: This should be a function of the header, not of the map self._validate_meta() if self.dtype == np.uint8: norm = None else: # Put import here to reduce sunpy.map import time from matplotlib import colors norm = colors.Normalize() # Visualization attributes self.plot_settings = {'cmap': 'gray', 'norm': norm, 'interpolation': 'nearest', 'origin': 'lower' } if plot_settings: self.plot_settings.update(plot_settings) # Try and set the colormap. This is not always possible if this method # is run before map sources fix some of their metadata, so # just ignore any exceptions raised. try: cmap = self._get_cmap_name() if cmap in sunpy_cm.cmlist: self.plot_settings['cmap'] = cmap except Exception: pass def __getitem__(self, key): """ This should allow indexing by physical coordinate """ raise NotImplementedError( "The ability to index Map by physical" " coordinate is not yet implemented.") def _text_summary(self): dt = self.exposure_time wave = self.wavelength measurement = self.measurement dt = 'Unknown' if dt is None else dt wave = 'Unknown' if wave is None else wave measurement = 'Unknown' if measurement is None else measurement return textwrap.dedent("""\ SunPy Map --------- Observatory:\t\t {obs} Instrument:\t\t {inst} Detector:\t\t {det} Measurement:\t\t {meas} Wavelength:\t\t {wave} Observation Date:\t {date} Exposure Time:\t\t {dt} Dimension:\t\t {dim} Coordinate System:\t {coord} Scale:\t\t\t {scale} Reference Pixel:\t {refpix} Reference Coord:\t {refcoord}\ """).format(obs=self.observatory, inst=self.instrument, det=self.detector, meas=measurement, wave=wave, date=self.date.strftime(TIME_FORMAT), dt=dt, dim=u.Quantity(self.dimensions), scale=u.Quantity(self.scale), coord=self._coordinate_frame_name, refpix=u.Quantity(self.reference_pixel), refcoord=u.Quantity((self._reference_longitude, self._reference_latitude)), tmf=TIME_FORMAT) def __str__(self): return f"{self._text_summary()}\n{self.data.__repr__()}" def __repr__(self): return f"{object.__repr__(self)}\n{self}" def _repr_html_(self, compute_dask=False): """ Produce an HTML summary with plots for use in Jupyter notebooks. """ # Convert the text repr to an HTML table partial_html = self._text_summary()[20:].replace('\n', '</td></tr><tr><th>')\ .replace(':\t', '</th><td>') text_to_table = textwrap.dedent(f"""\ <table style='text-align:left'> <tr><th>{partial_html}</td></tr> </table>""").replace('\n', '') # Handle bad values (infinite and NaN) in the data array finite_data = self.data[np.isfinite(self.data)] count_nan = np.isnan(self.data).sum() count_inf = np.isinf(self.data).sum() if DASK_INSTALLED and isinstance(finite_data, DaskArray): # This will fetch the entire data array into memory and only happens for the quicklook method if compute_dask: finite_data = finite_data.compute() else: dask_html = self.data._repr_html_() return textwrap.dedent(f"""\ <pre>{html.escape(object.__repr__(self))}</pre> <table> <tr> <td>{text_to_table}</td> <td> {dask_html} </td> </tr> <tr> </tr> </table>""") # Assemble an informational string with the counts of bad pixels bad_pixel_text = "" if count_nan + count_inf > 0: bad_pixel_text = "Bad pixels are shown in red: " text_list = [] if count_nan > 0: text_list.append(f"{count_nan} NaN") if count_inf > 0: text_list.append(f"{count_inf} infinite") bad_pixel_text += ", ".join(text_list) # Use a grayscale colormap with histogram equalization (and red for bad values) # Make a copy of the colormap to avoid modifying the matplotlib instance when # doing set_bad() (copy not needed when min mpl is 3.5, as already a copy) cmap = copy.copy(matplotlib.colormaps['gray']) cmap.set_bad(color='red') norm = ImageNormalize(stretch=HistEqStretch(finite_data)) # Plot the image in pixel space fig = Figure(figsize=(5.2, 4.8)) # Figure instances in matplotlib<3.1 do not create a canvas by default if fig.canvas is None: FigureCanvasBase(fig) ax = fig.subplots() ax.imshow(self.data, origin='lower', interpolation='nearest', cmap=cmap, norm=norm) ax.set_xlabel('X pixel') ax.set_ylabel('Y pixel') ax.set_title('In pixel space') pixel_src = _figure_to_base64(fig) bounds = ax.get_position().bounds # save these axes bounds for later use # Plot the image using WCS information, with the same axes bounds as above fig = Figure(figsize=(5.2, 4.8)) # Figure instances in matplotlib<3.1 do not create a canvas by default if fig.canvas is None: FigureCanvasBase(fig) # Create the WCSAxes manually because we need to avoid using pyplot ax = WCSAxes(fig, bounds, aspect='equal', wcs=self.wcs) fig.add_axes(ax) self.plot(axes=ax, cmap=cmap, norm=norm) ax.set_title('In coordinate space using WCS information') wcs_src = _figure_to_base64(fig) # Plot the histogram of pixel values fig = Figure(figsize=(4.8, 2.4), constrained_layout=True) # Figure instances in matplotlib<3.1 do not create a canvas by default if fig.canvas is None: FigureCanvasBase(fig) ax = fig.subplots() values, bins, patches = ax.hist(finite_data.ravel(), bins=100) norm_centers = norm(0.5 * (bins[:-1] + bins[1:])).data for c, p in zip(norm_centers, patches): plt.setp(p, "facecolor", cmap(c)) ax.plot(np.array([bins[:-1], bins[1:]]).T.ravel(), np.array([values, values]).T.ravel()) ax.set_facecolor('white') ax.semilogy() # Explicitly set the power limits for the X axis formatter to avoid text overlaps ax.xaxis.get_major_formatter().set_powerlimits((-3, 4)) ax.set_xlabel(f"Pixel value{' (' + str(self.unit) + ')' if self.unit else ''} in linear bins") ax.set_ylabel('# of pixels') ax.set_title('Distribution of pixel values [click for cumulative]') hist_src = _figure_to_base64(fig) # Plot the CDF of the pixel values using a symmetric-log horizontal scale fig = Figure(figsize=(4.8, 2.4), constrained_layout=True) # TODO: Figure instances in matplotlib<3.1 do not create a canvas by default if fig.canvas is None: FigureCanvasBase(fig) ax = fig.subplots() n_bins = 256 bins = norm.inverse(np.arange(n_bins + 1) / n_bins) values, _, patches = ax.hist(finite_data.ravel(), bins=bins, cumulative=True) for i, p in enumerate(patches): plt.setp(p, "facecolor", cmap((i + 0.5) / n_bins)) ax.plot(np.array([bins[:-1], bins[1:]]).T.ravel(), np.array([values, values]).T.ravel()) ax.set_facecolor('white') ax.set_xscale('symlog') ax.set_yscale('log') ax.set_xlabel(f"Pixel value{' (' + str(self.unit) + ')' if self.unit else ''} in equalized bins") ax.set_ylabel('Cumulative # of pixels') ax.set_title('Cumulative distribution of pixel values') cdf_src = _figure_to_base64(fig) return textwrap.dedent(f"""\ <pre>{html.escape(object.__repr__(self))}</pre> <table> <tr> <td>{text_to_table}</td> <td rowspan=3> <div align=center> Image colormap uses histogram equalization<br> Click on the image to toggle between units </div> <img src='data:image/png;base64,{wcs_src}' src2='data:image/png;base64,{pixel_src}' onClick='var temp = this.src; this.src = this.getAttribute("src2"); this.setAttribute("src2", temp)' /> <div align=center> {bad_pixel_text} </div> </td> </tr> <tr> </tr> <tr> <td><img src='data:image/png;base64,{hist_src}' src2='data:image/png;base64,{cdf_src}' onClick='var temp = this.src; this.src = this.getAttribute("src2"); this.setAttribute("src2", temp)' /> </td> </tr> </table>""")
[docs] def quicklook(self): """ Display a quicklook summary of the Map instance using the default web browser. Notes ----- The image colormap uses `histogram equalization <https://en.wikipedia.org/wiki/Histogram_equalization>`__. Clicking on the image to switch between pixel space and coordinate space requires Javascript support to be enabled in the web browser. Examples -------- >>> from sunpy.map import Map >>> import sunpy.data.sample # doctest: +REMOTE_DATA >>> smap = Map(sunpy.data.sample.AIA_171_IMAGE) # doctest: +REMOTE_DATA >>> smap.quicklook() # doctest: +SKIP (which will open the following content in the default web browser) .. generate:: html :html_border: from sunpy.map import Map import sunpy.data.sample smap = Map(sunpy.data.sample.AIA_171_IMAGE) print(smap._repr_html_()) """ with NamedTemporaryFile('w', delete=False, prefix='sunpy.map.', suffix='.html') as f: url = 'file://' + f.name f.write(textwrap.dedent(f"""\ <html> <title>Quicklook summary for {html.escape(object.__repr__(self))}</title> <body>{self._repr_html_(compute_dask=True)}</body> </html>""")) webbrowser.open_new_tab(url)
@classmethod def _new_instance(cls, data, meta, plot_settings=None, **kwargs): """ Instantiate a new instance of this class using given data. This is a shortcut for ``type(self)(data, meta, plot_settings)``. """ new_map = cls(data, meta, **kwargs) # plot_settings are set explicitly here as some map sources # explicitly set some of the plot_settings in the constructor # and we want to preserve the plot_settings of the previous # instance. if plot_settings is not None: new_map.plot_settings.update(plot_settings) return new_map def _get_lon_lat(self, frame): """ Given a coordinate frame, extract the lon and lat by casting to SphericalRepresentation first. """ r = frame.represent_as(UnitSphericalRepresentation) return r.lon.to(self.spatial_units[0]), r.lat.to(self.spatial_units[1]) @property def quantity(self): """Unitful representation of the map data.""" return u.Quantity(self.data, self.unit, copy=False) def _new_instance_from_op(self, new_data): """ Helper function for creating new map instances after arithmetic operations. """ new_meta = copy.deepcopy(self.meta) new_meta['bunit'] = new_data.unit.to_string('fits') return self._new_instance(new_data.value, new_meta, plot_settings=self.plot_settings) def __neg__(self): return self._new_instance(-self.data, self.meta, plot_settings=self.plot_settings) @check_arithmetic_compatibility def __pow__(self, value): new_data = self.quantity ** value return self._new_instance_from_op(new_data) @check_arithmetic_compatibility def __add__(self, value): new_data = self.quantity + value return self._new_instance_from_op(new_data) def __radd__(self, value): return self.__add__(value) def __sub__(self, value): return self.__add__(-value) def __rsub__(self, value): return self.__neg__().__add__(value) @check_arithmetic_compatibility def __mul__(self, value): new_data = self.quantity * value return self._new_instance_from_op(new_data) def __rmul__(self, value): return self.__mul__(value) def __truediv__(self, value): return self.__mul__(1/value) @check_arithmetic_compatibility def __rtruediv__(self, value): new_data = value / self.quantity return self._new_instance_from_op(new_data) @property def _meta_hash(self): return self.meta.item_hash() def _set_symmetric_vmin_vmax(self): """ Set symmetric vmin and vmax about zero """ threshold = np.nanmax(abs(self.data)) self.plot_settings['norm'].vmin = -threshold self.plot_settings['norm'].vmax = threshold @property @cached_property_based_on('_meta_hash') def wcs(self): """ The `~astropy.wcs.WCS` property of the map. """ w2 = astropy.wcs.WCS(naxis=2) # Add one to go from zero-based to one-based indexing w2.wcs.crpix = u.Quantity(self.reference_pixel) + 1 * u.pix # Make these a quantity array to prevent the numpy setting element of # array with sequence error. # Explicitly call ``.to()`` to check that scale is in the correct units w2.wcs.cdelt = u.Quantity([self.scale[0].to(self.spatial_units[0] / u.pix), self.scale[1].to(self.spatial_units[1] / u.pix)]) w2.wcs.crval = u.Quantity([self._reference_longitude, self._reference_latitude]) w2.wcs.ctype = self.coordinate_system w2.wcs.pc = self.rotation_matrix w2.wcs.set_pv(self._pv_values) # FITS standard doesn't allow both PC_ij *and* CROTA keywords w2.wcs.crota = (0, 0) w2.wcs.cunit = self.spatial_units w2.wcs.dateobs = self.date.isot w2.wcs.aux.rsun_ref = self.rsun_meters.to_value(u.m) # Set observer coordinate information except when we know it is not appropriate (e.g., HGS) sunpy_frame = sunpy.coordinates.wcs_utils._sunpy_frame_class_from_ctypes(w2.wcs.ctype) if sunpy_frame is None or hasattr(sunpy_frame, 'observer'): # Clear all the aux information that was set earlier. This is to avoid # issues with maps that store multiple observer coordinate keywords. # Note that we have to create a new WCS as it's not possible to modify # wcs.wcs.aux in place. header = w2.to_header() for kw in ['crln_obs', 'dsun_obs', 'hgln_obs', 'hglt_obs']: header.pop(kw, None) w2 = astropy.wcs.WCS(header) # Get observer coord, and set the aux information obs_coord = self.observer_coordinate sunpy.coordinates.wcs_utils._set_wcs_aux_obs_coord(w2, obs_coord) # Set the shape of the data array w2.array_shape = self.data.shape # Validate the WCS here. w2.wcs.set() return w2 @property def coordinate_frame(self): """ An `astropy.coordinates.BaseCoordinateFrame` instance created from the coordinate information for this Map, or None if the frame cannot be determined. """ try: return astropy.wcs.utils.wcs_to_celestial_frame(self.wcs) except ValueError as e: warn_user(f'Could not determine coordinate frame from map metadata.\n{e}') return None @property def _coordinate_frame_name(self): if self.coordinate_frame is None: return 'Unknown' return self.coordinate_frame.name def _as_mpl_axes(self): """ Compatibility hook for Matplotlib and WCSAxes. This functionality requires the WCSAxes package to work. The reason we include this here is that it allows users to use WCSAxes without having to explicitly import WCSAxes With this method, one can do:: import matplotlib.pyplot as plt import sunpy.map amap = sunpy.map.Map('filename.fits') fig = plt.figure() ax = plt.subplot(projection=amap) ... and this will generate a plot with the correct WCS coordinates on the axes. See <https://docs.astropy.org/en/stable/visualization/wcsaxes/index.html> for more information. """ # This code is reused from Astropy return WCSAxes, {'wcs': self.wcs} # Some numpy extraction @property def dimensions(self): """ The dimensions of the array (x axis first, y axis second). """ return PixelPair(*u.Quantity(np.flipud(self.data.shape), 'pixel')) @property def dtype(self): """ The `numpy.dtype` of the array of the map. """ return self.data.dtype @property def ndim(self): """ The value of `numpy.ndarray.ndim` of the data array of the map. """ return self.data.ndim
[docs] def std(self, *args, **kwargs): """ Calculate the standard deviation of the data array, ignoring NaNs. """ return np.nanstd(self.data, *args, **kwargs)
[docs] def mean(self, *args, **kwargs): """ Calculate the mean of the data array, ignoring NaNs. """ return np.nanmean(self.data, *args, **kwargs)
[docs] def min(self, *args, **kwargs): """ Calculate the minimum value of the data array, ignoring NaNs. """ return np.nanmin(self.data, *args, **kwargs)
[docs] def max(self, *args, **kwargs): """ Calculate the maximum value of the data array, ignoring NaNs. """ return np.nanmax(self.data, *args, **kwargs)
@staticmethod def _parse_fits_unit(unit_str): replacements = {'gauss': 'G', 'dn': 'ct', 'dn/s': 'ct/s', 'counts / pixel': 'ct/pix',} if unit_str.lower() in replacements: unit_str = replacements[unit_str.lower()] unit = u.Unit(unit_str, format='fits', parse_strict='silent') if isinstance(unit, u.UnrecognizedUnit): warn_metadata(f'Could not parse unit string "{unit_str}" as a valid FITS unit.\n' f'See {_META_FIX_URL} for how to fix metadata before loading it ' 'with sunpy.map.Map.\n' 'See https://fits.gsfc.nasa.gov/fits_standard.html for ' 'the FITS unit standards.') unit = None return unit @property def unit(self): """ Unit of the map data. This is taken from the 'BUNIT' FITS keyword. If no 'BUNIT' entry is present in the metadata then this returns `None`. If the 'BUNIT' value cannot be parsed into a unit a warning is raised, and `None` returned. """ unit_str = self.meta.get('bunit', None) if unit_str is None: return return self._parse_fits_unit(unit_str) # #### Keyword attribute and other attribute definitions #### # def _base_name(self): """Abstract the shared bit between name and latex_name""" if self.measurement is None: format_str = "{nickname} {date}" else: format_str = "{nickname} {{measurement}} {date}" return format_str.format(nickname=self.nickname, date=parse_time(self.date).strftime(TIME_FORMAT)) @property def name(self): """Human-readable description of the Map.""" return self._base_name().format(measurement=self.measurement) @property def latex_name(self): """LaTeX formatted description of the Map.""" if isinstance(self.measurement, u.Quantity): return self._base_name().format(measurement=self.measurement._repr_latex_()) else: return self.name @property def nickname(self): """An abbreviated human-readable description of the map-type; part of the Helioviewer data model.""" return self._nickname if self._nickname else self.detector @nickname.setter def nickname(self, n): self._nickname = n def _get_date(self, key): time = self.meta.get(key, None) if not time: return # Get the time scale if 'TAI' in time: # SDO specifies the 'TAI' scale in their time string, which is parsed # by parse_time(). If a different timescale is also present, warn the # user that it will be ignored. timesys = 'TAI' timesys_meta = self.meta.get('timesys', '').upper() if timesys_meta not in ('', 'TAI'): warn_metadata('Found "TAI" in time string, ignoring TIMESYS keyword ' f'which is set to "{timesys_meta}".') else: timesys = self._timesys return parse_time(time, scale=timesys.lower()) @property def _timesys(self): """ Time system. """ # UTC is the FITS standard default return self.meta.get('timesys', 'UTC') @property def date_start(self): """ Time of the beginning of the image acquisition. Taken from the DATE-BEG FITS keyword. """ return self._get_date('date-beg') @property def date_end(self): """ Time of the end of the image acquisition. Taken from the DATE-END FITS keyword. """ return self._get_date('date-end') @property def date_average(self): """ Average time of the image acquisition. Taken from the DATE-AVG FITS keyword if present, otherwise halfway between `date_start` and `date_end` if both pieces of metadata are present. """ avg = self._get_date('date-avg') if avg is None: start, end = self.date_start, self.date_end if start is not None and end is not None: avg = start + (end - start) / 2 return avg @property def _date_obs(self): # Get observation date from date-obs, falling back to date_obs time = self._get_date('date-obs') if is_time(self.meta.get('date_obs', None)): time = time or self._get_date('date_obs') return time @property def date(self): """ Image observation time. For different combinations of map metadata this can return either the start time, end time, or a time between these. It is recommended to use `~sunpy.map.GenericMap.date_average`, `~sunpy.map.GenericMap.date_start`, or `~sunpy.map.GenericMap.date_end` instead if you need one of these specific times. Taken from, in order of preference: 1. The DATE-OBS FITS keyword 2. `~sunpy.map.GenericMap.date_average` 3. `~sunpy.map.GenericMap.date_start` 4. `~sunpy.map.GenericMap.date_end` 5. The current time """ time = self._date_obs time = time or self.date_average time = time or self.date_start time = time or self.date_end if time is None: if self._default_time is None: warn_metadata("Missing metadata for observation time, " "setting observation time to current time. " "Set the 'DATE-AVG' FITS keyword to prevent this warning.") self._default_time = parse_time('now') time = self._default_time return time @property def detector(self): """ Detector name. This is taken from the 'DETECTOR' FITS keyword. """ return self.meta.get('detector', "") @property def timeunit(self): """ The `~astropy.units.Unit` of the exposure time of this observation. Taken from the "TIMEUNIT" FITS keyword, and defaults to seconds (as per) the FITS standard). """ return u.Unit(self.meta.get('timeunit', 's')) @property def exposure_time(self): """ Exposure time of the image. This is taken from the 'XPOSURE' keyword or the 'EXPTIME' FITS keyword, in that order. """ exptime = self.meta.get('xposure') or self.meta.get('exptime') if exptime is not None: return exptime * self.timeunit @property def instrument(self): """Instrument name.""" return self.meta.get('instrume', "").replace("_", " ") @property def measurement(self): """ Measurement wavelength. This is taken from the 'WAVELNTH' FITS keywords. If the keyword is not present, defaults to `None`. If 'WAVEUNIT' keyword isn't present, defaults to dimensionless units. """ return self.wavelength @property def waveunit(self): """ The `~astropy.units.Unit` of the wavelength of this observation. This is taken from the 'WAVEUNIT' FITS keyword. If the keyword is not present, defaults to `None` """ if 'waveunit' in self.meta: return u.Unit(self.meta['waveunit']) else: wunit = sunpy.io._fits.extract_waveunit(self.meta) if wunit is not None: return u.Unit(wunit) @property def wavelength(self): """ Wavelength of the observation. This is taken from the 'WAVELNTH' FITS keywords. If the keyword is not present, defaults to `None`. If 'WAVEUNIT' keyword isn't present, defaults to dimensionless units. """ if 'wavelnth' in self.meta: return u.Quantity(self.meta['wavelnth'], self.waveunit) @property def observatory(self): """ Observatory or Telescope name. This is taken from the 'OBSRVTRY' FITS keyword. """ return self.meta.get('obsrvtry', self.meta.get('telescop', "")).replace("_", " ") @property def processing_level(self): """ Returns the FITS processing level if present. This is taken from the 'LVL_NUM' FITS keyword. """ return self.meta.get('lvl_num', None) @property def bottom_left_coord(self): """ The physical coordinate at the center of the bottom left ([0, 0]) pixel. """ return self.wcs.pixel_to_world(0, 0) @property def top_right_coord(self): """ The physical coordinate at the center of the the top right ([-1, -1]) pixel. """ top_right = np.array([self.dimensions.x.value, self.dimensions.y.value]) - 1 return self.wcs.pixel_to_world(*top_right) @property def center(self): """ Return a coordinate object for the center pixel of the array. If the array has an even number of pixels in a given dimension, the coordinate returned lies on the edge between the two central pixels. """ center = (np.array([self.dimensions.x.value, self.dimensions.y.value]) - 1) / 2. return self.wcs.pixel_to_world(*center)
[docs] @u.quantity_input def shift_reference_coord(self, axis1: u.deg, axis2: u.deg): """ Returns a map shifted by a specified amount to, for example, correct for a bad map location. These values are applied directly to the `~sunpy.map.GenericMap.reference_coordinate`. To check how much the reference coordinate has been modified, see ``sunpy.map.GenericMap.meta.modified_items['CRVAL1']`` and ``sunpy.map.GenericMap.meta.modified_items['CRVAL2']``. Parameters ---------- axis1 : `~astropy.units.Quantity` The shift to apply to the Longitude (solar-x) coordinate. axis2 : `~astropy.units.Quantity` The shift to apply to the Latitude (solar-y) coordinate Returns ------- out : `~sunpy.map.GenericMap` or subclass A new shifted Map. """ new_meta = self.meta.copy() # Update crvals new_meta['crval1'] = ((self._reference_longitude + axis1).to(self.spatial_units[0])).value new_meta['crval2'] = ((self._reference_latitude + axis2).to(self.spatial_units[1])).value # Create new map with the modification new_map = self._new_instance(self.data, new_meta, self.plot_settings) return new_map
def _rsun_meters(self, dsun=None): """ This property exists to avoid circular logic in constructing the observer coordinate, by allowing a custom 'dsun' to be specified, instead of one extracted from the `.observer_coordinate` property. """ rsun = self.meta.get('rsun_ref', None) if rsun is not None: return rsun * u.m elif self._rsun_obs_no_default is not None: if dsun is None: dsun = self.dsun return sun._radius_from_angular_radius(self.rsun_obs, dsun) else: log.info("Missing metadata for solar radius: assuming " "the standard radius of the photosphere.") return constants.radius @property def rsun_meters(self): """ Assumed radius of observed emission from the Sun center. This is taken from the RSUN_REF FITS keyword, if present. If not, and angular radius metadata is present, it is calculated from `~sunpy.map.GenericMap.rsun_obs` and `~sunpy.map.GenericMap.dsun`. If neither pieces of metadata are present, defaults to the standard photospheric radius. """ return self._rsun_meters() @property def _rsun_obs_no_default(self): """ Get the angular radius value from FITS keywords without defaulting. Exists to avoid circular logic in `rsun_meters()` above. """ return self.meta.get('rsun_obs', self.meta.get('solar_r', self.meta.get('radius', None))) @property def rsun_obs(self): """ Angular radius of the observation from Sun center. This value is taken (in order of preference) from the 'RSUN_OBS', 'SOLAR_R', or 'RADIUS' FITS keywords. If none of these keys are present, the angular radius is calculated from `~sunpy.map.GenericMap.rsun_meters` and `~sunpy.map.GenericMap.dsun`. """ rsun_arcseconds = self._rsun_obs_no_default if rsun_arcseconds is not None: return rsun_arcseconds * u.arcsec else: return sun._angular_radius(self.rsun_meters, self.dsun) @property def coordinate_system(self): """ Coordinate system used for x and y axes (ctype1/2). If not present, defaults to (HPLN-TAN, HPLT-TAN), and emits a warning. """ ctype1 = self.meta.get('ctype1', None) if not ctype1: warn_metadata("Missing CTYPE1 from metadata, assuming CTYPE1 is HPLN-TAN") ctype1 = 'HPLN-TAN' ctype2 = self.meta.get('ctype2', None) if not ctype2: warn_metadata("Missing CTYPE2 from metadata, assuming CTYPE2 is HPLT-TAN") ctype2 = 'HPLT-TAN' # Astropy WCS does not understand the SOHO default of "solar-x" and # "solar-y" ctypes. This overrides the default assignment and # changes it to a ctype that is understood. See Thompson, 2006, A.&A., # 449, 791. if ctype1.lower() in ("solar-x", "solar_x"): ctype1 = 'HPLN-TAN' if ctype2.lower() in ("solar-y", "solar_y"): ctype2 = 'HPLT-TAN' return SpatialPair(ctype1, ctype2) @property def _supported_observer_coordinates(self): """ A list of supported coordinate systems. This is a list so it can easily maintain a strict order. The list of two element tuples, the first item in the tuple is the keys that need to be in the header to use this coordinate system and the second is the kwargs to SkyCoord. """ return [(('hgln_obs', 'hglt_obs', 'dsun_obs'), {'lon': self.meta.get('hgln_obs'), 'lat': self.meta.get('hglt_obs'), 'radius': self.meta.get('dsun_obs'), 'unit': (u.deg, u.deg, u.m), 'frame': "heliographic_stonyhurst"}), (('crln_obs', 'crlt_obs', 'dsun_obs'), {'lon': self.meta.get('crln_obs'), 'lat': self.meta.get('crlt_obs'), 'radius': self.meta.get('dsun_obs'), 'unit': (u.deg, u.deg, u.m), 'frame': "heliographic_carrington"}), ] @property def _default_observer_coordinate(self): """ The default observer coordinate to use when there is insufficient information in the metadata. This should be overridden by map sources as appropriate. """ def _remove_existing_observer_location(self): """ Remove all keys that this map might use for observer location. """ all_keys = expand_list([e[0] for e in self._supported_observer_coordinates]) for key in all_keys: self.meta.pop(key) @property @cached_property_based_on('_meta_hash') def observer_coordinate(self): """ The Heliographic Stonyhurst Coordinate of the observer. """ warning_message = [] for keys, kwargs in self._supported_observer_coordinates: missing_keys = set(keys) - self.meta.keys() if not missing_keys: sc = SkyCoord(obstime=self.date, **kwargs) # If the observer location is supplied in Carrington coordinates, # the coordinate's `observer` attribute should be set to "self" if isinstance(sc.frame, HeliographicCarrington): sc.frame._observer = "self" sc = sc.heliographic_stonyhurst # We set rsun after constructing the coordinate, as we need # the observer-Sun distance (sc.radius) to calculate this, which # may not be provided directly in metadata (if e.g. the # observer coordinate is specified in a cartesian # representation) return SkyCoord(sc.replicate(rsun=self._rsun_meters(sc.radius))) elif missing_keys != keys: frame = kwargs['frame'] if isinstance(kwargs['frame'], str) else kwargs['frame'].name warning_message.append(f"For frame '{frame}' the following metadata is missing: " f"{','.join(missing_keys)}") default = self._default_observer_coordinate if default is not None: # If a map source specifies a default observer, we log a message at the debug level warning_message = (["Missing metadata for observer: assuming custom default observer."] + warning_message) log.debug("\n".join(warning_message)) return default else: # If a map source does not specify a default observer, we assume Earth center and warn warning_message = (["Missing metadata for observer: assuming Earth-based observer."] + warning_message + [""]) warn_metadata("\n".join(warning_message), stacklevel=3) return get_earth(self.date) @property def heliographic_latitude(self): """Observer heliographic latitude.""" return self.observer_coordinate.lat @property def heliographic_longitude(self): """Observer heliographic longitude.""" return self.observer_coordinate.lon @property def carrington_latitude(self): """Observer Carrington latitude.""" hgc_frame = HeliographicCarrington(observer=self.observer_coordinate, obstime=self.date, rsun=self.rsun_meters) return self.observer_coordinate.transform_to(hgc_frame).lat @property def carrington_longitude(self): """Observer Carrington longitude.""" hgc_frame = HeliographicCarrington(observer=self.observer_coordinate, obstime=self.date, rsun=self.rsun_meters) return self.observer_coordinate.transform_to(hgc_frame).lon @property def dsun(self): """Observer distance from the center of the Sun.""" return self.observer_coordinate.radius.to('m') @property def _reference_longitude(self): """ FITS-WCS compatible longitude. Used in self.wcs and self.reference_coordinate. """ return self.meta.get('crval1', 0.) * self.spatial_units[0] @property def _reference_latitude(self): return self.meta.get('crval2', 0.) * self.spatial_units[1] @property def reference_coordinate(self): """Reference point WCS axes in data units (i.e. crval1, crval2). This value includes a shift if one is set.""" return SkyCoord(self._reference_longitude, self._reference_latitude, frame=self.coordinate_frame) @property def reference_pixel(self): """ Pixel of reference coordinate. The pixel returned uses zero-based indexing, so will be 1 pixel less than the FITS CRPIX values. """ naxis1 = self.meta.get('naxis1', self.data.shape[1]) naxis2 = self.meta.get('naxis2', self.data.shape[0]) return PixelPair((self.meta.get('crpix1', (naxis1 + 1) / 2.) - 1) * u.pixel, (self.meta.get('crpix2', (naxis2 + 1) / 2.) - 1) * u.pixel) @property def scale(self): """ Image scale along the x and y axes in units/pixel (i.e. cdelt1, cdelt2). If the CDij matrix is defined but no CDELTi values are explicitly defined, effective CDELTi values are constructed from the CDij matrix. The effective CDELTi values are chosen so that each row of the PCij matrix has unity norm. This choice is optimal if the PCij matrix is a pure rotation matrix, but may not be as optimal if the PCij matrix includes any skew. """ if 'cd1_1' in self.meta and 'cdelt1' not in self.meta and 'cdelt2' not in self.meta: cdelt1 = np.sqrt(self.meta['cd1_1']**2 + self.meta['cd1_2']**2) cdelt2 = np.sqrt(self.meta['cd2_1']**2 + self.meta['cd2_2']**2) else: cdelt1 = self.meta.get('cdelt1', 1.) cdelt2 = self.meta.get('cdelt2', 1.) return SpatialPair(cdelt1 * self.spatial_units[0] / u.pixel, cdelt2 * self.spatial_units[1] / u.pixel) @property def spatial_units(self): """ Image coordinate units along the x and y axes (i.e. cunit1, cunit2). """ units = self.meta.get('cunit1', None), self.meta.get('cunit2', None) units = [None if unit is None else u.Unit(unit.lower()) for unit in units] return SpatialPair(units[0], units[1]) @property def rotation_matrix(self): r""" Matrix describing the transformation needed to align the reference pixel with the coordinate axes. The order or precedence of FITS keywords which this is taken from is: - PC\*_\* - CD\*_\* - CROTA\* Notes ----- In many cases this is a simple rotation matrix, hence the property name. It general it does not have to be a pure rotation matrix, and can encode other transformations e.g., skews for non-orthgonal coordinate systems. """ if any(key in self.meta for key in ['PC1_1', 'PC1_2', 'PC2_1', 'PC2_2']): return np.array( [ [self.meta.get('PC1_1', 1), self.meta.get('PC1_2', 0)], [self.meta.get('PC2_1', 0), self.meta.get('PC2_2', 1)] ] ) elif any(key in self.meta for key in ['CD1_1', 'CD1_2', 'CD2_1', 'CD2_2']): cd = np.array( [ [self.meta.get('CD1_1', 0), self.meta.get('CD1_2', 0)], [self.meta.get('CD2_1', 0), self.meta.get('CD2_2', 0)] ] ) cdelt = u.Quantity(self.scale).value # Divide each row by each CDELT return cd / np.expand_dims(cdelt, axis=1) else: return self._rotation_matrix_from_crota() @staticmethod def _pc_matrix(lam, angle): """ Returns PC matrix from the scale ration (lam) and rotation angle in radians (angle). """ return np.array([[np.cos(angle), -1 * lam * np.sin(angle)], [1/lam * np.sin(angle), np.cos(angle)]]) def _rotation_matrix_from_crota(self, crota_key='CROTA2'): """ This method converts the deprecated CROTA FITS kwargs to the new PC rotation matrix. This method can be overridden if an instruments header does not use this conversion. Parameters ---------- crota_key : str, optional The key to use for CROTA2. Defaults to 'CROTA2'. Notes ----- If the specified key isn't present in the metadata, a default rotation of 0deg is returned. """ lam = self.scale[1] / self.scale[0] p = np.deg2rad(self.meta.get(crota_key, 0)) return self._pc_matrix(lam, p) @property def _pv_values(self): """ Return any PV values in the metadata. """ pattern = re.compile('pv[0-9]_[0-9]a', re.IGNORECASE) pv_keys = [k for k in self.meta.keys() if pattern.match(k)] pv_values = [] for k in pv_keys: i, m = int(k[2]), int(k[4]) pv_values.append((i, m, self.meta[k])) return pv_values @property def fits_header(self): """ A `~astropy.io.fits.Header` representation of the ``meta`` attribute. """ return sunpy.io._fits.header_to_fits(self.meta) # #### Miscellaneous #### # def _get_cmap_name(self): """Build the default color map name.""" cmap_string = (self.observatory + self.detector + str(int(self.wavelength.to('angstrom').value))) return cmap_string.lower() def _validate_meta(self): """ Validates some meta-information associated with a Map. This method includes very basic validation checks which apply to all of the kinds of files that sunpy can read. Datasource-specific validation should be handled in the relevant file in the sunpy.map.sources package. """ msg = ('Image coordinate units for axis {} not present in metadata.') err_message = [] for i in [0, 1]: if self.spatial_units[i] is None: err_message.append(msg.format(i+1, i+1)) if err_message: err_message.append( f'See {_META_FIX_URL} for instructions on how to add missing metadata.') raise MapMetaValidationError('\n'.join(err_message)) for meta_property in ('waveunit', ): if (self.meta.get(meta_property) and u.Unit(self.meta.get(meta_property), parse_strict='silent').physical_type == 'unknown'): warn_metadata(f"Unknown value for {meta_property.upper()}.") if (self.coordinate_system[0].startswith(('SOLX', 'SOLY')) or self.coordinate_system[1].startswith(('SOLX', 'SOLY'))): warn_user("sunpy Map does not support three dimensional data " "and therefore cannot represent heliocentric coordinates. Proceed at your own risk.") if not all(su.is_equivalent(u.arcsec) for su in self.spatial_units): units = [su.to_string() for su in self.spatial_units] raise MapMetaValidationError( 'Map only supports spherical coordinate systems with angular units ' f'(ie. equivalent to arcsec), but this map has units {units}') # #### Data conversion routines #### #
[docs] def world_to_pixel(self, coordinate): """ Convert a world (data) coordinate to a pixel coordinate. Parameters ---------- coordinate : `~astropy.coordinates.SkyCoord` or `~astropy.coordinates.BaseCoordinateFrame` The coordinate object to convert to pixel coordinates. Returns ------- x : `~astropy.units.Quantity` Pixel coordinate on the CTYPE1 axis. y : `~astropy.units.Quantity` Pixel coordinate on the CTYPE2 axis. """ x, y = self.wcs.world_to_pixel(coordinate) return PixelPair(x * u.pixel, y * u.pixel)
[docs] @u.quantity_input def pixel_to_world(self, x: u.pixel, y: u.pixel): """ Convert a pixel coordinate to a data (world) coordinate. Parameters ---------- x : `~astropy.units.Quantity` Pixel coordinate of the CTYPE1 axis. (Normally solar-x). y : `~astropy.units.Quantity` Pixel coordinate of the CTYPE2 axis. (Normally solar-y). Returns ------- coord : `astropy.coordinates.SkyCoord` A coordinate object representing the output coordinate. """ return self.wcs.pixel_to_world(x, y)
# #### I/O routines #### #
[docs] def save(self, filepath, filetype='auto', **kwargs): """ Save a map to a file. Parameters ---------- filepath : `str` Location to save file to. filetype : `str`, optional Any supported file extension, defaults to ``"auto"``. hdu_type : `~astropy.io.fits.hdu.base.ExtensionHDU` instance or class, optional By default, a FITS file is written with the map in its primary HDU. If a type is given, a new HDU of this type will be created. If a HDU instance is given, its data and header will be updated from the map. Then that HDU instance will be written to the file. The example below uses `astropy.io.fits.CompImageHDU` to compress the map. kwargs : Any additional keyword arguments are passed to `~sunpy.io.write_file`. Notes ----- Saving with the jp2 extension will write a modified version of the given data casted to uint8 values in order to support the JPEG2000 format. Examples -------- >>> from astropy.io.fits import CompImageHDU >>> from sunpy.map import Map >>> import sunpy.data.sample # doctest: +REMOTE_DATA >>> aia_map = Map(sunpy.data.sample.AIA_171_IMAGE) # doctest: +REMOTE_DATA >>> aia_map.save("aia171.fits", hdu_type=CompImageHDU) # doctest: +REMOTE_DATA """ io.write_file(filepath, self.data, self.meta, filetype=filetype, **kwargs)
# #### Image processing routines #### #
[docs] @u.quantity_input def resample(self, dimensions: u.pixel, method='linear'): """ Resample to new dimension sizes. Uses the same parameters and creates the same coordinate lookup points as IDL''s congrid routine, which apparently originally came from a VAX/VMS routine of the same name. Parameters ---------- dimensions : `~astropy.units.Quantity` Output pixel dimensions. The first argument corresponds to the 'x' axis and the second argument corresponds to the 'y' axis. method : str Method to use for resampling interpolation. * ``'nearest'`` and ``'linear'`` - Use n x 1-D interpolations using `scipy.interpolate.interp1d`. * ``'spline'`` - Uses piecewise polynomials (splines) for mapping the input array to new coordinates by interpolation using `scipy.ndimage.map_coordinates`. Returns ------- out : `~sunpy.map.GenericMap` or subclass Resampled map References ---------- `Rebinning <https://scipy-cookbook.readthedocs.io/items/Rebinning.html>`_ """ # Note: because the underlying ndarray is transposed in sense when # compared to the Map, the ndarray is transposed, resampled, then # transposed back # Note: "center" defaults to True in this function because data # coordinates in a Map are at pixel centers # Make a copy of the original data and perform resample new_data = sunpy_image_resample(self.data.copy().T, dimensions, method, center=True) new_data = new_data.T scale_factor_x = float(self.dimensions[0] / dimensions[0]) scale_factor_y = float(self.dimensions[1] / dimensions[1]) # Update image scale and number of pixels new_meta = self.meta.copy() # Update metadata for key in {'cdelt1', 'cd1_1', 'cd2_1'} & self.meta.keys(): new_meta[key] *= scale_factor_x for key in {'cdelt2', 'cd1_2', 'cd2_2'} & self.meta.keys(): new_meta[key] *= scale_factor_y if 'pc1_1' in self.meta: new_meta['pc1_2'] *= scale_factor_y / scale_factor_x new_meta['pc2_1'] *= scale_factor_x / scale_factor_y new_meta['crpix1'] = (self.reference_pixel.x.to_value(u.pix) + 0.5) / scale_factor_x + 0.5 new_meta['crpix2'] = (self.reference_pixel.y.to_value(u.pix) + 0.5) / scale_factor_y + 0.5 new_meta['naxis1'] = new_data.shape[1] new_meta['naxis2'] = new_data.shape[0] # Create new map instance new_map = self._new_instance(new_data, new_meta, self.plot_settings) return new_map
[docs] @add_common_docstring(rotation_function_names=_rotation_function_names) @u.quantity_input def rotate(self, angle: u.deg = None, rmatrix=None, order=3, scale=1.0, recenter=False, missing=np.nan, *, method='scipy', clip=True): """ Returns a new rotated and rescaled map. Specify either a rotation angle or a rotation matrix, but not both. If neither an angle or a rotation matrix are specified, the map will be rotated by the rotation information in the metadata, which should derotate the map so that the pixel axes are aligned with world-coordinate axes. Parameters ---------- angle : `~astropy.units.Quantity` The angle (degrees) to rotate counterclockwise. rmatrix : array-like 2x2 linear transformation rotation matrix. order : int Interpolation order to be used. The precise meaning depends on the rotation method specified by ``method``. Default: 3 scale : float A scale factor for the image, default is no scaling recenter : bool If `True`, position the reference coordinate at the center of the new map Default: `False` missing : float The value to use for pixels in the output map that are beyond the extent of the input map. Default: `numpy.nan` method : {{{rotation_function_names}}}, optional Rotation function to use. Defaults to ``'scipy'``. clip : `bool`, optional If `True`, clips the pixel values of the output image to the range of the input image (including the value of ``missing``, if used). Defaults to `True`. Returns ------- out : `~sunpy.map.GenericMap` or subclass A new Map instance containing the rotated and rescaled data of the original map. See Also -------- sunpy.image.transform.affine_transform : The routine this method calls for the rotation. Notes ----- The rotation information in the metadata of the new map is modified appropriately from that of the original map to account for the applied rotation. It will solely be in the form of a PCi_j matrix, even if the original map used the CROTA2 keyword or a CDi_j matrix. If the map does not already contain pixels with `numpy.nan` values, setting ``missing`` to an appropriate number for the data (e.g., zero) will reduce the computation time. For each NaN pixel in the input image, one or more pixels in the output image will be set to NaN, with the size of the pixel region affected depending on the interpolation order. All currently implemented rotation methods require a convolution step to handle image NaNs. This convolution normally uses :func:`scipy.signal.convolve2d`, but if `OpenCV <https://opencv.org>`__ is installed, the faster |cv2_filter2D|_ is used instead. See :func:`sunpy.image.transform.affine_transform` for details on each of the rotation functions. """ if angle is not None and rmatrix is not None: raise ValueError("You cannot specify both an angle and a rotation matrix.") elif angle is None and rmatrix is None: # Be aware that self.rotation_matrix may not actually be a pure rotation matrix rmatrix = self.rotation_matrix if order not in range(6): raise ValueError("Order must be between 0 and 5.") method = _get_transform_method(method) # The FITS-WCS transform is by definition defined around the # reference coordinate in the header. lon, lat = self._get_lon_lat(self.reference_coordinate.frame) rotation_center = u.Quantity([lon, lat]) # Copy meta data new_meta = self.meta.copy() if angle is not None: # Calculate the parameters for the affine_transform c = np.cos(np.deg2rad(angle)) s = np.sin(np.deg2rad(angle)) rmatrix = np.array([[c, -s], [s, c]]) # The data will be rotated by the inverse of the rotation matrix. Because rmatrix may not # actually be a pure rotation matrix, we calculate the inverse in a general manner. inv_rmatrix = np.linalg.inv(rmatrix) # Calculate the shape in pixels to contain all of the image data corners = itertools.product([-0.5, self.data.shape[1]-0.5], [-0.5, self.data.shape[0]-0.5]) rot_corners = np.vstack([rmatrix @ c for c in corners]) * scale extent = np.max(rot_corners, axis=0) - np.min(rot_corners, axis=0) # Calculate the needed padding or unpadding diff = np.asarray(np.ceil((extent - np.flipud(self.data.shape)) / 2), dtype=int) pad_x = np.max((diff[0], 0)) pad_y = np.max((diff[1], 0)) unpad_x = -np.min((diff[0], 0)) unpad_y = -np.min((diff[1], 0)) # Raise an informative error message if trying to pad an integer array with NaNs if (pad_x > 0 or pad_y > 0) and issubclass(self.data.dtype.type, numbers.Integral) and (missing % 1 != 0): raise ValueError("The underlying data is integers, but the fill value for missing " "pixels cannot be cast to an integer, which is the case for the " "default fill value of NaN. Set the `missing` keyword to an " "appropriate integer value for the data set.") new_data = np.pad(self.data, ((pad_y, pad_y), (pad_x, pad_x)), mode='constant', constant_values=(missing, missing)) # All of the following pixel calculations use a pixel origin of 0 pixel_array_center = (np.flipud(new_data.shape) - 1) / 2.0 pixel_rotation_center = u.Quantity(self.reference_pixel).value + [pad_x, pad_y] if recenter: pixel_center = pixel_rotation_center else: pixel_center = pixel_array_center # Apply the rotation to the image data new_data = affine_transform(new_data, np.asarray(inv_rmatrix), order=order, scale=scale, image_center=pixel_center, recenter=recenter, missing=missing, method=method, clip=clip) if recenter: new_reference_pixel = pixel_array_center else: # Calculate new pixel coordinates for the rotation center new_reference_pixel = pixel_center + np.dot(rmatrix * scale, pixel_rotation_center - pixel_center) new_reference_pixel = np.array(new_reference_pixel).ravel() # Define the new reference_pixel new_meta['crval1'] = rotation_center[0].value new_meta['crval2'] = rotation_center[1].value new_meta['crpix1'] = new_reference_pixel[0] + 1 # FITS pixel origin is 1 new_meta['crpix2'] = new_reference_pixel[1] + 1 # FITS pixel origin is 1 new_meta['NAXIS1'] = new_data.shape[1] new_meta['NAXIS2'] = new_data.shape[0] # Unpad the array if necessary if unpad_x > 0: new_data = new_data[:, unpad_x:-unpad_x] new_meta['crpix1'] -= unpad_x if unpad_y > 0: new_data = new_data[unpad_y:-unpad_y, :] new_meta['crpix2'] -= unpad_y # Calculate the new rotation matrix to store in the header by # "subtracting" the rotation matrix used in the rotate from the old one # That being calculate the dot product of the old header data with the # inverse of the rotation matrix. pc_C = np.dot(self.rotation_matrix, inv_rmatrix) new_meta['PC1_1'] = pc_C[0, 0] new_meta['PC1_2'] = pc_C[0, 1] new_meta['PC2_1'] = pc_C[1, 0] new_meta['PC2_2'] = pc_C[1, 1] # Update pixel size if image has been scaled. if scale != 1.0: new_meta['cdelt1'] = (self.scale[0] / scale).value new_meta['cdelt2'] = (self.scale[1] / scale).value # Remove old CROTA kwargs because we have saved a new PCi_j matrix. new_meta.pop('CROTA1', None) new_meta.pop('CROTA2', None) # Remove CDi_j header new_meta.pop('CD1_1', None) new_meta.pop('CD1_2', None) new_meta.pop('CD2_1', None) new_meta.pop('CD2_2', None) # Create new map with the modification new_map = self._new_instance(new_data, new_meta, self.plot_settings) return new_map
[docs] @u.quantity_input def submap(self, bottom_left, *, top_right=None, width: (u.deg, u.pix) = None, height: (u.deg, u.pix) = None): """ Returns a submap defined by a rectangle. Any pixels which have at least part of their area inside the rectangle are returned. If the rectangle is defined in world coordinates, the smallest array which contains all four corners of the rectangle as defined in world coordinates is returned. Parameters ---------- bottom_left : `astropy.units.Quantity` or `~astropy.coordinates.SkyCoord` The bottom-left coordinate of the rectangle. If a `~astropy.coordinates.SkyCoord` it can have shape ``(2,)`` and simultaneously define ``top_right``. If specifying pixel coordinates it must be given as an `~astropy.units.Quantity` object with units of pixels. top_right : `astropy.units.Quantity` or `~astropy.coordinates.SkyCoord`, optional The top-right coordinate of the rectangle. If ``top_right`` is specified ``width`` and ``height`` must be omitted. width : `astropy.units.Quantity`, optional The width of the rectangle. Required if ``top_right`` is omitted. height : `astropy.units.Quantity` The height of the rectangle. Required if ``top_right`` is omitted. Returns ------- out : `~sunpy.map.GenericMap` or subclass A new map instance is returned representing to specified sub-region. Notes ----- When specifying pixel coordinates, they are specified in Cartesian order not in numpy order. So, for example, the ``bottom_left=`` argument should be ``[left, bottom]``. Examples -------- >>> import astropy.units as u >>> from astropy.coordinates import SkyCoord >>> import sunpy.map >>> import sunpy.data.sample # doctest: +REMOTE_DATA >>> aia = sunpy.map.Map(sunpy.data.sample.AIA_171_IMAGE) # doctest: +REMOTE_DATA >>> bl = SkyCoord(-300*u.arcsec, -300*u.arcsec, frame=aia.coordinate_frame) # doctest: +REMOTE_DATA >>> tr = SkyCoord(500*u.arcsec, 500*u.arcsec, frame=aia.coordinate_frame) # doctest: +REMOTE_DATA >>> aia.submap(bl, top_right=tr) # doctest: +REMOTE_DATA <sunpy.map.sources.sdo.AIAMap object at ...> SunPy Map --------- Observatory: SDO Instrument: AIA 3 Detector: AIA Measurement: 171.0 Angstrom Wavelength: 171.0 Angstrom Observation Date: 2011-06-07 06:33:02 Exposure Time: 0.234256 s Dimension: [335. 335.] pix Coordinate System: helioprojective Scale: [2.402792 2.402792] arcsec / pix Reference Pixel: [126.5 125.5] pix Reference Coord: [3.22309951 1.38578135] arcsec ... >>> aia.submap([0,0]*u.pixel, top_right=[5,5]*u.pixel) # doctest: +REMOTE_DATA <sunpy.map.sources.sdo.AIAMap object at ...> SunPy Map --------- Observatory: SDO Instrument: AIA 3 Detector: AIA Measurement: 171.0 Angstrom Wavelength: 171.0 Angstrom Observation Date: 2011-06-07 06:33:02 Exposure Time: 0.234256 s Dimension: [6. 6.] pix Coordinate System: helioprojective Scale: [2.402792 2.402792] arcsec / pix Reference Pixel: [511.5 511.5] pix Reference Coord: [3.22309951 1.38578135] arcsec ... >>> width = 10 * u.arcsec >>> height = 10 * u.arcsec >>> aia.submap(bl, width=width, height=height) # doctest: +REMOTE_DATA <sunpy.map.sources.sdo.AIAMap object at ...> SunPy Map --------- Observatory: SDO Instrument: AIA 3 Detector: AIA Measurement: 171.0 Angstrom Wavelength: 171.0 Angstrom Observation Date: 2011-06-07 06:33:02 Exposure Time: 0.234256 s Dimension: [5. 5.] pix Coordinate System: helioprojective Scale: [2.402792 2.402792] arcsec / pix Reference Pixel: [125.5 125.5] pix Reference Coord: [3.22309951 1.38578135] arcsec ... >>> bottom_left_vector = SkyCoord([0, 10] * u.deg, [0, 10] * u.deg, frame='heliographic_stonyhurst') >>> aia.submap(bottom_left_vector) # doctest: +REMOTE_DATA <sunpy.map.sources.sdo.AIAMap object at ...> SunPy Map --------- Observatory: SDO Instrument: AIA 3 Detector: AIA Measurement: 171.0 Angstrom Wavelength: 171.0 Angstrom Observation Date: 2011-06-07 06:33:02 Exposure Time: 0.234256 s Dimension: [70. 69.] pix Coordinate System: helioprojective Scale: [2.402792 2.402792] arcsec / pix Reference Pixel: [1.5 0.5] pix Reference Coord: [3.22309951 1.38578135] arcsec ... """ # Check that we have been given a valid combination of inputs # [False, False, False] is valid if bottom_left contains the two corner coords if ([arg is not None for arg in (top_right, width, height)] not in [[True, False, False], [False, False, False], [False, True, True]]): raise ValueError("Either top_right alone or both width and height must be specified.") # parse input arguments pixel_corners = u.Quantity(self._parse_submap_input( bottom_left, top_right, width, height)).T # The pixel corners result is in Cartesian order, so the first index is # columns and the second is rows. bottom = np.min(pixel_corners[1]).to_value(u.pix) top = np.max(pixel_corners[1]).to_value(u.pix) left = np.min(pixel_corners[0]).to_value(u.pix) right = np.max(pixel_corners[0]).to_value(u.pix) # Round the lower left pixel to the nearest integer # We want 0.5 to be rounded up to 1, so use floor(x + 0.5) bottom = np.floor(bottom + 0.5) left = np.floor(left + 0.5) # Round the top right pixel to the nearest integer, then add 1 for array indexing # We want e.g. 2.5 to be rounded down to 2, so use ceil(x - 0.5) top = np.ceil(top - 0.5) + 1 right = np.ceil(right - 0.5) + 1 # Clip pixel values to max of array, prevents negative # indexing bottom = int(np.clip(bottom, 0, self.data.shape[0])) top = int(np.clip(top, 0, self.data.shape[0])) left = int(np.clip(left, 0, self.data.shape[1])) right = int(np.clip(right, 0, self.data.shape[1])) arr_slice = np.s_[bottom:top, left:right] # Get ndarray representation of submap new_data = self.data[arr_slice].copy() # Make a copy of the header with updated centering information new_meta = self.meta.copy() # Add one to go from zero-based to one-based indexing new_meta['crpix1'] = self.reference_pixel.x.to_value(u.pix) + 1 - left new_meta['crpix2'] = self.reference_pixel.y.to_value(u.pix) + 1 - bottom new_meta['naxis1'] = new_data.shape[1] new_meta['naxis2'] = new_data.shape[0] # Create new map instance if self.mask is not None: new_mask = self.mask[arr_slice].copy() # Create new map with the modification new_map = self._new_instance(new_data, new_meta, self.plot_settings, mask=new_mask) return new_map # Create new map with the modification new_map = self._new_instance(new_data, new_meta, self.plot_settings) return new_map
@seconddispatch def _parse_submap_input(self, bottom_left, top_right, width, height): """ Should take any valid input to submap() and return bottom_left and top_right in pixel coordinates. """ @_parse_submap_input.register(u.Quantity) def _parse_submap_quantity_input(self, bottom_left, top_right, width, height): if top_right is None and width is None: raise ValueError('Either top_right alone or both width and height must be specified ' 'when bottom_left is a Quantity') if bottom_left.shape != (2, ): raise ValueError('bottom_left must have shape (2, ) when specified as a Quantity') if top_right is not None: if top_right.shape != (2, ): raise ValueError('top_right must have shape (2, ) when specified as a Quantity') if not top_right.unit.is_equivalent(u.pix): raise TypeError("When bottom_left is a Quantity, top_right " "must be a Quantity in units of pixels.") # Have bottom_left and top_right in pixels already, so no need to do # anything else else: if not (width.unit.is_equivalent(u.pix) and height.unit.is_equivalent(u.pix)): raise TypeError("When bottom_left is a Quantity, width and height " "must be a Quantity in units of pixels.") # Add width and height to get top_right top_right = u.Quantity([bottom_left[0] + width, bottom_left[1] + height]) top_left = u.Quantity([top_right[0], bottom_left[1]]) bottom_right = u.Quantity([bottom_left[0], top_right[1]]) return bottom_left, top_left, top_right, bottom_right @_parse_submap_input.register(SkyCoord) @_parse_submap_input.register(BaseCoordinateFrame) def _parse_submap_coord_input(self, bottom_left, top_right, width, height): # Use helper function to get top_right as a SkyCoord bottom_left, top_right = get_rectangle_coordinates(bottom_left, top_right=top_right, width=width, height=height) if isinstance(bottom_left, SkyCoord): frame = bottom_left.frame frame = bottom_left left_lon, bottom_lat = self._get_lon_lat(bottom_left) right_lon, top_lat = self._get_lon_lat(top_right) corners = SkyCoord([left_lon, left_lon, right_lon, right_lon], [bottom_lat, top_lat, top_lat, bottom_lat], frame=frame) return tuple(u.Quantity(self.wcs.world_to_pixel(corners), u.pix).T)
[docs] @u.quantity_input def superpixel(self, dimensions: u.pixel, offset: u.pixel = (0, 0)*u.pixel, func=np.sum): """Returns a new map consisting of superpixels formed by applying 'func' to the original map data. Parameters ---------- dimensions : tuple One superpixel in the new map is equal to (dimension[0], dimension[1]) pixels of the original map. The first argument corresponds to the 'x' axis and the second argument corresponds to the 'y' axis. If non-integer values are provided, they are rounded using `int`. offset : tuple Offset from (0,0) in original map pixels used to calculate where the data used to make the resulting superpixel map starts. If non-integer value are provided, they are rounded using `int`. func Function applied to the original data. The function 'func' must take a numpy array as its first argument, and support the axis keyword with the meaning of a numpy axis keyword (see the description of `~numpy.sum` for an example.) The default value of 'func' is `~numpy.sum`; using this causes superpixel to sum over (dimension[0], dimension[1]) pixels of the original map. Returns ------- out : `~sunpy.map.GenericMap` or subclass A new Map which has superpixels of the required size. References ---------- | `Summarizing blocks of an array using a moving window <https://mail.scipy.org/pipermail/numpy-discussion/2010-July/051760.html>`_ """ # Note: because the underlying ndarray is transposed in sense when # compared to the Map, the ndarray is transposed, resampled, then # transposed back. # Note: "center" defaults to True in this function because data # coordinates in a Map are at pixel centers. if (offset.value[0] < 0) or (offset.value[1] < 0): raise ValueError("Offset is strictly non-negative.") # These are rounded by int() in reshape_image_to_4d_superpixel, # so round here too for use in constructing metadata later. dimensions = [int(dim) for dim in dimensions.to_value(u.pix)] offset = [int(off) for off in offset.to_value(u.pix)] # Make a copy of the original data, perform reshaping, and apply the # function. if self.mask is not None: data = np.ma.array(self.data.copy(), mask=self.mask) else: data = self.data.copy() reshaped = reshape_image_to_4d_superpixel(data, [dimensions[1], dimensions[0]], [offset[1], offset[0]]) new_array = func(func(reshaped, axis=3), axis=1) # Update image scale and number of pixels # create copy of new meta data new_meta = self.meta.copy() # Update metadata for key in {'cdelt1', 'cd1_1', 'cd2_1'} & self.meta.keys(): new_meta[key] *= dimensions[0] for key in {'cdelt2', 'cd1_2', 'cd2_2'} & self.meta.keys(): new_meta[key] *= dimensions[1] if 'pc1_1' in self.meta: new_meta['pc1_2'] *= dimensions[1] / dimensions[0] new_meta['pc2_1'] *= dimensions[0] / dimensions[1] new_meta['crpix1'] = ((self.reference_pixel.x.to_value(u.pix) + 0.5 - offset[0]) / dimensions[0]) + 0.5 new_meta['crpix2'] = ((self.reference_pixel.y.to_value(u.pix) + 0.5 - offset[1]) / dimensions[1]) + 0.5 new_meta['naxis1'] = new_array.shape[1] new_meta['naxis2'] = new_array.shape[0] # Create new map instance if self.mask is not None: new_data = np.ma.getdata(new_array) new_mask = np.ma.getmask(new_array) else: new_data = new_array new_mask = None # Create new map with the modified data new_map = self._new_instance(new_data, new_meta, self.plot_settings, mask=new_mask) return new_map
# #### Visualization #### # @property def cmap(self): """ Return the `matplotlib.colors.Colormap` instance this map uses. """ cmap = self.plot_settings['cmap'] if isinstance(cmap, str): cmap = plt.get_cmap(cmap) # Set the colormap to be this specific instance so we are not # returning a copy self.plot_settings['cmap'] = cmap return cmap
[docs] @u.quantity_input def draw_grid(self, axes=None, grid_spacing: u.deg = 15*u.deg, annotate=True, system='stonyhurst', **kwargs): """ Draws a coordinate overlay on the plot in the Heliographic Stonyhurst coordinate system. To overlay other coordinate systems see the `WCSAxes Documentation <https://docs.astropy.org/en/stable/visualization/wcsaxes/overlaying_coordinate_systems.html>`_ Parameters ---------- axes : `~matplotlib.axes` or `None` Axes to plot limb on, or `None` to use current axes. grid_spacing : `~astropy.units.Quantity` Spacing for longitude and latitude grid, if length two it specifies (lon, lat) spacing. annotate : `bool` Passing `False` disables the axes labels and the ticks on the top and right axes. system : str Coordinate system for the grid. Must be 'stonyhurst' or 'carrington'. kwargs : Additional keyword arguments are passed to `~sunpy.visualization.wcsaxes_compat.wcsaxes_heliographic_overlay`. Returns ------- overlay: `~astropy.visualization.wcsaxes.CoordinatesMap` The wcsaxes coordinate overlay instance. Notes ----- Keyword arguments are passed onto the `sunpy.visualization.wcsaxes_compat.wcsaxes_heliographic_overlay` function. """ axes = self._check_axes(axes) return wcsaxes_compat.wcsaxes_heliographic_overlay(axes, grid_spacing=grid_spacing, annotate=annotate, obstime=self.date, rsun=self.rsun_meters, observer=self.observer_coordinate, system=system, **kwargs)
[docs] def draw_limb(self, axes=None, *, resolution=1000, **kwargs): """ Draws the solar limb as seen by the map's observer. The limb is a circle for only the simplest plots. If the coordinate frame of the limb is different from 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 or ``None`` to use current axes. resolution : `int` The number of points to use to represent the limb. Returns ------- visible : `~matplotlib.patches.Polygon` or `~matplotlib.patches.Circle` 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 visible points (e.g., on a synoptic map any limb is fully) visible ``hidden`` 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. """ # Put imports here to reduce sunpy.map import time import sunpy.visualization.drawing axes = self._check_axes(axes) return sunpy.visualization.drawing.limb( axes, self.observer_coordinate, resolution=resolution, rsun=self.rsun_meters, **kwargs )
[docs] @u.quantity_input def draw_quadrangle(self, bottom_left, *, width: (u.deg, u.pix) = None, height: (u.deg, u.pix) = None, axes=None, top_right=None, **kwargs): """ Draw a quadrangle defined in world coordinates on the plot using Astropy's `~astropy.visualization.wcsaxes.Quadrangle`. This draws a quadrangle that has corners at ``(bottom_left, top_right)``, and has sides aligned with the coordinate axes of the frame of ``bottom_left``, which may be different from the coordinate axes of the map. If ``width`` and ``height`` are specified, they are respectively added to the longitude and latitude of the ``bottom_left`` coordinate to calculate a ``top_right`` coordinate. Parameters ---------- bottom_left : `~astropy.coordinates.SkyCoord` or `~astropy.units.Quantity` The bottom-left coordinate of the rectangle. If a `~astropy.coordinates.SkyCoord` it can have shape ``(2,)`` and simultaneously define ``top_right``. If specifying pixel coordinates it must be given as an `~astropy.units.Quantity` object with pixel units (e.g., ``pix``). top_right : `~astropy.coordinates.SkyCoord` or `~astropy.units.Quantity`, optional The top-right coordinate of the quadrangle. If ``top_right`` is specified ``width`` and ``height`` must be omitted. width : `astropy.units.Quantity`, optional The width of the quadrangle. Required if ``top_right`` is omitted. height : `astropy.units.Quantity` The height of the quadrangle. Required if ``top_right`` is omitted. axes : `matplotlib.axes.Axes` The axes on which to plot the quadrangle. Defaults to the current axes. Returns ------- quad : `~astropy.visualization.wcsaxes.Quadrangle` The added patch. Notes ----- Extra keyword arguments to this function are passed through to the `~astropy.visualization.wcsaxes.Quadrangle` instance. Examples -------- .. minigallery:: sunpy.map.GenericMap.draw_quadrangle """ axes = self._check_axes(axes) if isinstance(bottom_left, u.Quantity): anchor, _, top_right, _ = self._parse_submap_quantity_input(bottom_left, top_right, width, height) width, height = top_right - anchor transform = axes.get_transform(self.wcs if self.wcs is not axes.wcs else 'pixel') kwargs.update({"vertex_unit": u.pix}) else: bottom_left, top_right = get_rectangle_coordinates( bottom_left, top_right=top_right, width=width, height=height) width = Longitude(top_right.spherical.lon - bottom_left.spherical.lon) height = top_right.spherical.lat - bottom_left.spherical.lat anchor = self._get_lon_lat(bottom_left) transform = axes.get_transform(bottom_left.replicate_without_data()) kwergs = { "transform": transform, "edgecolor": "white", "fill": False, } kwergs.update(kwargs) quad = Quadrangle(anchor, width, height, **kwergs) axes.add_patch(quad) return quad
def _process_levels_arg(self, levels): """ Accept a percentage or dimensionless or map unit input for contours. """ levels = np.atleast_1d(levels) if not hasattr(levels, 'unit'): if self.unit is None: # No map units, so allow non-quantity through return levels else: raise TypeError("The levels argument has no unit attribute, " "it should be an Astropy Quantity object.") if levels.unit == u.percent: return 0.01 * levels.to_value('percent') * np.nanmax(self.data) elif self.unit is not None: return levels.to_value(self.unit) elif levels.unit.is_equivalent(u.dimensionless_unscaled): # Handle case where map data has no units return levels.to_value(u.dimensionless_unscaled) else: # Map data has no units, but levels doesn't have dimensionless units raise u.UnitsError("This map has no unit, so levels can only be specified in percent " "or in u.dimensionless_unscaled units.")
[docs] def draw_contours(self, levels, axes=None, *, fill=False, **contour_args): """ Draw contours of the data. Parameters ---------- levels : `~astropy.units.Quantity` A list of numbers indicating the contours to draw. These are given as a percentage of the maximum value of the map data, or in units equivalent to the `~sunpy.map.GenericMap.unit` attribute. axes : `matplotlib.axes.Axes` The axes on which to plot the contours. Defaults to the current axes. fill : `bool`, optional Determines the style of the contours: - If `False` (default), contours are drawn as lines using :meth:`~matplotlib.axes.Axes.contour`. - If `True`, contours are drawn as filled regions using :meth:`~matplotlib.axes.Axes.contourf`. Returns ------- cs : `list` The `~matplotlib.contour.QuadContourSet` object, after it has been added to ``axes``. Notes ----- Extra keyword arguments to this function are passed through to the corresponding matplotlib method. """ axes = self._check_axes(axes) levels = self._process_levels_arg(levels) # Pixel indices y, x = np.indices(self.data.shape) # Prepare a local variable in case we need to mask values data = self.data # Transform the indices if plotting to a different WCS # We do this instead of using the `transform` keyword argument so that Matplotlib does not # get confused about the bounds of the contours if self.wcs is not axes.wcs: transform = axes.get_transform(self.wcs) - axes.transData # pixel->pixel transform x_1d, y_1d = transform.transform(np.stack([x.ravel(), y.ravel()]).T).T x, y = np.reshape(x_1d, x.shape), np.reshape(y_1d, y.shape) # Mask out the data array anywhere the coordinate arrays are not finite data = np.ma.array(data, mask=~np.logical_and(np.isfinite(x), np.isfinite(y))) if fill: # Ensure we have more than one level if fill is True if len(levels) == 1: max_val = np.nanmax(self.data) # Ensure the existing level is less than max_val if levels[0] < max_val: levels = np.append(levels, max_val) else: raise ValueError( f"The provided level ({levels[0]}) is not smaller than the maximum data value ({max_val}). " "Contour level must be smaller than the maximum data value to use `fill=True`.") cs = axes.contourf(x, y, data, levels, **contour_args) else: cs = axes.contour(x, y, data, levels, **contour_args) return cs
[docs] @peek_show def peek(self, draw_limb=False, draw_grid=False, colorbar=True, **matplot_args): """ Displays a graphical overview of the data in this object for user evaluation. For the creation of plots, users should instead use the `~sunpy.map.GenericMap.plot` method and Matplotlib's pyplot framework. Parameters ---------- draw_limb : bool Whether the solar limb should be plotted. draw_grid : bool or `~astropy.units.Quantity` Whether solar meridians and parallels are plotted. If `~astropy.units.Quantity` then sets degree difference between parallels and meridians. colorbar : bool Whether to display a colorbar next to the plot. **matplot_args : dict Matplotlib Any additional imshow arguments that should be used when plotting. """ figure = plt.figure() axes = wcsaxes_compat.gca_wcs(self.wcs) im = self.plot(axes=axes, **matplot_args) grid_spacing = None # Handle case where draw_grid is actually the grid sapcing if isinstance(draw_grid, u.Quantity): grid_spacing = draw_grid draw_grid = True elif not isinstance(draw_grid, bool): raise TypeError("draw_grid should be a bool or an astropy Quantity.") if colorbar: if draw_grid: pad = 0.12 # Pad to compensate for ticks and axes labels else: pad = 0.05 # Default value for vertical colorbar colorbar_label = str(self.unit) if self.unit is not None else "" figure.colorbar(im, pad=pad).set_label(colorbar_label, rotation=0, labelpad=-50, y=-0.02, size=12) if draw_limb: self.draw_limb(axes=axes) if draw_grid: if grid_spacing is None: self.draw_grid(axes=axes) else: self.draw_grid(axes=axes, grid_spacing=grid_spacing) return figure
[docs] @u.quantity_input def plot(self, annotate=True, axes=None, title=True, autoalign=False, clip_interval: u.percent = None, **imshow_kwargs): """ Plots the map object using matplotlib, in a method equivalent to :meth:`~matplotlib.axes.Axes.imshow` using nearest neighbor interpolation. Parameters ---------- annotate : `bool`, optional If `True`, the data is plotted at its natural scale; with title and axis labels. axes : `~matplotlib.axes.Axes` or None If provided the image will be plotted on the given axes. Else the current Matplotlib axes will be used. title : `str`, `bool`, optional The plot title. If `True`, uses the default title for this map. clip_interval : two-element `~astropy.units.Quantity`, optional If provided, the data will be clipped to the percentile interval bounded by the two numbers. autoalign : `bool` or `str`, optional If other than `False`, the plotting accounts for any difference between the WCS of the map and the WCS of the `~astropy.visualization.wcsaxes.WCSAxes` axes (e.g., a difference in rotation angle). If ``pcolormesh``, this method will use :meth:`~matplotlib.axes.Axes.pcolormesh` instead of the default :meth:`~matplotlib.axes.Axes.imshow`. Specifying `True` is equivalent to specifying ``pcolormesh``. **imshow_kwargs : `dict` Any additional imshow arguments are passed to :meth:`~matplotlib.axes.Axes.imshow`. Examples -------- >>> # Simple Plot with color bar >>> aia.plot() # doctest: +SKIP >>> plt.colorbar() # doctest: +SKIP >>> # Add a limb line and grid >>> aia.plot() # doctest: +SKIP >>> aia.draw_limb() # doctest: +SKIP >>> aia.draw_grid() # doctest: +SKIP Notes ----- The ``autoalign`` functionality is computationally intensive. If the plot will be interactive, the alternative approach of preprocessing the map (e.g., de-rotating it) to match the desired axes will result in better performance. When combining ``autoalign`` functionality with `~sunpy.coordinates.Helioprojective` coordinates, portions of the map that are beyond the solar disk may not appear, which may also inhibit Matplotlib's autoscaling of the plot limits. The plot limits can be set manually. To preserve the off-disk parts of the map, using the :meth:`~sunpy.coordinates.Helioprojective.assume_spherical_screen` context manager may be appropriate. """ # Users sometimes assume that the first argument is `axes` instead of `annotate` if not isinstance(annotate, bool): raise TypeError("You have provided a non-boolean value for the `annotate` parameter. " "If you are specifying the axes, use `axes=...` to pass it in.") # Set the default approach to autoalignment if autoalign not in [False, True, 'pcolormesh']: raise ValueError("The value for `autoalign` must be False, True, or 'pcolormesh'.") if autoalign is True: autoalign = 'pcolormesh' axes = self._check_axes(axes, warn_different_wcs=autoalign is False) # Normal plot plot_settings = copy.deepcopy(self.plot_settings) if 'title' in plot_settings: plot_settings_title = plot_settings.pop('title') else: plot_settings_title = self.latex_name # Anything left in plot_settings is given to imshow imshow_args = plot_settings if annotate: if title is True: title = plot_settings_title if title: axes.set_title(title) # WCSAxes has unit identifiers on the tick labels, so no need # to add unit information to the label spatial_units = [None, None] ctype = axes.wcs.wcs.ctype axes.set_xlabel(axis_labels_from_ctype(ctype[0], spatial_units[0])) axes.set_ylabel(axis_labels_from_ctype(ctype[1], spatial_units[1])) # Take a deep copy here so that a norm in imshow_kwargs doesn't get modified # by setting it's vmin and vmax imshow_args.update(copy.deepcopy(imshow_kwargs)) if clip_interval is not None: if len(clip_interval) == 2: clip_percentages = clip_interval.to('%').value vmin, vmax = AsymmetricPercentileInterval(*clip_percentages).get_limits(self.data) else: raise ValueError("Clip percentile interval must be specified as two numbers.") imshow_args['vmin'] = vmin imshow_args['vmax'] = vmax msg = ('Cannot manually specify {0}, as the norm ' 'already has {0} set. To prevent this error set {0} on ' '`m.plot_settings["norm"]` or the norm passed to `m.plot`.') if (norm := imshow_args.get('norm', None)) is not None: if 'vmin' in imshow_args: if norm.vmin is not None: raise ValueError(msg.format('vmin')) norm.vmin = imshow_args.pop('vmin') if 'vmax' in imshow_args: if norm.vmax is not None: raise ValueError(msg.format('vmax')) norm.vmax = imshow_args.pop('vmax') if self.mask is None: data = self.data else: data = np.ma.array(np.asarray(self.data), mask=self.mask) if autoalign == 'pcolormesh': # We have to handle an `aspect` keyword separately axes.set_aspect(imshow_args.get('aspect', 1)) # pcolormesh does not do interpolation if imshow_args.get('interpolation', None) not in [None, 'none', 'nearest']: warn_user("The interpolation keyword argument is ignored when using autoalign " "functionality.") # Remove imshow keyword arguments that are not accepted by pcolormesh for item in ['aspect', 'extent', 'interpolation', 'origin']: if item in imshow_args: del imshow_args[item] imshow_args.setdefault('transform', axes.get_transform(self.wcs)) # The quadrilaterals of pcolormesh can slightly overlap, which creates the appearance # of a grid pattern when alpha is not 1. These settings minimize the overlap. if imshow_args.get('alpha', 1) != 1: imshow_args.setdefault('antialiased', True) imshow_args.setdefault('linewidth', 0) ret = axes.pcolormesh(np.arange(data.shape[1] + 1) - 0.5, np.arange(data.shape[0] + 1) - 0.5, data, **imshow_args) else: ret = axes.imshow(data, **imshow_args) wcsaxes_compat.default_wcs_grid(axes) # Set current axes/image if pyplot is being used (makes colorbar work) for i in plt.get_fignums(): if axes in plt.figure(i).axes: plt.sca(axes) plt.sci(ret) return ret
[docs] def contour(self, level, **kwargs): """ Returns coordinates of the contours for a given level value. For details of the contouring algorithm see `skimage.measure.find_contours`. Parameters ---------- level : float, astropy.units.Quantity Value along which to find contours in the array. If the map unit attribute is not `None`, this must be a `~astropy.units.Quantity` with units equivalent to the map data units. kwargs : Additional keyword arguments are passed to `skimage.measure.find_contours`. Returns ------- contours: list of (n,2) `~astropy.coordinates.SkyCoord` Coordinates of each contour. Examples -------- >>> import astropy.units as u >>> import sunpy.map >>> import sunpy.data.sample # doctest: +REMOTE_DATA >>> aia = sunpy.map.Map(sunpy.data.sample.AIA_171_IMAGE) # doctest: +REMOTE_DATA >>> contours = aia.contour(50000 * u.ct) # doctest: +REMOTE_DATA >>> print(contours[0]) # doctest: +REMOTE_DATA <SkyCoord (Helioprojective: obstime=2011-06-07T06:33:02.770, rsun=696000.0 km, observer=<HeliographicStonyhurst Coordinate (obstime=2011-06-07T06:33:02.770, rsun=696000.0 km): (lon, lat, radius) in (deg, deg, m) (-0.00406308, 0.04787238, 1.51846026e+11)>): (Tx, Ty) in arcsec [(719.59798458, -352.60839064), (717.19243987, -353.75348121), ... See Also -------- skimage.measure.find_contours """ from skimage import measure level = self._process_levels_arg(level) if level.size != 1: raise ValueError("level must be a single scalar value") else: # _process_levels_arg converts level to a 1D array, but # find_contours expects a scalar below level = level[0] contours = measure.find_contours(self.data, level=level, **kwargs) contours = [self.wcs.array_index_to_world(c[:, 0], c[:, 1]) for c in contours] return contours
def _check_axes(self, axes, warn_different_wcs=False): """ - If axes is None, get the current Axes object. - Error if not a WCSAxes. - Return axes. Parameters ---------- axes : matplotlib.axes.Axes Axes to validate. warn_different_wcs : bool If `True`, warn if the Axes WCS is different from the Map WCS. This is only used for `.plot()`, and can be removed once support is added for plotting a map on a different WCSAxes. """ if not axes: axes = wcsaxes_compat.gca_wcs(self.wcs) if not wcsaxes_compat.is_wcsaxes(axes): raise TypeError("The axes need to be an instance of WCSAxes. " "To fix this pass set the `projection` keyword " "to this map when creating the axes.") elif warn_different_wcs and not axes.wcs.wcs.compare(self.wcs.wcs, tolerance=0.01): warn_user('The map world coordinate system (WCS) is different from the axes WCS. ' 'The map data axes may not correctly align with the coordinate axes. ' 'To automatically transform the data to the coordinate axes, specify ' '`autoalign=True`.') return axes
[docs] def reproject_to(self, target_wcs, *, algorithm='interpolation', return_footprint=False, **reproject_args): """ Reproject the map to a different world coordinate system (WCS) .. note:: This method requires the optional package `reproject` to be installed. Additional keyword arguments are passed through to the reprojection function. Parameters ---------- target_wcs : `dict` or `~astropy.wcs.WCS` The destination FITS WCS header or WCS instance algorithm : `str` One of the supported `reproject` algorithms (see below) return_footprint : `bool` If ``True``, the footprint is returned in addition to the new map. Defaults to ``False``. Returns ------- outmap : `~sunpy.map.GenericMap` The reprojected map footprint : `~numpy.ndarray` Footprint of the input arary in the output array. Values of 0 indicate no coverage or valid values in the input image, while values of 1 indicate valid values. Intermediate values indicate partial coverage. Only returned if ``return_footprint`` is ``True``. Notes ----- The reprojected map does not preserve any metadata beyond the WCS-associated metadata. The supported `reproject` algorithms are: * 'interpolation' for :func:`~reproject.reproject_interp` * 'adaptive' for :func:`~reproject.reproject_adaptive` * 'exact' for :func:`~reproject.reproject_exact` See the respective documentation for these functions for additional keyword arguments that are allowed. .. minigallery:: sunpy.map.GenericMap.reproject_to """ try: import reproject except ImportError as exc: raise ImportError("This method requires the optional package `reproject`.") from exc if not isinstance(target_wcs, astropy.wcs.WCS): target_wcs = astropy.wcs.WCS(target_wcs) # Select the desired reprojection algorithm functions = {'interpolation': reproject.reproject_interp, 'adaptive': reproject.reproject_adaptive, 'exact': reproject.reproject_exact} if algorithm not in functions: raise ValueError(f"The specified algorithm must be one of: {list(functions.keys())}") func = functions[algorithm] # reproject does not automatically grab the array shape from the WCS instance if target_wcs.array_shape is not None: reproject_args.setdefault('shape_out', target_wcs.array_shape) # Reproject the array output_array = func(self, target_wcs, return_footprint=return_footprint, **reproject_args) if return_footprint: output_array, footprint = output_array # Create and return a new GenericMap outmap = GenericMap(output_array, target_wcs.to_header(), plot_settings=self.plot_settings) # Check rsun mismatch if self.rsun_meters != outmap.rsun_meters: warn_user("rsun mismatch detected: " f"{self.name}.rsun_meters={self.rsun_meters}; {outmap.name}.rsun_meters={outmap.rsun_meters}. " "This might cause unexpected results during reprojection.") if return_footprint: return outmap, footprint return outmap
GenericMap.__doc__ += textwrap.indent(_notes_doc, " ") class InvalidHeaderInformation(ValueError): """Exception to raise when an invalid header tag value is encountered for a FITS/JPEG 2000 file."""