Plotting in SunPy¶
SunPy makes use of matplotlib for all of its plotting - as such, it tries to follow the matplotlib plotting philosophy. It is therefore useful to go over how matplotlib works as background.
1. Matplotlib Tutorial¶
The tutorial provided here is a summary of one that can be found in the matplotlib usage documentation.
Matplotlib provides two main pathways for plotting. One is meant for interactive use (e.g. command-line) and the other for non-interactive use (e.g. modules). It is important to recognize though that the interactive-use pathway (referred to as pyplot) just provides shortcuts for doing many of the more advanced non-interactive functions in the background. It is therefore possible to switch between the two as necessary and it is possible to use pyplot in a non-interactive way. In this manner pyplot is just a shortcut to making it quicker to set up plot axes and figures. In order to get access to the full interactive capabilities of pyplot it is necessary to turn this feature on. Pylab is another matplotlib usage scenario but it is essentially just pyplot with the interactive capabilities turned on and numpy and matplotlib imported into the main namespace.
Here is a simple example of pyplot usage.
import matplotlib.pyplot as plt plt.plot(range(10), range(10)) plt.title("A simple Plot") plt.show()
show command opens a plot on the screen and blocks
execution until the plot window is closed. The
command only works once. If you were to call
after the above code is executed nothing happens. This confusing behavior
is something that the matplotlib devs get complaints about often and so this may change.
A discussion about this can be found here.
Don’t be confused by another command called
This is only used while in interactive mode.
To turn on interactivity for pyplot use the command
In interactive mode, the plot will appear at the first
command and most commands will update the plot as you call them. Here is some
>>> plt.plot(range(10), range(10)) >>> plt.title("Simple Plot")
In this example, you’ll see that the title appears right on the plot when you call it.
Note that in this case the
show command is useless as the
plot shows up right when you create it. Also note that some commands will not
automatically update the plot and you have to use the
command. The following command
turns off interactivity.
3. Advanced Pyplot¶
If you need more fine-grained control over plots the recommended path is to use pyplot and access the figures and axes objects. This is shown in the following example.
import matplotlib.pyplot as plt import numpy as np x = np.arange(0, 10, 0.2) y = np.sin(x) fig = plt.figure() ax = plt.subplot() ax.plot(x, y) ax.set_xlabel('x') plt.show()
Figure is the top-level container for all plot elements and
Axes is the top-level container for a particular plot. So the above example,
creates a figure then creates an axes and populates the plot in
ax. With this method you
now have your hands on the
Axes object so you can do things
like change the labels on the x and y axes or add a legend.
In the previous section, pyplot took care of creating these
objects for you so you didn’t have to worry about creating them yourself.
4. SunPy Plotting Convention¶
To be consistent with matplotlib, SunPy has developed a standard plotting policy which supports both simple and advanced matplotlib usage. The following examples focus on the map object but they should be applicable across all of the data objects.
For quick and easy access to a plot all SunPy base objects (i.e. maps, spectra,
timeseries) define their own
peek command which
will create a plot for you and show it without you having to deal with any
matplotlib setup. This is so that it is easy to take a quick look at your data.
For example you can make the following plot.
import sunpy.map import sunpy.data.sample smap = sunpy.map.Map(sunpy.data.sample.AIA_171_IMAGE) smap.peek(draw_limb=True)
This creates a plot window with all axes defined, a plot title, and the image of the map data defined by the contents of the map. In non-interactive mode the plot window blocks the command line terminal and must be closed before doing anything else.
For more advanced plotting the base SunPy objects also provide a
plot command. This command is similar to the
imshow command in that it will create a figure and
axes object for you if you haven’t already.
plot it is possible to customise the look of the
plot by combining SunPy and matplotlib commands, for example you can over plot
contours on the Map:
import matplotlib.pyplot as plt import astropy.units as u import sunpy.map import sunpy.data.sample aia_map = sunpy.map.Map(sunpy.data.sample.AIA_171_IMAGE) aia_map.plot() aia_map.draw_limb() # let's add contours as well aia_map.draw_contours([10,20,30,40,50,60,70,80,90] * u.percent) plt.colorbar() plt.show()
import matplotlib.pyplot as plt import astropy.units as u from astropy.coordinates import SkyCoord import sunpy.map import sunpy.data.sample smap = sunpy.map.Map(sunpy.data.sample.AIA_171_IMAGE) fig = plt.figure() # Provide the Map as a projection, which creates a WCSAxes object ax = plt.subplot(projection=smap) im = smap.plot() # Prevent the image from being re-scaled while overplotting. ax.set_autoscale_on(False) xc = [0,100,1000] * u.arcsec yc = [0,100,1000] * u.arcsec coords = SkyCoord(xc, yc, frame=smap.coordinate_frame) p = ax.plot_coord(coords, 'o') # Set title. ax.set_title('Custom plot with WCSAxes') plt.colorbar() plt.show()
It is possible to create the same plot, explicitly not using
wcsaxes, however, this will not have the features of
wcsaxes which include correct representation of
rotation and plotting in different coordinate systems.
import matplotlib.pyplot as plt import astropy.units as u import sunpy.map import sunpy.data.sample smap = sunpy.map.Map(sunpy.data.sample.AIA_171_IMAGE) fig, ax = plt.subplots() im = smap.plot() # Prevent the image from being re-scaled while overplotting. ax.set_autoscale_on(False) xc = [0,100,1000] * u.arcsec yc = [0,100,1000] * u.arcsec p = plt.plot(xc, yc, 'o') # Set title. ax.set_title('Custom plot without WCSAxes') plt.colorbar() plt.show()
Plotting Maps with wcsaxes¶
By default SunPy map uses the
astropy.visualization.wcsaxes module to improve
the representation of world coordinates, and calling
peek() will use wcsaxes
for plotting. Unless a standard
matplotlib.axes.Axes object is explicitly
To explicitly create a
WCSAxes instance do the
>>> fig = plt.figure() >>> ax = plt.subplot(projection=smap)
when plotting on an
WCSAxes axes, it will by
default plot in pixel coordinates, you can override this behavior and plot in
‘world’ coordinates by getting the transformation from the axes with
ax.get_transform('world'). Note: World coordinates are always in degrees
so you will have to convert to degrees.:
>>> smap.plot() >>> ax.plot((100*u.arcsec).to(u.deg), (500*u.arcsec).to(u.deg), ... transform=ax.get_transform('world'))
Finally, here is a more complex example using SunPy maps, wcsaxes and Astropy units to plot a AIA image and a zoomed in view of an active region.
import matplotlib.pyplot as plt from matplotlib import patches import astropy.units as u from astropy.coordinates import SkyCoord import sunpy.map import sunpy.data.sample # Define a region of interest length = 250 * u.arcsec x0 = -100 * u.arcsec y0 = -400 * u.arcsec # Create a SunPy Map, and a second submap over the region of interest. smap = sunpy.map.Map(sunpy.data.sample.AIA_171_IMAGE) bottom_left = SkyCoord(x0 - length, y0 - length, frame=smap.coordinate_frame) top_right = SkyCoord(x0 + length, y0 + length, frame=smap.coordinate_frame) submap = smap.submap(bottom_left, top_right=top_right) # Create a new matplotlib figure, larger than default. fig = plt.figure(figsize=(5, 12)) # Add a first Axis, using the WCS from the map. ax1 = fig.add_subplot(2, 1, 1, projection=smap) # Plot the Map on the axes with default settings. smap.plot() # Draw a box on the image smap.draw_rectangle(bottom_left, height=length * 2, width=length * 2) # Create a second axis on the plot. ax2 = fig.add_subplot(2, 1, 2, projection=submap) submap.plot() # Add a overlay grid. submap.draw_grid(grid_spacing=10*u.deg) # Change the title. ax2.set_title('Zoomed View') plt.show()