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 astropy.units as u
import matplotlib.pyplot as plt
from sunpy.net import Fido, attrs as a
import sunpy.map
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("jsoc@cadair.com"))
files = Fido.fetch(result)

Out:

Export request pending. [id=JSOC_20210120_1142_X_IN, status=2]
Waiting for 0 seconds...
1 URLs found for download. Full request totalling 4MB

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

hmi.synoptic_mr_polfil_720s.CAR_ROT.Mr_polfil.fits:   0%|          | 0.00/4.26M [00:00<?, ?B/s]

hmi.synoptic_mr_polfil_720s.CAR_ROT.Mr_polfil.fits:   0%|          | 100/4.26M [00:00<1:34:04, 755B/s]

hmi.synoptic_mr_polfil_720s.CAR_ROT.Mr_polfil.fits:   0%|          | 20.9k/4.26M [00:00<00:44, 94.3kB/s]

hmi.synoptic_mr_polfil_720s.CAR_ROT.Mr_polfil.fits:   1%|1         | 55.1k/4.26M [00:00<00:24, 172kB/s] 

hmi.synoptic_mr_polfil_720s.CAR_ROT.Mr_polfil.fits:   2%|2         | 105k/4.26M [00:00<00:16, 257kB/s] 

hmi.synoptic_mr_polfil_720s.CAR_ROT.Mr_polfil.fits:   4%|3         | 162k/4.26M [00:00<00:12, 325kB/s]

hmi.synoptic_mr_polfil_720s.CAR_ROT.Mr_polfil.fits:   6%|5         | 242k/4.26M [00:00<00:09, 427kB/s]

hmi.synoptic_mr_polfil_720s.CAR_ROT.Mr_polfil.fits:   8%|8         | 358k/4.26M [00:00<00:06, 580kB/s]

hmi.synoptic_mr_polfil_720s.CAR_ROT.Mr_polfil.fits:  12%|#1        | 499k/4.26M [00:01<00:05, 747kB/s]

hmi.synoptic_mr_polfil_720s.CAR_ROT.Mr_polfil.fits:  15%|#4        | 638k/4.26M [00:01<00:04, 847kB/s]

hmi.synoptic_mr_polfil_720s.CAR_ROT.Mr_polfil.fits:  19%|#9        | 826k/4.26M [00:01<00:03, 1.04MB/s]

hmi.synoptic_mr_polfil_720s.CAR_ROT.Mr_polfil.fits:  24%|##4       | 1.03M/4.26M [00:01<00:02, 1.21MB/s]

hmi.synoptic_mr_polfil_720s.CAR_ROT.Mr_polfil.fits:  29%|##9       | 1.25M/4.26M [00:01<00:02, 1.38MB/s]

hmi.synoptic_mr_polfil_720s.CAR_ROT.Mr_polfil.fits:  35%|###5      | 1.50M/4.26M [00:01<00:01, 1.51MB/s]

hmi.synoptic_mr_polfil_720s.CAR_ROT.Mr_polfil.fits:  42%|####1     | 1.77M/4.26M [00:01<00:01, 1.67MB/s]

hmi.synoptic_mr_polfil_720s.CAR_ROT.Mr_polfil.fits:  48%|####7     | 2.04M/4.26M [00:01<00:01, 1.83MB/s]

hmi.synoptic_mr_polfil_720s.CAR_ROT.Mr_polfil.fits:  55%|#####5    | 2.36M/4.26M [00:02<00:00, 1.94MB/s]

hmi.synoptic_mr_polfil_720s.CAR_ROT.Mr_polfil.fits:  63%|######3   | 2.68M/4.26M [00:02<00:00, 2.11MB/s]

hmi.synoptic_mr_polfil_720s.CAR_ROT.Mr_polfil.fits:  72%|#######1  | 3.06M/4.26M [00:02<00:00, 2.33MB/s]

hmi.synoptic_mr_polfil_720s.CAR_ROT.Mr_polfil.fits:  82%|########1 | 3.47M/4.26M [00:02<00:00, 2.59MB/s]

