Overlaying Two Maps

This example demonstrates how to draw contours of one map on top of another to compare features.

import matplotlib.colors
import matplotlib.pyplot as plt

import astropy.units as u
from astropy.coordinates import SkyCoord

import sunpy.data.sample
from sunpy.map import Map

We start with the sample data. HMI shows the line-of-sight magnetic field at the photosphere while AIA 171 images show the resulting magnetic fields filled with hot plasma above, in the corona. We want to see what coronal features overlap with regions of strong line-of-sight magnetic fields.

aia_map = Map(sunpy.data.sample.AIA_171_IMAGE)
hmi_map = Map(sunpy.data.sample.HMI_LOS_IMAGE)

bottom_left = [0, 0] * u.arcsec
top_right = [800, 800] * u.arcsec
aia_smap = aia_map.submap(SkyCoord(*bottom_left,
                                   frame=aia_map.coordinate_frame),
                          top_right=SkyCoord(*top_right,
                                             frame=aia_map.coordinate_frame))
hmi_smap = hmi_map.submap(SkyCoord(*bottom_left,
                                   frame=hmi_map.coordinate_frame),
                          top_right=SkyCoord(*top_right,
                                             frame=hmi_map.coordinate_frame))

Let’s set the contours of the HMI map from a few hundred to a thousand Gauss, which is the typical field strength associated with umbral regions of pores and sunspots.

levels = [-1000, -500, -250, 250, 500, 1000] * u.G

Now let’s look at the result. Notice that we can see the coronal structures present on the AIA image and how they correspond to the line of sight magnetic field.

fig = plt.figure()
ax = fig.add_subplot(projection=aia_smap)
aia_smap.plot(axes=ax)
hmi_smap.draw_contours(axes=ax, levels=levels, cmap="viridis")
norm = matplotlib.colors.Normalize(vmin=-1000, vmax=1000)
cbar = fig.colorbar(plt.cm.ScalarMappable(norm=norm, cmap="viridis"), ax=ax)
cbar.ax.set_title("Gauss", fontsize=10)
ax.set_title("AIA & HMI Composite Plot")

plt.show()