1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
|
"""
========================
Reading ADAPT FITS Files
========================
This example demonstrates how to load data from the
Air Force Data Assimilative Photospheric Flux Transport (ADAPT) model into a list of `sunpy.map.Map` objects.
"""
import matplotlib.pyplot as plt
import astropy.units as u
from astropy.io import fits
import sunpy.map
from sunpy.coordinates.sun import carrington_rotation_time
from sunpy.net import Fido
from sunpy.net import attrs as a
###############################################################################
# First we will download an ADAPT FITS file.
# To do this we will use the `sunpy.net.Fido` search interface to search for
# ADAPT data for Carrington Rotation 2193, and the first longitude type (0)
# which means the data will be in a Carrington coordinate frame.
date_start = carrington_rotation_time(2193)
date_end = date_start + 6*u.h
result = Fido.search(a.Time(date_start, date_end), a.Instrument('adapt'), a.adapt.ADAPTLonType("0"))
print(result)
downloaded_file = Fido.fetch(result)
###############################################################################
# ADAPT FITS files contain 12 realizations of synoptic magnetogram
# output as a result of varying model assumptions. `This is explained in detail in this talk
# <https://www.swpc.noaa.gov/sites/default/files/images/u33/SWW_2012_Talk_04_27_2012_Arge.pdf>`__.
#
# Because the array in the FITS file is 3D, it cannot be passed directly to `sunpy.map.Map`,
# as this will take the first slice only, ignoring the other realizations.
# Instead, we'll open the FITS file with `astropy.io.fits` and manually pull out each
# model realization.
adapt_fits = fits.open(downloaded_file[0])
###############################################################################
# ``adapt_fits`` is a list of ``HDU`` objects. The first of these contains
# the 12 realizations of the data and a header with sufficient information to build
# a list of maps. We unpack this information into a list of
# ``(data, header)`` tuples where ``data`` are the different adapt realizations.
data_header_pairs = [(map_slice, adapt_fits[0].header) for map_slice in adapt_fits[0].data]
###############################################################################
# Next, pass this list of tuples to `sunpy.map.Map` to make a list of maps:
adapt_maps = sunpy.map.Map(data_header_pairs)
###############################################################################
# ``adapt_maps`` is now a list of our individual adapt realizations.
# Here, we generate a static plot accessing a subset of the individual maps in turn.
fig = plt.figure(figsize=(5, 10), layout='constrained')
for i, a_map in enumerate(adapt_maps[::4]):
ax = fig.add_subplot(3, 1, i+1, projection=a_map)
ax.tick_params(axis='x', labelsize=6)
im = a_map.plot(axes=ax, cmap='bwr', vmin=-50, vmax=50, title=f"Realization {1+i*4:02d}")
fig.colorbar(im, label='Magnetic Field Strength [G]', orientation='horizontal')
plt.show()
|