hmi.synoptic_mr_polfil_720s.CAR_ROT.Mr_polfil.fits:  89%|########9 | 3.79M/4.26M [00:02<00:00, 2.74MB/s]

hmi.synoptic_mr_polfil_720s.CAR_ROT.Mr_polfil.fits:  96%|#########5| 4.07M/4.26M [00:02<00:00, 2.68MB/s]
Files Downloaded: 100%|##########| 1/1 [00:02<00:00,  2.90s/file]
Files Downloaded: 100%|##########| 1/1 [00:02<00:00,  2.90s/file]


                                                                                                        

Read in a file. This will read in the first file downloaded to a sunpy Map object. Note that HMI maps have several bits of metadata that do not comply to the FITS standard, so we need to fix them first.

hmi_map = sunpy.map.Map(files[0])
pfsspy.utils.fix_hmi_meta(hmi_map)
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:

/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/0.6.5/lib/python3.7/site-packages/astropy/wcs/wcs.py:484: FITSFixedWarning: CRDER1 = 'nan '
a floating-point value was expected.
  colsel=colsel, hdulist=fobj)
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/0.6.5/lib/python3.7/site-packages/astropy/wcs/wcs.py:484: FITSFixedWarning: CRDER2 = 'nan '
a floating-point value was expected.
  colsel=colsel, hdulist=fobj)
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/0.6.5/lib/python3.7/site-packages/astropy/wcs/wcs.py:709: FITSFixedWarning: 'unitfix' made the change 'Changed units:
  'degree' -> 'deg'.
  FITSFixedWarning)
New shape:  (180, 360)

Now calculate the PFSS solution

nrho = 35
rss = 2.5
input = pfsspy.Input(hmi_map, nrho, rss)
output = pfsspy.pfss(input)

Out:

/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/0.6.5/lib/python3.7/site-packages/pfsspy/input.py:40: UserWarning: Input data has a non-zero mean. pfsspy will ignore this non-zero monopole term when calculating the PFSS solution.
  warnings.warn('Input data has a non-zero mean. '

Using the Output object we can plot the source surface field, and the polarity inversion line.

ss_br = output.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(output.source_surface_pils[0])
# Plot formatting
plt.colorbar()
ax.set_title('Source surface magnetic field')

plt.show()
Source surface magnetic field

Out:

/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/0.6.5/lib/python3.7/site-packages/astropy/wcs/wcs.py:484: FITSFixedWarning: CRDER1 = 'nan '
a floating-point value was expected.
  colsel=colsel, hdulist=fobj)
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/0.6.5/lib/python3.7/site-packages/astropy/wcs/wcs.py:484: FITSFixedWarning: CRDER2 = 'nan '
a floating-point value was expected.
  colsel=colsel, hdulist=fobj)
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/0.6.5/lib/python3.7/site-packages/astropy/wcs/wcs.py:709: FITSFixedWarning: 'unitfix' made the change 'Changed units:
  'degree' -> 'deg'.
  FITSFixedWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/0.6.5/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:211: MatplotlibDeprecationWarning: Passing parameters norm and vmin/vmax simultaneously is deprecated since 3.3 and will become an error two minor releases later. Please pass vmin/vmax directly to the norm when creating it.
  return super().imshow(X, *args, origin=origin, **kwargs)
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/0.6.5/lib/python3.7/site-packages/astropy/wcs/wcs.py:484: FITSFixedWarning: CRDER1 = 'nan '
a floating-point value was expected.
  colsel=colsel, hdulist=fobj)
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/0.6.5/lib/python3.7/site-packages/astropy/wcs/wcs.py:484: FITSFixedWarning: CRDER2 = 'nan '
a floating-point value was expected.
  colsel=colsel, hdulist=fobj)
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/0.6.5/lib/python3.7/site-packages/astropy/wcs/wcs.py:709: FITSFixedWarning: 'unitfix' made the change 'Changed units:
  'degree' -> 'deg'.
  FITSFixedWarning)

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

Gallery generated by Sphinx-Gallery