Creating a Full Sun Map with AIA and EUVI

With SDO/AIA and STEREO/A and STEREO/B, it is now possible (given specific dates) to combine combine three EUV images from these satellites to produce a full latitude / longitude map of the Sun.

You will need an active internet connection as well as reproject v0.6 or higher installed.

# sphinx_gallery_thumbnail_number = 4

import matplotlib.pyplot as plt
import numpy as np
from reproject import reproject_interp
from reproject.mosaicking import reproject_and_coadd

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

import sunpy.map
import sunpy.sun
from sunpy.coordinates import get_body_heliographic_stonyhurst
from sunpy.net import Fido
from sunpy.net import attrs as a

To get started, let’s download the data:

stereo = (a.Instrument('EUVI') &
          a.Time('2011-11-01', '2011-11-01T00:10:00'))

aia = (a.Instrument('AIA') &
       a.vso.Sample(24 * u.hour) &
       a.Time('2011-11-01', '2011-11-02'))

wave = a.Wavelength(19.5 * u.nm, 19.5 * u.nm)

res = Fido.search(wave, aia | stereo)

files = Fido.fetch(res)

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/3 [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/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)
/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)


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

aia_lev1_193a_2011_11_01t00_00_07_84z_image_lev1.fits:   0%|          | 48.8k/12.7M [00:00<00:34, 371kB/s]



secchi_l0_b_img_euvi_20111101_20111101_000530_n4eub.fts:   0%|          | 0.00/8.41M [00:00<?, ?B/s]


secchi_l0_a_img_euvi_20111101_20111101_000530_n4eua.fts:   0%|          | 0.00/8.41M [00:00<?, ?B/s]

aia_lev1_193a_2011_11_01t00_00_07_84z_image_lev1.fits:   1%|1         | 133k/12.7M [00:00<00:29, 424kB/s] 


secchi_l0_a_img_euvi_20111101_20111101_000530_n4eua.fts:   2%|1         | 131k/8.41M [00:00<00:08, 956kB/s]

aia_lev1_193a_2011_11_01t00_00_07_84z_image_lev1.fits:   2%|1         | 245k/12.7M [00:00<00:25, 498kB/s]


secchi_l0_a_img_euvi_20111101_20111101_000530_n4eua.fts:   8%|7         | 656k/8.41M [00:00<00:06, 1.27MB/s]

aia_lev1_193a_2011_11_01t00_00_07_84z_image_lev1.fits:   3%|3         | 391k/12.7M [00:00<00:20, 596kB/s]


secchi_l0_a_img_euvi_20111101_20111101_000530_n4eua.fts:  16%|#6        | 1.38M/8.41M [00:00<00:04, 1.68MB/s]



secchi_l0_b_img_euvi_20111101_20111101_000530_n4eub.fts:   2%|2         | 180k/8.41M [00:00<00:19, 419kB/s]


secchi_l0_a_img_euvi_20111101_20111101_000530_n4eua.fts:  26%|##5       | 2.15M/8.41M [00:00<00:02, 2.17MB/s]

aia_lev1_193a_2011_11_01t00_00_07_84z_image_lev1.fits:   4%|4         | 571k/12.7M [00:00<00:16, 718kB/s]



secchi_l0_b_img_euvi_20111101_20111101_000530_n4eub.fts:   8%|7         | 656k/8.41M [00:00<00:13, 575kB/s]


secchi_l0_a_img_euvi_20111101_20111101_000530_n4eua.fts:  33%|###3      | 2.79M/8.41M [00:00<00:02, 2.70MB/s]

aia_lev1_193a_2011_11_01t00_00_07_84z_image_lev1.fits:   6%|6         | 782k/12.7M [00:00<00:13, 858kB/s]



secchi_l0_b_img_euvi_20111101_20111101_000530_n4eub.fts:  11%|#         | 886k/8.41M [00:00<00:10, 742kB/s]


secchi_l0_a_img_euvi_20111101_20111101_000530_n4eua.fts:  46%|####5     | 3.84M/8.41M [00:00<00:01, 3.45MB/s]

aia_lev1_193a_2011_11_01t00_00_07_84z_image_lev1.fits:   8%|8         | 1.04M/12.7M [00:00<00:11, 1.03MB/s]



