Examples using Sunpy for querying and plotting¶

Querying data¶

In [10]:

import numpy as np

import matplotlib.pyplot as plt
import matplotlib.animation as animation
import matplotlib.colors as colors

import astropy.units as u

import sunpy.map
from sunpy.net import Fido
from sunpy.net import attrs as a
import sunpy.visualization.colormaps as cm

import numpy as np import matplotlib.pyplot as plt import matplotlib.animation as animation import matplotlib.colors as colors import astropy.units as u import sunpy.map from sunpy.net import Fido from sunpy.net import attrs as a import sunpy.visualization.colormaps as cm

To download the required data, we use sunpy.net.Fido, a downloader client, to query the Virtual Solar Observatory to acquire AIA and HMI data.

The following is the data available:

HMI Bands (A): 171 (gold), 193 (bronze), 304 (red), 211 (purple), 131 (teal), 335 (blue), 094 (green), 1600 (yellow/green), 1700 (pink)

HMI Products: Magnetogram, Intensitygram, Dopplergram

In [11]:

# data parameters
cadence = a.Sample(24*u.hour)  # querying cadence
start_date = '2021-1-1T12:00:00'  # start date of query
end_date = '2021-01-5T12:00:00'  # end date of query

# AIA specific parameters
aia_inst = a.Instrument.aia
# list of all available wavelengths for reference
wavelength_list = [a.Wavelength(94 * u.angstrom), a.Wavelength(131 * u.angstrom), a.Wavelength(171 * u.angstrom),
                   a.Wavelength(193 * u.angstrom), a.Wavelength(211 * u.angstrom), a.Wavelength(304 * u.angstrom),
                   a.Wavelength(335 * u.angstrom), a.Wavelength(1600 * u.angstrom), a.Wavelength(1700 * u.angstrom)]

# HMI specific parameters
hmi_inst = a.Instrument.hmi
# list of physical observables available
physobs_list = [a.Physobs.los_velocity, a.Physobs.los_magnetic_field, a.Physobs.intensity]

# data parameters cadence = a.Sample(24*u.hour) # querying cadence start_date = '2021-1-1T12:00:00' # start date of query end_date = '2021-01-5T12:00:00' # end date of query # AIA specific parameters aia_inst = a.Instrument.aia # list of all available wavelengths for reference wavelength_list = [a.Wavelength(94 * u.angstrom), a.Wavelength(131 * u.angstrom), a.Wavelength(171 * u.angstrom), a.Wavelength(193 * u.angstrom), a.Wavelength(211 * u.angstrom), a.Wavelength(304 * u.angstrom), a.Wavelength(335 * u.angstrom), a.Wavelength(1600 * u.angstrom), a.Wavelength(1700 * u.angstrom)] # HMI specific parameters hmi_inst = a.Instrument.hmi # list of physical observables available physobs_list = [a.Physobs.los_velocity, a.Physobs.los_magnetic_field, a.Physobs.intensity]

Now we can look at querying for specific information. We can look for results from multiple instruments, wavelengths, and physical observables at once.

It is useful to search for results before download to ensure you don't download excessive amounts of data.

Let's say we only want these specific results from AIA. We search for images at a 24 hour cadence between our start and end times at specified wavelengths.

In [12]:

aia_result = Fido.search(a.Time(start_date, end_date),
                     aia_inst, wavelength_list[0]  | wavelength_list[2] | wavelength_list[3]
                     | wavelength_list[6], cadence)

aia_result = Fido.search(a.Time(start_date, end_date), aia_inst, wavelength_list[0] | wavelength_list[2] | wavelength_list[3] | wavelength_list[6], cadence)

We can do the same with HMI data products.

In [13]:

hmi_result = Fido.search(a.Time(start_date, end_date),
                     a.Instrument.hmi, physobs_list[0] | physobs_list[1] | physobs_list[2], cadence)

hmi_result = Fido.search(a.Time(start_date, end_date), a.Instrument.hmi, physobs_list[0] | physobs_list[1] | physobs_list[2], cadence)

You can also query for both instruments simultaneously.

In [14]:

result = Fido.search(a.Time(start_date, end_date),
                     aia_inst | hmi_inst, wavelength_list[0] | wavelength_list[1] | wavelength_list[2] | wavelength_list[3]
                     | wavelength_list[4] | wavelength_list[5] | wavelength_list[6] | wavelength_list[7] | wavelength_list[8]
                     | physobs_list[0] | physobs_list[1] | physobs_list[2],
                     cadence)

result = Fido.search(a.Time(start_date, end_date), aia_inst | hmi_inst, wavelength_list[0] | wavelength_list[1] | wavelength_list[2] | wavelength_list[3] | wavelength_list[4] | wavelength_list[5] | wavelength_list[6] | wavelength_list[7] | wavelength_list[8] | physobs_list[0] | physobs_list[1] | physobs_list[2], cadence)

/Users/tervin/opt/anaconda3/envs/neid_data/lib/python3.7/site-packages/astropy/table/table.py:3197: FutureWarning: elementwise == comparison failed and returning scalar instead; this will raise an error or perform elementwise comparison in the future.
  result = self.as_array() == other

