Aligning AIA and HMI Data with Reproject

A common case when using data from multiple sources, is aligning so they have the same reference frame. It is also possible to use reproject to align data, by reprojecting one image to the WCS of another. This is a very generic way of aligning data, and can be very accurate.

You will need reproject v0.6 or higher installed.

import matplotlib.pyplot as plt
from reproject import reproject_interp

import astropy.units as u

import sunpy.map
from sunpy.net import Fido
from sunpy.net import attrs as a

In this example we are going to make a lot of side by side figures, so let’s change the default figure size.

plt.rcParams['figure.figsize'] = (16, 8)

We are going to download one AIA and one HMI magnetogram image.

time = (a.vso.Sample(24 * u.hour) &
        a.Time('2010-08-19', '2010-08-19T00:10:00', '2010-08-19') &
        a.vso.Extent(0, 0, 0, 0, "FULLDISK"))
aia = a.Instrument('AIA') & a.Wavelength(17 * u.nm, 18 * u.nm)
hmi = a.Instrument('HMI') & a.vso.Physobs("LOS_magnetic_field")


res = Fido.search(time, aia | hmi)

files = Fido.fetch(res[:, 0])

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/parfive/downloader.py:86: DeprecationWarning: The loop argument is deprecated since Python 3.8, and scheduled for removal in Python 3.10.
  self.http_queue = asyncio.Queue(loop=self.loop)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/asyncio/queues.py:48: DeprecationWarning: The loop argument is deprecated since Python 3.8, and scheduled for removal in Python 3.10.
  self._finished = locks.Event(loop=loop)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/parfive/downloader.py:87: DeprecationWarning: The loop argument is deprecated since Python 3.8, and scheduled for removal in Python 3.10.
  self.http_tokens = asyncio.Queue(maxsize=self.max_conn, loop=self.loop)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/asyncio/queues.py:48: DeprecationWarning: The loop argument is deprecated since Python 3.8, and scheduled for removal in Python 3.10.
  self._finished = locks.Event(loop=loop)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/parfive/downloader.py:88: DeprecationWarning: The loop argument is deprecated since Python 3.8, and scheduled for removal in Python 3.10.
  self.ftp_queue = asyncio.Queue(loop=self.loop)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/asyncio/queues.py:48: DeprecationWarning: The loop argument is deprecated since Python 3.8, and scheduled for removal in Python 3.10.
  self._finished = locks.Event(loop=loop)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/parfive/downloader.py:89: DeprecationWarning: The loop argument is deprecated since Python 3.8, and scheduled for removal in Python 3.10.
  self.ftp_tokens = asyncio.Queue(maxsize=self.max_conn, loop=self.loop)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/asyncio/queues.py:48: DeprecationWarning: The loop argument is deprecated since Python 3.8, and scheduled for removal in Python 3.10.
  self._finished = locks.Event(loop=loop)

