Let's take as an example a recent observation of the Sun made by C. Buil, le 29 november 2020 with a Sol'Ex 1 :
http://www.astrosurf.com/topic/142367-solex-un-instrument-solaire-%C3%A0-petit-prix/
It can be interesting to compare and illustrate the results of the observation in other wavelengths or types of visualization.
The examples presented here are based on the official documentation of the module : https://sunpy.org/.
To display the initial image, we can use the Pillow library which, as its description indicates, allows us to manipulate and process images.
Pillow : https://pillow.readthedocs.io/en/stable/
The Python Imaging Library adds image processing capabilities to your Python interpreter."
#Install Pillow Lib
!pip install --upgrade Pillow;
from PIL import Image
import requests #for download image
#retrieve image
url = 'http://www.astrosurf.com/uploads/monthly_2020_12/_obs_29112022-41.jpg.3bf652c08d329145898043ed0eb4ccb7.jpg'
img = Image.open(requests.get(url, stream=True).raw)
#show image
img
Depending on the use, it is necessary to install the full version of the library, as below.
#Requirements
#python >= 3.6
#astropy >= 3.2
import warnings
warnings.filterwarnings('ignore') # disable warning installation
#Sunpy installation (all packages, including optional dependencies)
!pip install "sunpy[all]"
#imports
#units
import astropy.units as u
#sunpy
from sunpy.net import Fido # => downloader client
from sunpy.net import attrs as a
For a first use, a data sample package is provided with the library to discover its functioning and to carry out tests. We can import this data sample with "import sunpy.data.sample". Thus, we have for example a data set of the AIA instrument for the wavelength 171 with sunpy.data.sample.AIA_171_IMAGE
import sunpy.data.sample.AIA_171_IMAGE #import data sample
To display this data, Sunpy provides a map object. These maps are matrices of spatial data, which are thus declined for a 2D image, a temporal series of 2D images, or 2D images aligned in time. (https://docs.sunpy.org/en/stable/guide/tour.html)
We will be able to create a new map by passing a file name or a list of file names as parameters. As always, we have to start with the import.
import sunpy.data.sample
import sunpy.map
aia = sunpy.map.Map(sunpy.data.sample.AIA_171_IMAGE)
aia.peek()
The call of the peek() method applicable on the aia object allows here to have a fast and easy graphical display. This saves us from having to prepare and configure a matplotlib graph. The latter is accessible on all Sunpy base objects.
However, what we are interested in is to access the data of a certain date.
In order to check if there are data available at the date we are interested in and to download them, we can use Fido. The latter allows us to search and retrieve data, providing a unified interface and allowing the interrogation of several sources simultaneously (https://iopscience.iop.org/article/10.3847/1538-4357/ab4f7a/pdf).
The fido object takes various parameters, for example in our case :
#Get data from AIA with FIDO
result = Fido.search(a.Time('2020/11/29 12:05', '2020/11/29 12:06'),
a.Instrument.aia,
a.Wavelength(171*u.angstrom))
print(result)
Results from 1 Provider: 5 Results from the VSOClient: Start Time ... ... ----------------------- ... 2020-11-29 12:05:09.000 ... 2020-11-29 12:05:21.000 ... 2020-11-29 12:05:33.000 ... 2020-11-29 12:05:45.000 ... 2020-11-29 12:05:57.000 ...
#download files
downloaded_files = Fido.fetch(result)
print(downloaded_files) # output cells are hide here -> too long output. Let's try it.
As mentioned before, we can display the data thanks to the map object of the library, which takes as parameters the elements we have just downloaded. This time by preparing a graph using matplotlib.
import sunpy.map
import matplotlib.pyplot as plt
# prepare map
map = sunpy.map.Map(downloaded_files[0])
#prepare figure
fig = plt.figure(figsize=(16,9))
map.plot()
#show plot
plt.show()
We can go even further by downloading images from two different instruments to superimpose them, with AIA and the HMI (Helioseismic and Magnetic Imager) instrument which maps the magnetic field and the velocity of the surface of the Sun.
To perform this overlay, we need the Reproject library.
The reproject library, when displaying two FITS images, allows to crop the display with the same celestial projection, based on the WCS data of the header, with the possibility to change the orientation, the coordinate system, etc.
%%capture
!pip install reproject
# Imports
from reproject import reproject_interp
# Configuration de la taille de la figure
plt.rcParams['figure.figsize'] = (16, 8)
%%capture
# Configuration of time sample
time = (a.Sample(24 * u.hour) &
a.Time('2020-11-29', '2020-11-29T10:00:00', '2020-11-29') &
a.vso.Extent(0, 0, 0, 0, "FULLDISK"))
# Configuration for aia instrument data
aia = a.Instrument.aia & a.Wavelength(17 * u.nm, 18 * u.nm)
# Configuration for HMI instrument data
hmi = a.Instrument.hmi & a.Physobs.los_magnetic_field
# launch search
res = Fido.search(time, aia | hmi)
# get data
files = Fido.fetch(res[:, 0])
map_aia, map_hmi = [m.resample((1024, 1024)*u.pix) for m in sunpy.map.Map(sorted(files))]
map_hmi.plot_settings['cmap'] = "hmimag"
map_hmi.plot_settings['norm'] = plt.Normalize(-2000, 2000)
fig = plt.figure()
# Atmospheric Imaging Assembly (AIA) Plot
ax1 = fig.add_subplot(1, 2, 1, projection=map_aia)
map_aia.plot(axes=ax1)
# Helioseismic and Magnetic Imager (HMI) Plot
ax2 = fig.add_subplot(1, 2, 2, projection=map_hmi)
map_hmi.plot(axes=ax2)
<matplotlib.image.AxesImage at 0x140577970>
The function used here for reprojection, reproject_interp, takes as parameter in our case:
output, footprint = reproject_interp(map_hmi, map_aia.wcs, map_aia.data.shape)
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()
# plot 1
ax1 = fig.add_subplot(1, 2, 1, projection=map_aia)
map_aia.plot(axes=ax1)
# plot 2
ax2 = fig.add_subplot(1, 2, 2, projection=out_hmi)
out_hmi.plot(axes=ax2)
<matplotlib.image.AxesImage at 0x140678bb0>
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()
# Show
im
This example is only an introduction, the Sunpy library contains many features and data types. Here is for example the list of available instruments, knowing that each instrument has several accessible configurations !
from sunpy.net import Fido, attrs as a
print(a.Instrument)
sunpy.net.attrs.Instrument Specifies the Instrument name for the search. Attribute Name ... --------------------------- ... aia ... bbi ... bcs ... be_continuum ... be_halpha ... bic_hifi ... bigbear ... caii ... cds ... celias ... cerrotololo ... chp ... chrotel ... climso ... cooke ... costep ... cp ... dpm ... eis ... eit ... elteide ... erne ... eve ... eve ... film ... five_12_channelmagnetograph ... foxsi ... gbm ... gfpi ... goes ... golf ... gong ... gris ... ha2 ... hi_c ... hi_c21 ... hmi ... hxeclipse ... hxt ... imax ... impact ... iris ... isoon ... iss ... ivm ... k_cor ... kpdc ... lasco ... learmonth ... longwave_lobe_06 ... longwave_lobe_07 ... longwave_slit_06 ... longwave_slit_07 ... lyra ... maunaloa ... mdi ... mees ... mergedgong ... meudonspectroheliograph ... mk4 ... noaa_indices ... noaa_predict ... norh ... ovsa ... phoenix ... phoka ... plastic ... ptmc ... rhessi ... rhessi ... secchi ... shortwave_lobe_06 ... shortwave_lobe_07 ... shortwave_slit_06 ... shortwave_slit_07 ... six_0_ftshg ... sj ... solarftsspectrometer ... soon ... sot ... sp1 ... sp2 ... spectroheliograph ... spectromagnetograph ... srs_table ... sumer ... suvi ... suvi ... swan ... swap ... swaves ... sxi_0 ... sxt ... tm_1001 ... tm_1010 ... trace ... udaipur ... uvcs ... vault_1999 ... vault_2002 ... vault_2014 ... virgo ... vsm ... wbs ... wispr ... x123 ... xrs ... xrt ... zimpol ...