Next we download the actual data found from our search.

In [15]:

file_download = Fido.fetch(result)

file_download = Fido.fetch(result)

Files Downloaded:   0%|          | 0/54 [00:00<?, ?file/s]

After downloading our data, we want to sort it by date. We then build a map sequence of our data.

In [16]:

map_seq = sunpy.map.Map(sorted(file_download))

map_seq = sunpy.map.Map(sorted(file_download))

We then split our results by data product.

In [17]:

aia_094, aia_131, aia_171, aia_193, aia_211, aia_304, aia_335, aia_1600, aia_1700, hmi_vel, hmi_mag, hmi_int \
    = [], [], [], [], [], [], [], [], [], [], [], []

for i, map_obj in enumerate(map_seq):
    if map_obj.meta['wavelnth'] == 94:
        aia_094.append(map_obj)
    elif map_obj.meta['wavelnth'] == 131:
        aia_131.append(map_obj)
    elif map_obj.meta['wavelnth'] == 171:
        aia_171.append(map_obj)
    elif map_obj.meta['wavelnth'] == 193:
        aia_193.append(map_obj)
    elif map_obj.meta['wavelnth'] == 211:
        aia_211.append(map_obj)
    elif map_obj.meta['wavelnth'] == 304:
        aia_304.append(map_obj)
    elif map_obj.meta['wavelnth'] == 335:
        aia_335.append(map_obj)
    elif map_obj.meta['wavelnth'] == 1600:
        aia_1600.append(map_obj)
    elif map_obj.meta['wavelnth'] == 1700:
        aia_1700.append(map_obj)
    elif map_obj.meta['wavelnth'] == 6173 and map_obj.meta['content'] == 'DOPPLERGRAM':
        hmi_vel.append(map_obj)
    elif map_obj.meta['wavelnth'] == 6173 and map_obj.meta['content'] == 'MAGNETOGRAM':
        hmi_mag.append(map_obj)
    elif map_obj.meta['wavelnth'] == 6173 and map_obj.meta['content'] == 'CONTINUUM INTENSITY':
        hmi_int.append(map_obj)

aia_094, aia_131, aia_171, aia_193, aia_211, aia_304, aia_335, aia_1600, aia_1700, hmi_vel, hmi_mag, hmi_int \ = [], [], [], [], [], [], [], [], [], [], [], [] for i, map_obj in enumerate(map_seq): if map_obj.meta['wavelnth'] == 94: aia_094.append(map_obj) elif map_obj.meta['wavelnth'] == 131: aia_131.append(map_obj) elif map_obj.meta['wavelnth'] == 171: aia_171.append(map_obj) elif map_obj.meta['wavelnth'] == 193: aia_193.append(map_obj) elif map_obj.meta['wavelnth'] == 211: aia_211.append(map_obj) elif map_obj.meta['wavelnth'] == 304: aia_304.append(map_obj) elif map_obj.meta['wavelnth'] == 335: aia_335.append(map_obj) elif map_obj.meta['wavelnth'] == 1600: aia_1600.append(map_obj) elif map_obj.meta['wavelnth'] == 1700: aia_1700.append(map_obj) elif map_obj.meta['wavelnth'] == 6173 and map_obj.meta['content'] == 'DOPPLERGRAM': hmi_vel.append(map_obj) elif map_obj.meta['wavelnth'] == 6173 and map_obj.meta['content'] == 'MAGNETOGRAM': hmi_mag.append(map_obj) elif map_obj.meta['wavelnth'] == 6173 and map_obj.meta['content'] == 'CONTINUUM INTENSITY': hmi_int.append(map_obj)

Plotting images¶

You can plot individual images using inherent Sunpy features.

In [19]:

aia_211[0].peek()  # this plots the aia image with a colorbar, axes labels, and title

aia_211[0].peek() # this plots the aia image with a colorbar, axes labels, and title

In [20]:

aia_211[0].plot()  # this just plots the data, allows you to customize labels, colormap, and title
plt.show()

aia_211[0].plot() # this just plots the data, allows you to customize labels, colormap, and title plt.show()

Plotting HMI images requires us to rotate the image to match the orientation of AIA images.

The same peek and plot methods work for HMI images.

In [21]:

hmi_vel_rot = hmi_vel[0].rotate(order=3)
hmi_vel_rot.peek()

hmi_vel_rot = hmi_vel[0].rotate(order=3) hmi_vel_rot.peek()

HMI images look nicer when customized a bit.

In [22]:

# plotting intensity map

int_map = hmi_int[0].rotate(order=3)
int_map.plot_settings['cmap'] = plt.get_cmap('hinodesotintensity')
int_map.plot()
plt.show()

# plotting intensity map int_map = hmi_int[0].rotate(order=3) int_map.plot_settings['cmap'] = plt.get_cmap('hinodesotintensity') int_map.plot() plt.show()

In [23]:

# plotting magnetic map

mag_map = hmi_mag[0].rotate(order=3)
mag_map.plot_settings['norm'] = plt.Normalize(-100, 100)
mag_map.plot()
plt.show()