secchi_l0_b_img_euvi_20111101_20111101_000530_n4eub.fts:  17%|#6        | 1.43M/8.41M [00:00<00:06, 999kB/s]


secchi_l0_a_img_euvi_20111101_20111101_000530_n4eua.fts:  61%|######1   | 5.16M/8.41M [00:00<00:00, 4.44MB/s]



secchi_l0_b_img_euvi_20111101_20111101_000530_n4eub.fts:  22%|##1       | 1.82M/8.41M [00:00<00:05, 1.28MB/s]

aia_lev1_193a_2011_11_01t00_00_07_84z_image_lev1.fits:  10%|#         | 1.33M/12.7M [00:01<00:09, 1.23MB/s]


secchi_l0_a_img_euvi_20111101_20111101_000530_n4eua.fts:  71%|#######1  | 6.01M/8.41M [00:00<00:00, 4.51MB/s]



secchi_l0_b_img_euvi_20111101_20111101_000530_n4eub.fts:  31%|###       | 2.59M/8.41M [00:00<00:03, 1.71MB/s]

aia_lev1_193a_2011_11_01t00_00_07_84z_image_lev1.fits:  13%|#3        | 1.68M/12.7M [00:01<00:07, 1.44MB/s]


secchi_l0_a_img_euvi_20111101_20111101_000530_n4eua.fts:  85%|########5 | 7.15M/8.41M [00:01<00:00, 5.51MB/s]



secchi_l0_b_img_euvi_20111101_20111101_000530_n4eub.fts:  42%|####1     | 3.53M/8.41M [00:01<00:02, 2.26MB/s]


                                                                                                             
Files Downloaded:  33%|###3      | 1/3 [00:02<00:04,  2.21s/file]

aia_lev1_193a_2011_11_01t00_00_07_84z_image_lev1.fits:  17%|#6        | 2.11M/12.7M [00:01<00:06, 1.73MB/s]



secchi_l0_b_img_euvi_20111101_20111101_000530_n4eub.fts:  55%|#####5    | 4.64M/8.41M [00:01<00:01, 2.97MB/s]

aia_lev1_193a_2011_11_01t00_00_07_84z_image_lev1.fits:  20%|##        | 2.61M/12.7M [00:01<00:04, 2.07MB/s]



secchi_l0_b_img_euvi_20111101_20111101_000530_n4eub.fts:  64%|######3   | 5.38M/8.41M [00:01<00:00, 3.47MB/s]

aia_lev1_193a_2011_11_01t00_00_07_84z_image_lev1.fits:  25%|##5       | 3.21M/12.7M [00:01<00:03, 2.58MB/s]



secchi_l0_b_img_euvi_20111101_20111101_000530_n4eub.fts:  72%|#######2  | 6.07M/8.41M [00:01<00:00, 4.07MB/s]

aia_lev1_193a_2011_11_01t00_00_07_84z_image_lev1.fits:  29%|##9       | 3.70M/12.7M [00:01<00:03, 2.90MB/s]



secchi_l0_b_img_euvi_20111101_20111101_000530_n4eub.fts:  80%|########  | 6.76M/8.41M [00:01<00:00, 4.43MB/s]

aia_lev1_193a_2011_11_01t00_00_07_84z_image_lev1.fits:  35%|###4      | 4.41M/12.7M [00:01<00:02, 3.52MB/s]



secchi_l0_b_img_euvi_20111101_20111101_000530_n4eub.fts:  92%|#########2| 7.77M/8.41M [00:01<00:00, 5.33MB/s]



                                                                                                             
Files Downloaded:  67%|######6   | 2/3 [00:02<00:01,  1.70s/file]

aia_lev1_193a_2011_11_01t00_00_07_84z_image_lev1.fits:  41%|####      | 5.22M/12.7M [00:01<00:01, 4.11MB/s]

aia_lev1_193a_2011_11_01t00_00_07_84z_image_lev1.fits:  49%|####8     | 6.20M/12.7M [00:02<00:01, 4.98MB/s]

aia_lev1_193a_2011_11_01t00_00_07_84z_image_lev1.fits:  57%|#####7    | 7.29M/12.7M [00:02<00:00, 5.72MB/s]