Files Downloaded:   0%|          | 0/2 [00:00<?, ?file/s]/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/aiohttp/connector.py:964: DeprecationWarning: The loop argument is deprecated since Python 3.8, and scheduled for removal in Python 3.10.
  hosts = await asyncio.shield(self._resolve_host(
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/aiohttp/connector.py:964: DeprecationWarning: The loop argument is deprecated since Python 3.8, and scheduled for removal in Python 3.10.
  hosts = await asyncio.shield(self._resolve_host(
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/aiohttp/locks.py:21: DeprecationWarning: The loop argument is deprecated since Python 3.8, and scheduled for removal in Python 3.10.
  self._event = asyncio.Event(loop=loop)



hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:   0%|          | 0.00/15.0M [00:00<?, ?B/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:   0%|          | 30.9k/15.0M [00:00<01:22, 182kB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:   1%|          | 86.4k/15.0M [00:00<01:08, 216kB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:   0%|          | 0.00/12.3M [00:00<?, ?B/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:   1%|1         | 165k/15.0M [00:00<00:56, 263kB/s] 

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:   0%|          | 30.3k/12.3M [00:00<00:58, 208kB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:   2%|1         | 267k/15.0M [00:00<00:45, 323kB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:   1%|          | 80.4k/12.3M [00:00<00:51, 236kB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:   3%|2         | 402k/15.0M [00:00<00:36, 394kB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:   1%|1         | 146k/12.3M [00:00<00:44, 275kB/s] 


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:   4%|3         | 564k/15.0M [00:00<00:29, 489kB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:   2%|1         | 223k/12.3M [00:00<00:37, 321kB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:   5%|4         | 748k/15.0M [00:01<00:23, 599kB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:   3%|2         | 318k/12.3M [00:00<00:31, 379kB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:   6%|5         | 871k/15.0M [00:01<00:19, 708kB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:   4%|3         | 441k/12.3M [00:00<00:26, 453kB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:   7%|7         | 1.11M/15.0M [00:01<00:16, 852kB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:   5%|4         | 589k/12.3M [00:01<00:21, 543kB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:   9%|9         | 1.38M/15.0M [00:01<00:13, 1.01MB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:   6%|6         | 758k/12.3M [00:01<00:17, 646kB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:  11%|#1        | 1.72M/15.0M [00:01<00:10, 1.22MB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:   8%|7         | 944k/12.3M [00:01<00:15, 758kB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:  14%|#3        | 2.05M/15.0M [00:01<00:09, 1.38MB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:   9%|9         | 1.16M/12.3M [00:01<00:12, 886kB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:  16%|#6        | 2.45M/15.0M [00:01<00:07, 1.62MB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:  11%|#1        | 1.41M/12.3M [00:01<00:10, 1.04MB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:  19%|#9        | 2.89M/15.0M [00:02<00:06, 1.89MB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:  21%|##1       | 3.17M/15.0M [00:02<00:05, 2.09MB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:  14%|#3        | 1.70M/12.3M [00:01<00:08, 1.21MB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:  24%|##4       | 3.61M/15.0M [00:02<00:04, 2.40MB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:  17%|#6        | 2.05M/12.3M [00:01<00:07, 1.42MB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:  26%|##6       | 3.93M/15.0M [00:02<00:04, 2.58MB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:  20%|#9        | 2.44M/12.3M [00:02<00:05, 1.65MB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:  30%|###       | 4.51M/15.0M [00:02<00:03, 2.94MB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:  23%|##3       | 2.89M/12.3M [00:02<00:04, 1.92MB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:  35%|###4      | 5.20M/15.0M [00:02<00:02, 3.43MB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:  27%|##7       | 3.37M/12.3M [00:02<00:04, 2.19MB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:  38%|###7      | 5.61M/15.0M [00:02<00:02, 3.39MB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:  31%|###1      | 3.86M/12.3M [00:02<00:03, 2.44MB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:  42%|####2     | 6.33M/15.0M [00:02<00:02, 3.75MB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:  35%|###5      | 4.36M/12.3M [00:02<00:02, 2.67MB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:  47%|####7     | 7.09M/15.0M [00:03<00:01, 4.06MB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:  51%|#####     | 7.56M/15.0M [00:03<00:01, 4.20MB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:  40%|###9      | 4.91M/12.3M [00:02<00:02, 2.92MB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:  56%|#####6    | 8.43M/15.0M [00:03<00:01, 4.82MB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:  45%|####4     | 5.53M/12.3M [00:02<00:02, 3.30MB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:  60%|######    | 9.01M/15.0M [00:03<00:01, 5.01MB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:  50%|#####     | 6.17M/12.3M [00:03<00:01, 3.85MB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:  67%|######7   | 10.1M/15.0M [00:03<00:00, 5.69MB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:  54%|#####4    | 6.71M/12.3M [00:03<00:01, 4.04MB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:  74%|#######3  | 11.0M/15.0M [00:03<00:00, 6.44MB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:  61%|######    | 7.48M/12.3M [00:03<00:01, 4.71MB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:  80%|#######9  | 11.9M/15.0M [00:03<00:00, 7.08MB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:  66%|######6   | 8.18M/12.3M [00:03<00:00, 4.95MB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:  85%|########5 | 12.7M/15.0M [00:03<00:00, 6.92MB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:  74%|#######3  | 9.10M/12.3M [00:03<00:00, 5.76MB/s]


hmi_m_45s_2010_08_19_00_01_30_tai_magnetogram.fits:  94%|#########4| 14.1M/15.0M [00:03<00:00, 7.85MB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:  80%|#######9  | 9.82M/12.3M [00:03<00:00, 5.65MB/s]


                                                                                                        
Files Downloaded:  50%|#####     | 1/2 [00:05<00:05,  5.20s/file]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits:  89%|########9 | 11.0M/12.3M [00:03<00:00, 6.22MB/s]

aia_lev1_171a_2010_08_19t00_00_00_34z_image_lev1.fits: 100%|#########9| 12.3M/12.3M [00:03<00:00, 6.79MB/s]

                                                                                                           
Files Downloaded: 100%|##########| 2/2 [00:05<00:00,  3.72s/file]
Files Downloaded: 100%|##########| 2/2 [00:05<00:00,  2.73s/file]

We create a map for each image and resample each one just to reduce the computation time.

map_aia, map_hmi = [m.resample((1024, 1024)*u.pix) for m in sunpy.map.Map(sorted(files))]
# Why do we have to do this?
map_hmi.plot_settings['cmap'] = "hmimag"
map_hmi.plot_settings['norm'] = plt.Normalize(-2000, 2000)

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/astropy/wcs/wcs.py:464: FITSFixedWarning: CRDER1 = 'nan '
a floating-point value was expected.
  wcsprm = _wcs.Wcsprm(header=header_bytes, key=key,
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/astropy/wcs/wcs.py:464: FITSFixedWarning: CRDER2 = 'nan '
a floating-point value was expected.
  wcsprm = _wcs.Wcsprm(header=header_bytes, key=key,
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/astropy/wcs/wcs.py:464: FITSFixedWarning: CRDER1 = 'nan '
a floating-point value was expected.
  wcsprm = _wcs.Wcsprm(header=header_bytes, key=key,
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/astropy/wcs/wcs.py:464: FITSFixedWarning: CRDER2 = 'nan '
a floating-point value was expected.
  wcsprm = _wcs.Wcsprm(header=header_bytes, key=key,

Plot both images side by side.

fig = plt.figure()

ax1 = fig.add_subplot(1, 2, 1, projection=map_aia)
map_aia.plot(axes=ax1)

ax2 = fig.add_subplot(1, 2, 2, projection=map_hmi)
map_hmi.plot(axes=ax2)
../../../_images/sphx_glr_reprojection_align_aia_hmi_001.png

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/astropy/wcs/wcs.py:464: FITSFixedWarning: CRDER1 = 'nan '
a floating-point value was expected.
  wcsprm = _wcs.Wcsprm(header=header_bytes, key=key,
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/astropy/wcs/wcs.py:464: FITSFixedWarning: CRDER2 = 'nan '
a floating-point value was expected.
  wcsprm = _wcs.Wcsprm(header=header_bytes, key=key,

<matplotlib.image.AxesImage object at 0x7f38946b5580>

We can now reproject the HMI image to the WCS of the AIA image. We are using the fast reproject_interp, however the slower but most accurate reproject_exact would also work well here. The reproject_exact function only works when reprojecting between two WCSes with the same observer, which makes it well suited to aligning data.

output, footprint = reproject_interp(map_hmi, map_aia.wcs, map_aia.data.shape)

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/astropy/wcs/wcs.py:464: FITSFixedWarning: CRDER1 = 'nan '
a floating-point value was expected.
  wcsprm = _wcs.Wcsprm(header=header_bytes, key=key,
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/astropy/wcs/wcs.py:464: FITSFixedWarning: CRDER2 = 'nan '
a floating-point value was expected.
  wcsprm = _wcs.Wcsprm(header=header_bytes, key=key,

Construct an output map and set some nice plotting defaults.

out_hmi = sunpy.map.Map(output, map_aia.wcs)
out_hmi.plot_settings['cmap'] = "hmimag"
out_hmi.plot_settings['norm'] = plt.Normalize(-1500, 1500)

fig = plt.figure()
ax1 = fig.add_subplot(1, 2, 1, projection=map_aia)
map_aia.plot(axes=ax1)
ax2 = fig.add_subplot(1, 2, 2, projection=out_hmi)
out_hmi.plot(axes=ax2)
../../../_images/sphx_glr_reprojection_align_aia_hmi_002.png

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/astropy/wcs/wcs.py:684: FITSFixedWarning: 'datfix' made the change 'Set DATE-REF to '1858-11-17' from MJD-REF'.
  warnings.warn(

<matplotlib.image.AxesImage object at 0x7f3894934a30>

As both of these images are now on the same pixel grid we can directly plot them over one another, by setting the transparency of the HMI plot.

fig = plt.figure(figsize=(20, 15))
ax1 = fig.add_subplot(1, 1, 1, projection=map_aia)
map_aia.plot(axes=ax1)
out_hmi.plot(axes=ax1, alpha=0.5)

plt.show()
../../../_images/sphx_glr_reprojection_align_aia_hmi_003.png

Total running time of the script: ( 0 minutes 22.521 seconds)

Gallery generated by Sphinx-Gallery