# plotting magnetic map mag_map = hmi_mag[0].rotate(order=3) mag_map.plot_settings['norm'] = plt.Normalize(-100, 100) mag_map.plot() plt.show()

In [24]:

# plotting doppler map

vel_map = hmi_vel[0].rotate(order=3)
vel_map.plot_settings['norm'] = plt.Normalize(-5500, 5500)
fig = plt.figure()
vel_map.plot()
plt.show()

# plotting doppler map vel_map = hmi_vel[0].rotate(order=3) vel_map.plot_settings['norm'] = plt.Normalize(-5500, 5500) fig = plt.figure() vel_map.plot() plt.show()

We can now make plots of multiple data types. Don't forget to rotate HMI images such that they are oriented the same as AIA.

In [25]:

fig, axs = plt.subplots(2, 2)
fig.suptitle("Images %s" % (str(aia_171[0].meta['t_obs'])))
axs[0, 0].set_title("%s %s" % (str(aia_171[0].meta['telescop']), str(aia_171[0].meta['wavelnth'])))
axs[0, 1].set_title("%s %s" % (str(aia_211[0].meta['telescop']), str(aia_211[0].meta['wavelnth'])))
axs[1, 0].set_title("%s %s" % (str(hmi_vel[0].meta['telescop']), str(hmi_vel[0].meta['content'])))
axs[1, 1].set_title("%s %s" % (str(hmi_mag[0].meta['telescop']), str(hmi_mag[0].meta['content'])))

for j in range(0, 2):
    for k in range(0, 2):
        axs[j, k].set_xticks([])
        axs[j, k].set_yticks([])

aia_171[0].plot(axes=axs[0, 0], annotate=False)
aia_211[0].plot(axes=axs[0, 1], annotate=False)

# make sure you rotate HMI images
hmi_vel_rot = hmi_vel[0].rotate(order=3)
hmi_vel_rot.plot(axes=axs[1, 0], annotate=False)
hmi_mag_rot = hmi_mag[0].rotate(order=3)
hmi_mag_rot.plot(axes=axs[1, 1], annotate=False)

plt.show()

fig, axs = plt.subplots(2, 2) fig.suptitle("Images %s" % (str(aia_171[0].meta['t_obs']))) axs[0, 0].set_title("%s %s" % (str(aia_171[0].meta['telescop']), str(aia_171[0].meta['wavelnth']))) axs[0, 1].set_title("%s %s" % (str(aia_211[0].meta['telescop']), str(aia_211[0].meta['wavelnth']))) axs[1, 0].set_title("%s %s" % (str(hmi_vel[0].meta['telescop']), str(hmi_vel[0].meta['content']))) axs[1, 1].set_title("%s %s" % (str(hmi_mag[0].meta['telescop']), str(hmi_mag[0].meta['content']))) for j in range(0, 2): for k in range(0, 2): axs[j, k].set_xticks([]) axs[j, k].set_yticks([]) aia_171[0].plot(axes=axs[0, 0], annotate=False) aia_211[0].plot(axes=axs[0, 1], annotate=False) # make sure you rotate HMI images hmi_vel_rot = hmi_vel[0].rotate(order=3) hmi_vel_rot.plot(axes=axs[1, 0], annotate=False) hmi_mag_rot = hmi_mag[0].rotate(order=3) hmi_mag_rot.plot(axes=axs[1, 1], annotate=False) plt.show()

WARNING: SunpyUserWarning: The axes of this map are not aligned to the pixel grid. Plot axes may be incorrect. [sunpy.map.mapbase]
WARNING: SunpyDeprecationWarning: WCSAxes not being used as the axes object for this plot. Support for this is deprecated, and will be removed in sunpy 3.1. To fix this pass set the `projection` keyword to this map when creating the axes. [sunpy.map.mapbase]
WARNING: SunpyUserWarning: The axes of this map are not aligned to the pixel grid. Plot axes may be incorrect. [sunpy.map.mapbase]
WARNING: SunpyDeprecationWarning: WCSAxes not being used as the axes object for this plot. Support for this is deprecated, and will be removed in sunpy 3.1. To fix this pass set the `projection` keyword to this map when creating the axes. [sunpy.map.mapbase]
WARNING: SunpyUserWarning: The axes of this map are not aligned to the pixel grid. Plot axes may be incorrect. [sunpy.map.mapbase]
WARNING: SunpyDeprecationWarning: WCSAxes not being used as the axes object for this plot. Support for this is deprecated, and will be removed in sunpy 3.1. To fix this pass set the `projection` keyword to this map when creating the axes. [sunpy.map.mapbase]
WARNING: SunpyUserWarning: The axes of this map are not aligned to the pixel grid. Plot axes may be incorrect. [sunpy.map.mapbase]
WARNING: SunpyDeprecationWarning: WCSAxes not being used as the axes object for this plot. Support for this is deprecated, and will be removed in sunpy 3.1. To fix this pass set the `projection` keyword to this map when creating the axes. [sunpy.map.mapbase]

In [ ]: