Note
Click here to download the full example code
HMI PFSS solutions
Calculating a PFSS solution from a HMI synoptic map.
This example shows how to calcualte a PFSS solution from a HMI synoptic map. There are a couple of important things that this example shows:
HMI maps have non-standard metadata, so this needs to be fixed
HMI synoptic maps are very big (1440 x 3600), so need to be downsampled in order to calculate the PFSS solution in a reasonable time.
First import the required modules
import os
import astropy.units as u
import matplotlib.pyplot as plt
import sunpy.map
from sunpy.net import Fido
from sunpy.net import attrs as a
import pfsspy
import pfsspy.utils
Set up the search.
Note that for SunPy versions earlier than 2.0, a time attribute is needed to do the search, even if (in this case) it isn’t used, as the synoptic maps are labelled by Carrington rotation number instead of time
time = a.Time('2010/01/01', '2010/01/01')
series = a.jsoc.Series('hmi.synoptic_mr_polfil_720s')
crot = a.jsoc.PrimeKey('CAR_ROT', 2210)
Do the search.
If you use this code, please replace this email address with your own one, registered here: http://jsoc.stanford.edu/ajax/register_email.html
result = Fido.search(time, series, crot,
a.jsoc.Notify(os.environ["JSOC_EMAIL"]))
files = Fido.fetch(result)
Out:
Export request pending. [id=JSOC_20211018_1007_X_IN, status=2]
Waiting for 0 seconds...
1 URLs found for download. Full request totalling 6MB
Files Downloaded: 0%| | 0/1 [00:00<?, ?file/s]
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 0%| | 0.00/6.17M [00:00<?, ?B/s]
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 0%| | 100/6.17M [00:00<2:02:54, 836B/s]
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 0%| | 21.2k/6.17M [00:00<00:58, 105kB/s]
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 1%|1 | 74.7k/6.17M [00:00<00:22, 265kB/s]
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 4%|3 | 220k/6.17M [00:00<00:09, 644kB/s]
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 8%|8 | 515k/6.17M [00:00<00:04, 1.32MB/s]
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 17%|#6 | 1.04M/6.17M [00:00<00:02, 2.37MB/s]
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 28%|##8 | 1.75M/6.17M [00:00<00:01, 3.76MB/s]
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 45%|####4 | 2.75M/6.17M [00:00<00:00, 5.59MB/s]
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 59%|#####8 | 3.61M/6.17M [00:01<00:00, 6.47MB/s]
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 75%|#######5 | 4.63M/6.17M [00:01<00:00, 7.41MB/s]
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 92%|#########2| 5.68M/6.17M [00:01<00:00, 8.09MB/s]
Files Downloaded: 100%|##########| 1/1 [00:01<00:00, 1.40s/file]
Files Downloaded: 100%|##########| 1/1 [00:01<00:00, 1.40s/file]
Read in a file. This will read in the first file downloaded to a sunpy Map object
hmi_map = sunpy.map.Map(files[0])
print('Data shape: ', hmi_map.data.shape)
Out:
Data shape: (1440, 3600)
Since this map is far to big to calculate a PFSS solution quickly, lets resample it down to a smaller size.
hmi_map = hmi_map.resample([360, 180] * u.pix)
print('New shape: ', hmi_map.data.shape)
Out:
New shape: (180, 360)
Now calculate the PFSS solution
nrho = 35
rss = 2.5
pfss_in = pfsspy.Input(hmi_map, nrho, rss)
pfss_out = pfsspy.pfss(pfss_in)
Using the Output object we can plot the source surface field, and the polarity inversion line.
ss_br = pfss_out.source_surface_br
# Create the figure and axes
fig = plt.figure()
ax = plt.subplot(projection=ss_br)
# Plot the source surface map
ss_br.plot()
# Plot the polarity inversion line
ax.plot_coord(pfss_out.source_surface_pils[0])
# Plot formatting
plt.colorbar()
ax.set_title('Source surface magnetic field')
plt.show()
Out:
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/1.0.0/lib/python3.7/site-packages/pfsspy/output.py:95: UserWarning: Could not parse unit string "Mx/cm^2" as a valid FITS unit.
See https://fits.gsfc.nasa.gov/fits_standard.html for the FITS unit standards.
warnings.warn(f'Could not parse unit string "{unit_str}" as a valid FITS unit.\n'
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/1.0.0/lib/python3.7/site-packages/astropy/wcs/wcs.py:710: FITSFixedWarning: 'datfix' made the change 'Invalid DATE-OBS format '2018.11.09_12:30:52_TAI''.
FITSFixedWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/1.0.0/lib/python3.7/site-packages/astropy/wcs/wcs.py:710: FITSFixedWarning: 'unitfix' made the change 'Changed units:
'degree' -> 'deg'.
FITSFixedWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/1.0.0/lib/python3.7/site-packages/sunpy/util/decorators.py:378: SunpyMetadataWarning: Missing metadata for observer: assuming Earth-based observer.
new_val = prop(instance)
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/1.0.0/lib/python3.7/site-packages/pfsspy/output.py:95: UserWarning: Could not parse unit string "Mx/cm^2" as a valid FITS unit.
See https://fits.gsfc.nasa.gov/fits_standard.html for the FITS unit standards.
warnings.warn(f'Could not parse unit string "{unit_str}" as a valid FITS unit.\n'
Total running time of the script: ( 0 minutes 16.177 seconds)