HMI PFSS solutions

Calculating a PFSS solution from a HMI synoptic map.

This example shows how to calculate 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.

import os

import matplotlib.pyplot as plt

import astropy.units as u

import sunpy.map
from sunpy.net import Fido
from sunpy.net import attrs as a

from sunkit_magex import pfss

Set up the search.

The synoptic maps are labelled by Carrington rotation number instead of time.

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(series, crot, a.jsoc.Notify(os.environ["JSOC_EMAIL"]))
files = Fido.fetch(result)

2026-05-20 08:52:44 - drms - INFO: Export request pending. [id=JSOC_20260519_015245_X_IN, status=2]
2026-05-20 08:52:44 - drms - INFO: Waiting for 5 seconds...
INFO: 1 URLs found for download. Full request totaling 4MB [sunpy.net.jsoc.jsoc]
2026-05-20 08:52:50 - sunpy - INFO: 1 URLs found for download. Full request totaling 4MB

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

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

hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits:   0%|          | 1.02k/4.26M [00:00<28:31, 2.49kB/s]

hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits:   2%|▏         | 82.6k/4.26M [00:00<00:21, 197kB/s]

hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits:  10%|▉         | 410k/4.26M [00:00<00:03, 967kB/s]

hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits:  36%|███▌      | 1.54M/4.26M [00:00<00:00, 3.48MB/s]

hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits:  84%|████████▍ | 3.57M/4.26M [00:00<00:00, 7.90MB/s]


Files Downloaded: 100%|██████████| 1/1 [00:01<00:00,  1.51s/file]
Files Downloaded: 100%|██████████| 1/1 [00:01<00:00,  1.51s/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])

Since this map is far to big to calculate a PFSS solution quickly, lets resample it down to a smaller size.

print('Old shape: ', hmi_map.data.shape)
hmi_map = hmi_map.resample([360, 180] * u.pix)
print('New shape: ', hmi_map.data.shape)

Old shape:  (1440, 3600)
New shape:  (180, 360)

Now calculate the PFSS solution

nrho = 35
rss = 2.5
pfss_in = pfss.Input(hmi_map, nrho, rss)
pfss_out = pfss.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

fig = plt.figure()
ax = fig.add_subplot(projection=ss_br)

ss_br.plot(axes=ax)
# Plot the polarity inversion line
ax.plot_coord(pfss_out.source_surface_pils[0])
plt.colorbar()
ax.set_title('Source surface magnetic field')

plt.show()

INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
2026-05-20 08:53:05 - sunpy - INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere.
/home/docs/checkouts/readthedocs.org/user_builds/sunkit-magex/conda/stable/lib/python3.12/site-packages/sunpy/map/mapbase.py:671: SunpyMetadataWarning: Missing metadata for observer: assuming Earth-based observer.
For frame 'heliographic_stonyhurst' the following metadata is missing: dsun_obs,hglt_obs,hgln_obs
For frame 'heliographic_carrington' the following metadata is missing: dsun_obs,crln_obs,crlt_obs

  obs_coord = self.observer_coordinate