aia_lev1_193a_2011_11_01t00_00_07_84z_image_lev1.fits:  67%|######6   | 8.51M/12.7M [00:02<00:00, 6.80MB/s]

aia_lev1_193a_2011_11_01t00_00_07_84z_image_lev1.fits:  77%|#######6  | 9.77M/12.7M [00:02<00:00, 7.60MB/s]

aia_lev1_193a_2011_11_01t00_00_07_84z_image_lev1.fits:  88%|########7 | 11.2M/12.7M [00:02<00:00, 8.79MB/s]

aia_lev1_193a_2011_11_01t00_00_07_84z_image_lev1.fits:  99%|#########9| 12.6M/12.7M [00:02<00:00, 10.0MB/s]

                                                                                                           
Files Downloaded: 100%|##########| 3/3 [00:03<00:00,  1.41s/file]
Files Downloaded: 100%|##########| 3/3 [00:03<00:00,  1.15s/file]

Next we create a SunPy map for each of the files.

maps = sunpy.map.Map(sorted(files))

To reduce memory consumption we also downsample these maps before continuing, you can disable this.

maps = [m.resample((1024, 1024)*u.pix) for m in maps]

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/astropy/wcs/wcs.py:464: FITSFixedWarning: CROTA = -2.59161103276
keyword looks very much like CROTAn but isn't.
  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: CROTA = -2.59161103276
keyword looks very much like CROTAn but isn't.
  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: CROTA = 4.34749588558
keyword looks very much like CROTAn but isn't.
  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: CROTA = 4.34749588558
keyword looks very much like CROTAn but isn't.
  wcsprm = _wcs.Wcsprm(header=header_bytes, key=key,

When combining these images all three need to assume the same radius of the Sun for the data. The AIA images specify a slightly different value than the IAU 2015 constant. To avoid coordinate transformation issues we reset this here.

maps[0].meta['rsun_ref'] = sunpy.sun.constants.radius.to_value(u.m)

Next we will plot the locations of the three spacecraft with respect to the Sun so we can easily see the relative separations.

earth = get_body_heliographic_stonyhurst('earth', maps[0].date)

plt.figure(figsize=(8, 8))
r_unit = u.AU

ax = plt.subplot(projection='polar')
circle = plt.Circle((0.0, 0.0), (10*u.Rsun).to_value(r_unit),
                    transform=ax.transProjectionAffine + ax.transAxes, color="yellow",
                    alpha=1, label="Sun")
ax.add_artist(circle)
ax.text(earth.lon.to_value("rad")+0.05, earth.radius.to_value(r_unit), "Earth")

for this_satellite, this_coord in [(m.observatory, m.observer_coordinate) for m in maps]:
    plt.polar(this_coord.lon.to('rad'), this_coord.radius.to(r_unit), 'o', label=this_satellite)

ax.set_theta_zero_location("S")
ax.set_rlim(0, 1.3)

plt.legend()
plt.show()
../../../_images/sphx_glr_reprojection_aia_euvi_mosaic_001.png

The next step is to calculate the output coordinate system for the combined map. We select a heliographic Stonyhurst frame, and a Plate Carree (CAR) projection, and generate a header using sunpy.map.make_fitswcs_header and then construct a World Coordinate System (WCS) object for that header.

shape_out = (180, 360)  # This is set deliberately low to reduce memory consumption

header = sunpy.map.make_fitswcs_header(shape_out,
                                       SkyCoord(0, 0, unit=u.deg,
                                                frame="heliographic_stonyhurst",
                                                obstime=maps[0].date),
                                       scale=[180 / shape_out[0],
                                              360 / shape_out[1]] * u.deg / u.pix,
                                       wavelength=int(maps[0].meta['wavelnth']) * u.AA,
                                       projection_code="CAR")
out_wcs = WCS(header)

Next we call the reproject.mosaicking.reproject_and_coadd function, which takes a list of maps, and the desired output WCS and array shape.

array, footprint = reproject_and_coadd(maps, out_wcs, shape_out,
                                       reproject_function=reproject_interp)

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/astropy/wcs/wcs.py:464: FITSFixedWarning: CROTA = -2.59161103276
keyword looks very much like CROTAn but isn't.
  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: CROTA = 4.34749588558
keyword looks very much like CROTAn but isn't.
  wcsprm = _wcs.Wcsprm(header=header_bytes, key=key,

To display the output we construct a new map using the new array and our generated header. We also borrow the plot settings from the AIA map.

outmap = sunpy.map.Map((array, header))
outmap.plot_settings = maps[0].plot_settings

outmap.plot()
plt.show()
../../../_images/sphx_glr_reprojection_aia_euvi_mosaic_002.png

Improving the Output

As you can see this leaves a little to be desired. To reduce the obvious warping towards the points which are close to the limb in the input images, we can define a set of weights to use when co-adding the output arrays. To reduce this warping we want to calculate an set of weights which highly weigh points close to the centre of the disk in the input image.

We can achieve this by using SunPy’s coordinate framework. First we calculate all the world coordinates for all the pixels in all three input maps.

coordinates = tuple(map(sunpy.map.all_coordinates_from_map, maps))

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/astropy/wcs/wcs.py:464: FITSFixedWarning: CROTA = -2.59161103276
keyword looks very much like CROTAn but isn't.
  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: CROTA = -2.59161103276
keyword looks very much like CROTAn but isn't.
  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: CROTA = 4.34749588558
keyword looks very much like CROTAn but isn't.
  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: CROTA = 4.34749588558
keyword looks very much like CROTAn but isn't.
  wcsprm = _wcs.Wcsprm(header=header_bytes, key=key,

To get a weighting which is high close to disk centre and low towards the limb, we can use the Z coordinate in the heliocentric frame. This coordinate is the distance of the sphere from the centre of the Sun towards the observer.

weights = [coord.transform_to("heliocentric").z.value for coord in coordinates]

These weights are good, but they are better if the ramp down is a little smoother, and more biased to the centre. Also we can scale them to the range 0-1, and set any off disk (NaN) regions to 0.

weights = [(w / np.nanmax(w)) ** 3 for w in weights]
for w in weights:
    w[np.isnan(w)] = 0

plt.figure()
plt.imshow(weights[0])
plt.colorbar()
plt.show()
../../../_images/sphx_glr_reprojection_aia_euvi_mosaic_003.png

Now we can rerun the reprojection. This time we also set match_background=True which scales the images by a single scaling factor so they are of similar brightness. We also set background_reference=0 which uses the AIA map as the reference for the background scaling.

Here we are using the fastest but least accurate method of reprojection, reproject.reproject_interp, a more accurate but slower method is reproject.reproject_adaptive.

array, _ = reproject_and_coadd(maps, out_wcs, shape_out,
                               input_weights=weights,
                               reproject_function=reproject_interp,
                               match_background=True,
                               background_reference=0)

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.8/site-packages/astropy/wcs/wcs.py:464: FITSFixedWarning: CROTA = -2.59161103276
keyword looks very much like CROTAn but isn't.
  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: CROTA = 4.34749588558
keyword looks very much like CROTAn but isn't.
  wcsprm = _wcs.Wcsprm(header=header_bytes, key=key,

Once again we create a new map, and this time we customise the plot a little.

outmap = sunpy.map.Map((array, header))
outmap.plot_settings = maps[0].plot_settings
outmap.nickname = 'AIA + EUVI/A + EUVI/B'

plt.figure(figsize=(10, 5))
ax = plt.subplot(projection=out_wcs)
im = outmap.plot(vmin=400)

lon, lat = ax.coords

lon.set_coord_type("longitude")
lon.coord_wrap = 180
lon.set_format_unit(u.deg)
lat.set_coord_type("latitude")
lat.set_format_unit(u.deg)

lon.set_axislabel('Heliographic Longitude', minpad=0.8)
lat.set_axislabel('Heliographic Latitude', minpad=0.9)
lon.set_ticks(spacing=25*u.deg, color='k')
lat.set_ticks(spacing=15*u.deg, color='k')

plt.colorbar(im, ax=ax)

# Reset the view to pixel centers
_ = ax.axis((0, shape_out[1], 0, shape_out[0]))

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

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

Gallery generated by Sphinx-Gallery