Note
Go to the end to download the full example code.
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)
Export request pending. [id=JSOC_20240513_000389_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.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<21:57, 3.23kB/s]
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 2%|▏ | 99.0k/4.26M [00:00<00:13, 304kB/s]
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 10%|█ | 443k/4.26M [00:00<00:03, 1.25MB/s]
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 42%|████▏ | 1.80M/4.26M [00:00<00:00, 4.93MB/s]
Files Downloaded: 100%|██████████| 1/1 [00:01<00:00, 1.02s/file]
Files Downloaded: 100%|██████████| 1/1 [00:01<00:00, 1.02s/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)
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)
New shape: (180, 360)
Now calculate the PFSS solution
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()
INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
/home/docs/checkouts/readthedocs.org/user_builds/sunkit-magex/conda/latest/lib/python3.10/site-packages/sunpy/map/mapbase.py:633: SunpyMetadataWarning: Missing metadata for observer: assuming Earth-based observer.
For frame 'heliographic_stonyhurst' the following metadata is missing: hgln_obs,dsun_obs,hglt_obs
For frame 'heliographic_carrington' the following metadata is missing: crln_obs,dsun_obs,crlt_obs
obs_coord = self.observer_coordinate
Total running time of the script: (0 minutes 14.519 seconds)