pfsspy is a python package for carrying out Potential Field Source Surface modelling, a commonly used magnetic field model of the Sun and other stars.
Note
Until pfsspy 1.0 is released, elements of the API are liable to change between versions. A full changelog that lists breaking changes, and how to adapt your code for them can be found at Changelog.
Note
If you find any bugs or have any suggestions for improvement, please raise an issue here: https://github.com/dstansby/pfsspy/issues
Installing¶
pfsspy can be installed from PyPi using
pip install pfsspy
This will install pfsspy and all of its dependencies. In addition to the core dependencies, there are two optional dependencies (numba, streamtracer) that improve code performance. These can be installed with
pip install pfsspy[performance]
Citing¶
If you use pfsspy in work that results in publication, please cite the Journal of Open Source Software paper at https://doi.org/10.21105/joss.02732. A ready made bibtex entry is
@article{Stansby2020,
doi = {10.21105/joss.02732},
url = {https://doi.org/10.21105/joss.02732},
year = {2020},
publisher = {The Open Journal},
volume = {5},
number = {54},
pages = {2732},
author = {David Stansby and Anthony Yeates and Samuel T. Badman},
title = {pfsspy: A Python package for potential field source surface modelling},
journal = {Journal of Open Source Software}
}
Contents¶
Examples¶
Using pfsspy¶
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_20210131_501_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][A
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 0%| | 100/4.26M [00:00<1:33:28, 759B/s][A
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 0%| | 19.4k/4.26M [00:00<00:48, 87.8kB/s][A
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 1%|1 | 50.9k/4.26M [00:00<00:26, 158kB/s] [A
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 2%|2 | 95.1k/4.26M [00:00<00:18, 231kB/s][A
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 4%|3 | 151k/4.26M [00:00<00:13, 303kB/s] [A
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 5%|5 | 224k/4.26M [00:00<00:10, 392kB/s][A
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 7%|7 | 319k/4.26M [00:00<00:07, 507kB/s][A
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 10%|# | 431k/4.26M [00:01<00:06, 621kB/s][A
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 13%|#3 | 569k/4.26M [00:01<00:04, 762kB/s][A
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 17%|#7 | 732k/4.26M [00:01<00:03, 916kB/s][A
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 21%|##1 | 899k/4.26M [00:01<00:03, 1.04MB/s][A
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 26%|##6 | 1.11M/4.26M [00:01<00:02, 1.22MB/s][A
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 31%|###1 | 1.33M/4.26M [00:01<00:02, 1.38MB/s][A
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 37%|###6 | 1.57M/4.26M [00:01<00:01, 1.52MB/s][A
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 44%|####3 | 1.87M/4.26M [00:01<00:01, 1.76MB/s][A
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 50%|##### | 2.14M/4.26M [00:02<00:01, 1.87MB/s][A
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 57%|#####6 | 2.41M/4.26M [00:02<00:00, 1.88MB/s][A
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 63%|######2 | 2.67M/4.26M [00:02<00:00, 2.05MB/s][A
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 69%|######8 | 2.93M/4.26M [00:02<00:00, 2.15MB/s][A
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 76%|#######6 | 3.25M/4.26M [00:02<00:00, 2.36MB/s][A
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 85%|########4 | 3.62M/4.26M [00:02<00:00, 2.67MB/s][A
hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits: 94%|#########3| 4.00M/4.26M [00:02<00:00, 2.99MB/s][A
Files Downloaded: 100%|##########| 1/1 [00:02<00:00, 2.96s/file]
Files Downloaded: 100%|##########| 1/1 [00:02<00:00, 2.96s/file]
[A
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.6/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.6/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.6/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.6/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()
Out:
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/0.6.6/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.6/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.6/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.6/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.6/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.6/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.6/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.173 seconds)
Note
Click here to download the full example code
Magnetic field along a field line¶
How to get the value of the magnetic field along a field line traced through the PFSS solution.
First, import required modules
import astropy.constants as const
import astropy.units as u
from astropy.coordinates import SkyCoord
import matplotlib.pyplot as plt
import numpy as np
import sunpy.map
import pfsspy
from pfsspy import tracing
from pfsspy.sample_data import get_gong_map
Load a GONG magnetic field map. If ‘gong.fits’ is present in the current directory, just use that, otherwise download a sample GONG map.
gong_fname = get_gong_map()
Out:
Files Downloaded: 0%| | 0/1 [00:00<?, ?file/s]
mrzqs200901t1304c2234_022.fits.gz: 0%| | 0.00/242k [00:00<?, ?B/s][A
mrzqs200901t1304c2234_022.fits.gz: 0%| | 100/242k [00:00<07:29, 536B/s][A
mrzqs200901t1304c2234_022.fits.gz: 83%|########3 | 201k/242k [00:00<00:00, 849kB/s][A
Files Downloaded: 100%|##########| 1/1 [00:00<00:00, 1.71file/s]
Files Downloaded: 100%|##########| 1/1 [00:00<00:00, 1.71file/s]
[A
We can now use SunPy to load the GONG fits file, and extract the magnetic field data.
The mean is subtracted to remove the monopole component
gong_map = sunpy.map.Map(gong_fname)
# Remove the mean
gong_map = sunpy.map.Map(gong_map.data - np.mean(gong_map.data), gong_map.meta)
The PFSS solution is calculated on a regular 3D grid in (phi, s, rho), where rho = ln(r), and r is the standard spherical radial coordinate. We need to define the number of rho grid points, and the source surface radius.
nrho = 35
rss = 2.5
From the boundary condition, number of radial grid points, and source surface, we now construct an Input object that stores this information
input = pfsspy.Input(gong_map, nrho, rss)
output = pfsspy.pfss(input)
Now take a seed point, and trace a magnetic field line through the PFSS solution from this point
tracer = tracing.PythonTracer()
r = 1.2 * const.R_sun
lat = 70 * u.deg
lon = 0 * u.deg
seeds = SkyCoord(lon, lat, r, frame=output.coordinate_frame)
field_lines = tracer.trace(seeds, output)
From this field line we can extract the coordinates, and then use
.b_along_fline
to get the components of the magnetic field along the
field line.
From the plot we can see that the non-radial component of the mangetic field goes to zero at the source surface, as expected.
field_line = field_lines[0]
B = field_line.b_along_fline
r = field_line.coords.radius
fig, ax = plt.subplots()
ax.plot(r.to(const.R_sun), B[:, 0], label=r'$B_{r}$')
ax.plot(r.to(const.R_sun), B[:, 1], label=r'$B_{\theta}$')
ax.plot(r.to(const.R_sun), B[:, 2], label=r'$B_{\phi}$')
ax.legend()
ax.set_xlabel(r'r / $R_{\odot}$')
ax.set_ylabel(f'B ({B.unit})')
plt.show()
Total running time of the script: ( 0 minutes 4.958 seconds)
Note
Click here to download the full example code
Open/closed field map¶
Creating an open/closed field map on the solar surface.
First, import required modules
import astropy.units as u
import astropy.constants as const
from astropy.coordinates import SkyCoord
import matplotlib.pyplot as plt
import matplotlib.colors as mcolor
import numpy as np
import sunpy.map
import pfsspy
from pfsspy import tracing
from pfsspy.sample_data import get_gong_map
Load a GONG magnetic field map. If ‘gong.fits’ is present in the current directory, just use that, otherwise download a sample GONG map.
gong_fname = get_gong_map()
We can now use SunPy to load the GONG fits file, and extract the magnetic field data.
The mean is subtracted to enforce div(B) = 0 on the solar surface: n.b. it is not obvious this is the correct way to do this, so use the following lines at your own risk!
gong_map = sunpy.map.Map(gong_fname)
# Remove the mean
gong_map = sunpy.map.Map(gong_map.data - np.mean(gong_map.data), gong_map.meta)
Set the model parameters
nrho = 60
rss = 2.5
Construct the input, and calculate the output solution
input = pfsspy.Input(gong_map, nrho, rss)
output = pfsspy.pfss(input)
Finally, using the 3D magnetic field solution we can trace some field lines. In this case a grid of 90 x 180 points equally gridded in theta and phi are chosen and traced from the source surface outwards.
First, set up the tracing seeds
r = const.R_sun
# Number of steps in cos(latitude)
nsteps = 90
lon_1d = np.linspace(0, 2 * np.pi, nsteps * 2 + 1)
lat_1d = np.arcsin(np.linspace(-1, 1, nsteps + 1))
lon, lat = np.meshgrid(lon_1d, lat_1d, indexing='ij')
lon, lat = lon*u.rad, lat*u.rad
seeds = SkyCoord(lon.ravel(), lat.ravel(), r, frame=output.coordinate_frame)
Trace the field lines
print('Tracing field lines...')
tracer = tracing.FortranTracer(max_steps=2000)
field_lines = tracer.trace(seeds, output)
print('Finished tracing field lines')
Plot the result. The to plot is the input magnetogram, and the bottom plot shows a contour map of the the footpoint polarities, which are +/- 1 for open field regions and 0 for closed field regions.
fig = plt.figure()
m = input.map
ax = fig.add_subplot(2, 1, 1, projection=m)
m.plot()
ax.set_title('Input GONG magnetogram')
ax = fig.add_subplot(2, 1, 2)
cmap = mcolor.ListedColormap(['tab:red', 'black', 'tab:blue'])
norm = mcolor.BoundaryNorm([-1.5, -0.5, 0.5, 1.5], ncolors=3)
pols = field_lines.polarities.reshape(2 * nsteps + 1, nsteps + 1).T
ax.contourf(np.rad2deg(lon_1d), np.sin(lat_1d), pols, norm=norm, cmap=cmap)
ax.set_ylabel('sin(latitude)')
ax.set_title('Open (blue/red) and closed (black) field')
ax.set_aspect(0.5 * 360 / 2)
plt.show()
Total running time of the script: ( 0 minutes 0.000 seconds)
Note
Click here to download the full example code
Dipole source solution¶
A simple example showing how to use pfsspy to compute the solution to a dipole source field.
First, import required modules
import astropy.constants as const
import astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.time import Time
import matplotlib.pyplot as plt
import matplotlib.patches as mpatch
import numpy as np
import sunpy.map
import pfsspy
import pfsspy.coords as coords
To start with we need to construct an input for the PFSS model. To do this, first set up a regular 2D grid in (phi, s), where s = cos(theta) and (phi, theta) are the standard spherical coordinate system angular coordinates. In this case the resolution is (360 x 180).
nphi = 360
ns = 180
phi = np.linspace(0, 2 * np.pi, nphi)
s = np.linspace(-1, 1, ns)
s, phi = np.meshgrid(s, phi)
Now we can take the grid and calculate the boundary condition magnetic field.
def dipole_Br(r, s):
return 2 * s / r**3
br = dipole_Br(1, s)
The PFSS solution is calculated on a regular 3D grid in (phi, s, rho), where rho = ln(r), and r is the standard spherical radial coordinate. We need to define the number of rho grid points, and the source surface radius.
nrho = 30
rss = 2.5
From the boundary condition, number of radial grid points, and source surface, we now construct an Input object that stores this information
header = pfsspy.utils.carr_cea_wcs_header(Time('2020-1-1'), br.shape)
input_map = sunpy.map.Map((br.T, header))
input = pfsspy.Input(input_map, nrho, rss)
Using the Input object, plot the input field
m = input.map
fig = plt.figure()
ax = plt.subplot(projection=m)
m.plot()
plt.colorbar()
ax.set_title('Input dipole field')
Out:
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/0.6.6/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)
Text(0.5, 1.0, 'Input dipole field')
Now calculate the PFSS solution.
output = pfsspy.pfss(input)
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')
Out:
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/0.6.6/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)
Text(0.5, 1.0, 'Source surface magnetic field')
Finally, using the 3D magnetic field solution we can trace some field lines. In this case 32 points equally spaced in theta are chosen and traced from the source surface outwards.
fig, ax = plt.subplots()
ax.set_aspect('equal')
# Take 32 start points spaced equally in theta
r = 1.01 * const.R_sun
lon = np.pi / 2 * u.rad
lat = np.linspace(-np.pi / 2, np.pi / 2, 33) * u.rad
seeds = SkyCoord(lon, lat, r, frame=output.coordinate_frame)
tracer = pfsspy.tracing.PythonTracer()
field_lines = tracer.trace(seeds, output)
for field_line in field_lines:
coords = field_line.coords
coords.representation_type = 'cartesian'
color = {0: 'black', -1: 'tab:blue', 1: 'tab:red'}.get(field_line.polarity)
ax.plot(coords.y / const.R_sun,
coords.z / const.R_sun, color=color)
# Add inner and outer boundary circles
ax.add_patch(mpatch.Circle((0, 0), 1, color='k', fill=False))
ax.add_patch(mpatch.Circle((0, 0), input.grid.rss, color='k', linestyle='--',
fill=False))
ax.set_title('PFSS solution for a dipole source field')
plt.show()
Total running time of the script: ( 0 minutes 5.353 seconds)
Note
Click here to download the full example code
GONG PFSS extrapolation¶
Calculating PFSS solution for a GONG synoptic magnetic field map.
First, import required modules
import astropy.constants as const
import astropy.units as u
from astropy.coordinates import SkyCoord
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import sunpy.map
import pfsspy
from pfsspy import coords
from pfsspy import tracing
from pfsspy.sample_data import get_gong_map
Load a GONG magnetic field map. If ‘gong.fits’ is present in the current directory, just use that, otherwise download a sample GONG map.
gong_fname = get_gong_map()
We can now use SunPy to load the GONG fits file, and extract the magnetic field data.
The mean is subtracted to enforce div(B) = 0 on the solar surface: n.b. it is not obvious this is the correct way to do this, so use the following lines at your own risk!
gong_map = sunpy.map.Map(gong_fname)
# Remove the mean
gong_map = sunpy.map.Map(gong_map.data - np.mean(gong_map.data), gong_map.meta)
The PFSS solution is calculated on a regular 3D grid in (phi, s, rho), where rho = ln(r), and r is the standard spherical radial coordinate. We need to define the number of rho grid points, and the source surface radius.
nrho = 35
rss = 2.5
From the boundary condition, number of radial grid points, and source surface, we now construct an Input object that stores this information
input = pfsspy.Input(gong_map, nrho, rss)
def set_axes_lims(ax):
ax.set_xlim(0, 360)
ax.set_ylim(0, 180)
Using the Input object, plot the input field
m = input.map
fig = plt.figure()
ax = plt.subplot(projection=m)
m.plot()
plt.colorbar()
ax.set_title('Input field')
set_axes_lims(ax)
Out:
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/0.6.6/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)
Now calculate the PFSS solution
output = pfsspy.pfss(input)
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')
set_axes_lims(ax)
Out:
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/0.6.6/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)
It is also easy to plot the magnetic field at an arbitrary height within the PFSS solution.
# Get the radial magnetic field at a given height
ridx = 15
br = output.bc[0][:, :, ridx]
# Create a sunpy Map object using output WCS
br = sunpy.map.Map(br.T, output.source_surface_br.wcs)
# Get the radial coordinate
r = np.exp(output.grid.rc[ridx])
# Create the figure and axes
fig = plt.figure()
ax = plt.subplot(projection=br)
# Plot the source surface map
br.plot(cmap='RdBu')
# Plot formatting
plt.colorbar()
ax.set_title('$B_{r}$ ' + f'at r={r:.2f}' + '$r_{\\odot}$')
set_axes_lims(ax)
Finally, using the 3D magnetic field solution we can trace some field lines. In this case 64 points equally gridded in theta and phi are chosen and traced from the source surface outwards.
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
tracer = tracing.PythonTracer()
r = 1.2 * const.R_sun
lat = np.linspace(-np.pi / 2, np.pi / 2, 8, endpoint=False)
lon = np.linspace(0, 2 * np.pi, 8, endpoint=False)
lat, lon = np.meshgrid(lat, lon, indexing='ij')
lat, lon = lat.ravel() * u.rad, lon.ravel() * u.rad
seeds = SkyCoord(lon, lat, r, frame=output.coordinate_frame)
field_lines = tracer.trace(seeds, output)
for field_line in field_lines:
color = {0: 'black', -1: 'tab:blue', 1: 'tab:red'}.get(field_line.polarity)
coords = field_line.coords
coords.representation_type = 'cartesian'
ax.plot(coords.x / const.R_sun,
coords.y / const.R_sun,
coords.z / const.R_sun,
color=color, linewidth=1)
ax.set_title('PFSS solution')
plt.show()
# sphinx_gallery_thumbnail_number = 4
Total running time of the script: ( 0 minutes 8.560 seconds)
Note
Click here to download the full example code
Overplotting field lines on AIA maps¶
This example shows how to take a PFSS solution, trace some field lines, and overplot the traced field lines on an AIA 193 map.
First, we import the required modules
from datetime import datetime
import os
import astropy.constants as const
import astropy.units as u
from astropy.coordinates import SkyCoord
import matplotlib.pyplot as plt
import numpy as np
import sunpy.map
import sunpy.io.fits
import pfsspy
import pfsspy.tracing as tracing
from pfsspy.sample_data import get_gong_map
Load a GONG magnetic field map. If ‘gong.fits’ is present in the current directory, just use that, otherwise download a sample GONG map.
gong_fname = get_gong_map()
We can now use SunPy to load the GONG fits file, and extract the magnetic field data.
The mean is subtracted to enforce div(B) = 0 on the solar surface: n.b. it is not obvious this is the correct way to do this, so use the following lines at your own risk!
gong_map = sunpy.map.Map(gong_fname)
# Remove the mean
gong_map = sunpy.map.Map(gong_map.data - np.mean(gong_map.data), gong_map.meta)
Load the corresponding AIA 193 map
if not os.path.exists('aia_map.fits'):
import urllib.request
urllib.request.urlretrieve(
'http://jsoc2.stanford.edu/data/aia/synoptic/2020/09/01/H1300/AIA20200901_1300_0193.fits',
'aia_map.fits')
aia = sunpy.map.Map('aia_map.fits')
dtime = aia.date
The PFSS solution is calculated on a regular 3D grid in (phi, s, rho), where rho = ln(r), and r is the standard spherical radial coordinate. We need to define the number of grid points in rho, and the source surface radius.
nrho = 25
rss = 2.5
From the boundary condition, number of radial grid points, and source
surface, we now construct an Input
object that stores this information
input = pfsspy.Input(gong_map, nrho, rss)
Using the Input
object, plot the input photospheric magnetic field
m = input.map
fig = plt.figure()
ax = plt.subplot(projection=m)
m.plot()
plt.colorbar()
ax.set_title('Input field')
Out:
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/0.6.6/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)
Text(0.5, 1.0, 'Input field')
We can also plot the AIA map to give an idea of the global picture. There is a nice active region in the top right of the AIA plot, that can also be seen in the top left of the photospheric field plot above.
ax = plt.subplot(1, 1, 1, projection=aia)
aia.plot(ax)
Out:
<matplotlib.image.AxesImage object at 0x7f10c85b78d0>
Now we construct a 5 x 5 grid of footpoitns to trace some magnetic field lines from. These coordinates are defined in the native Carrington coordinates of the input magnetogram.
# Create 5 points spaced between sin(lat)={0.35, 0.55}
s = np.linspace(0.35, 0.55, 5)
# Create 5 points spaced between long={60, 100} degrees
phi = np.linspace(60, 100, 5)
print(f's = {s}')
print(f'phi = {phi}')
# Make a 2D grid from these 1D points
s, phi = np.meshgrid(s, phi)
# Now convert the points to a coordinate object
lat = np.arcsin(s) * u.rad
lon = phi * u.deg
seeds = SkyCoord(lon.ravel(), lat.ravel(), 1.01 * const.R_sun,
frame=gong_map.coordinate_frame)
Out:
s = [0.35 0.4 0.45 0.5 0.55]
phi = [ 60. 70. 80. 90. 100.]
Plot the magnetogram and the seed footpoints The footpoints are centered around the active region metnioned above.
m = input.map
fig = plt.figure()
ax = plt.subplot(projection=m)
m.plot()
plt.colorbar()
ax.plot_coord(seeds, color='black', marker='o', linewidth=0, markersize=2)
# Set the axes limits. These limits have to be in pixel values
# ax.set_xlim(0, 180)
# ax.set_ylim(45, 135)
ax.set_title('Field line footpoints')
ax.set_ylim(bottom=0)
Out:
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/0.6.6/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)
(0.0, 179.5)
Compute the PFSS solution from the GONG magnetic field input
output = pfsspy.pfss(input)
Trace field lines from the footpoints defined above.
tracer = tracing.PythonTracer()
flines = tracer.trace(seeds, output)
Plot the input GONG magnetic field map, along with the traced mangetic field lines.
m = input.map
fig = plt.figure()
ax = plt.subplot(projection=m)
m.plot()
plt.colorbar()
for fline in flines:
ax.plot_coord(fline.coords, color='black', linewidth=1)
# Set the axes limits. These limits have to be in pixel values
# ax.set_xlim(0, 180)
# ax.set_ylim(45, 135)
ax.set_title('Photospheric field and traced field lines')
Out:
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/0.6.6/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)
Text(0.5, 1.0, 'Photospheric field and traced field lines')
Plot the AIA map, along with the traced magnetic field lines. Inside the loop the field lines are converted to the AIA observer coordinate frame, and then plotted on top of the map.
fig = plt.figure()
ax = plt.subplot(1, 1, 1, projection=aia)
aia.plot(ax)
for fline in flines:
ax.plot_coord(fline.coords, alpha=0.8, linewidth=1, color='white')
# ax.set_xlim(500, 900)
# ax.set_ylim(400, 800)
plt.show()
# sphinx_gallery_thumbnail_number = 5
Total running time of the script: ( 0 minutes 9.922 seconds)
Finding data¶
Examples showing how to find, download, and load magnetograms.
Note
Click here to download the full example code
HMI data¶
How to search for HMI data.
This example shows how to search for, download, and load HMI data, using the
sunpy.net.Fido
interface. HMI data is available via. the Joint Stanford
Operations Center (JSOC).
The polar filled radial magnetic field synoptic maps are obtained using the ‘hmi.synoptic_mr_polfil_720s’ series keyword. Note that they are large (1440 x 720), so you will probably want to downsample them to a smaller resolution to use them to calculate PFSS solutions.
For more information on the maps, see the synoptic maps page on the JSOC site.
First import the required modules
from sunpy.net import Fido, attrs as a
import sunpy.map
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')
Do the search. This will return all the maps in the ‘hmi_mrsynop_small_720s series.’
result = Fido.search(time, series)
print(result)
Out:
Results from 1 Provider:
143 Results from the JSOCClient:
T_REC TELESCOP INSTRUME WAVELNTH CAR_ROT
--------------- -------- --------- -------- -------
Invalid KeyLink SDO/HMI HMI_SIDE1 6173.0 2097
Invalid KeyLink SDO/HMI HMI_SIDE1 6173.0 2098
Invalid KeyLink SDO/HMI HMI_SIDE1 6173.0 2099
Invalid KeyLink SDO/HMI HMI_SIDE1 6173.0 2100
Invalid KeyLink SDO/HMI HMI_SIDE1 6173.0 2101
Invalid KeyLink SDO/HMI HMI_SIDE1 6173.0 2102
Invalid KeyLink SDO/HMI HMI_SIDE1 6173.0 2103
Invalid KeyLink SDO/HMI HMI_SIDE1 6173.0 2104
Invalid KeyLink SDO/HMI HMI_SIDE1 6173.0 2105
Invalid KeyLink SDO/HMI HMI_SIDE1 6173.0 2106
... ... ... ... ...
Invalid KeyLink SDO/HMI HMI_SIDE1 6173.0 2229
Invalid KeyLink SDO/HMI HMI_SIDE1 6173.0 2230
Invalid KeyLink SDO/HMI HMI_SIDE1 6173.0 2231
Invalid KeyLink SDO/HMI HMI_SIDE1 6173.0 2232
Invalid KeyLink SDO/HMI HMI_SIDE1 6173.0 2233
Invalid KeyLink SDO/HMI HMI_SIDE1 6173.0 2234
Invalid KeyLink SDO/HMI HMI_SIDE1 6173.0 2235
Invalid KeyLink SDO/HMI HMI_SIDE1 6173.0 2236
Invalid KeyLink SDO/HMI HMI_SIDE1 6173.0 2237
Invalid KeyLink SDO/HMI HMI_SIDE1 6173.0 2238
Invalid KeyLink SDO/HMI HMI_SIDE1 6173.0 2239
Length = 143 rows
If we just want to download a specific map, we can specify a Carrington rotation number. In addition, downloading files from JSOC requires a notification email. If you use this code, please replace this email address with your own one, registered here: http://jsoc.stanford.edu/ajax/register_email.html
crot = a.jsoc.PrimeKey('CAR_ROT', 2210)
result = Fido.search(time, series, crot, a.jsoc.Notify("jsoc@cadair.com"))
print(result)
Out:
Results from 1 Provider:
1 Results from the JSOCClient:
T_REC TELESCOP INSTRUME WAVELNTH CAR_ROT
--------------- -------- --------- -------- -------
Invalid KeyLink SDO/HMI HMI_SIDE1 6173.0 2210
Download the files. This downloads files to the default sunpy download directory.
files = Fido.fetch(result)
print(files)
Out:
Export request pending. [id=JSOC_20210131_501_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]
Files Downloaded: 100%|##########| 1/1 [00:00<00:00, 7.13file/s]
Files Downloaded: 100%|##########| 1/1 [00:00<00:00, 7.11file/s]
['/home/docs/sunpy/data/hmi.synoptic_mr_polfil_720s.2210.Mr_polfil.fits']
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)
hmi_map.peek()
Out:
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/0.6.6/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.6/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.6/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 12.490 seconds)
Note
Click here to download the full example code
Parsing ADAPT Ensemble .fits files¶
Parse an ADAPT FITS file into a sunpy.map.MapSequence
.
Necessary imports
import sunpy.map
import sunpy.io
import matplotlib.pyplot as plt
from matplotlib import gridspec
from pfsspy.sample_data import get_adapt_map
Load an example ADAPT fits file, utility stored in adapt_helpers.py
adapt_fname = get_adapt_map()
ADAPT synoptic magnetograms contain 12 realizations of synoptic magnetograms output as a result of varying model assumptions. See [here](https://www.swpc.noaa.gov/sites/default/files/images/u33/SWW_2012_Talk_04_27_2012_Arge.pdf))
Because the fits data is 3D, it cannot be passed directly to sunpy.map.Map
,
because this will take the first slice only and the other realizations are
lost. We want to end up with a sunpy.map.MapSequence
containing all these
realiations as individual maps. These maps can then be individually accessed
and PFSS solutions generated from them.
We first read in the fits file using sunpy.io
:
adapt_fits = sunpy.io.fits.read(adapt_fname)
adapt_fits
is a list of HDPair
objects. The first of these contains
the 12 realizations data and a header with sufficient information to build
the MapSequence
. We unpack this HDPair
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 as the argument to sunpy.map.Map
to create
the map sequence :
adaptMapSequence = sunpy.map.Map(data_header_pairs, sequence=True)
adapt_map_sequence
is now a list of our individual adapt realizations.
Note the .peek()` and ``.plot()
methods of MapSequence
returns instances of
sunpy.visualization.MapSequenceAnimator
and
matplotlib.animation.FuncAnimation1
. Here, we generate a static
plot accessing the individual maps in turn :
fig = plt.figure(figsize=(7, 8))
gs = gridspec.GridSpec(4, 3, figure=fig)
for ii, aMap in enumerate(adaptMapSequence):
ax = fig.add_subplot(gs[ii], projection=aMap)
aMap.plot(axes=ax, cmap='bwr', vmin=-2, vmax=2, title=f"Realization {1+ii:02d}")
plt.tight_layout(pad=5, h_pad=2)
plt.show()
Total running time of the script: ( 0 minutes 0.000 seconds)
pfsspy utilities¶
Note
Click here to download the full example code
Re-projecting from CAR to CEA¶
The pfsspy solver takes a cylindrical-equal-area (CEA) projected magnetic field map as input, which is equally spaced in sin(latitude). Some synoptic field maps are equally spaced in latitude however, which is a plate carée (CAR) projection.
This example shows how to use the pfsspy.utils.car_to_cea
function to
reproject a CAR projection to a CEA projection that pfsspy can take as input.
from pfsspy import sample_data
from pfsspy import utils
import matplotlib.pyplot as plt
Load a sample ADAPT map, which has a CAR projection
adapt_maps = utils.load_adapt(sample_data.get_adapt_map())
adapt_map_car = adapt_maps[0]
Out:
Files Downloaded: 0%| | 0/1 [00:00<?, ?file/s]
adapt40311_03k012_202001010000_i00005600n1.fts.gz: 0%| | 0.00/3.11M [00:00<?, ?B/s][A
adapt40311_03k012_202001010000_i00005600n1.fts.gz: 0%| | 100/3.11M [00:00<1:38:29, 526B/s][A
adapt40311_03k012_202001010000_i00005600n1.fts.gz: 8%|8 | 249k/3.11M [00:00<00:02, 1.06MB/s][A
adapt40311_03k012_202001010000_i00005600n1.fts.gz: 27%|##7 | 841k/3.11M [00:00<00:00, 2.90MB/s][A
adapt40311_03k012_202001010000_i00005600n1.fts.gz: 48%|####7 | 1.48M/3.11M [00:00<00:00, 4.13MB/s][A
adapt40311_03k012_202001010000_i00005600n1.fts.gz: 73%|#######3 | 2.28M/3.11M [00:00<00:00, 5.42MB/s][A
Files Downloaded: 100%|##########| 1/1 [00:00<00:00, 1.15file/s]
Files Downloaded: 100%|##########| 1/1 [00:00<00:00, 1.14file/s]
Re-project into a CEA projection
adapt_map_cea = utils.car_to_cea(adapt_map_car)
Out:
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/0.6.6/lib/python3.7/site-packages/sunpy/map/mapbase.py:871: SunpyUserWarning: Missing metadata for observer: assuming Earth-based observer.
warnings.warn(warning_message, SunpyUserWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/0.6.6/lib/python3.7/site-packages/sunpy/map/mapbase.py:871: SunpyUserWarning: Missing metadata for observer: assuming Earth-based observer.
warnings.warn(warning_message, SunpyUserWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/0.6.6/lib/python3.7/site-packages/sunpy/map/mapbase.py:871: SunpyUserWarning: Missing metadata for observer: assuming Earth-based observer.
warnings.warn(warning_message, SunpyUserWarning)
Plot the original map and the reprojected map
plt.figure()
adapt_map_car.plot()
plt.figure()
adapt_map_cea.plot()
plt.show()
Out:
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/0.6.6/lib/python3.7/site-packages/sunpy/map/mapbase.py:871: SunpyUserWarning: Missing metadata for observer: assuming Earth-based observer.
warnings.warn(warning_message, SunpyUserWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pfsspy/envs/0.6.6/lib/python3.7/site-packages/sunpy/map/mapbase.py:871: SunpyUserWarning: Missing metadata for observer: assuming Earth-based observer.
warnings.warn(warning_message, SunpyUserWarning)
Total running time of the script: ( 0 minutes 2.545 seconds)
pfsspy information¶
Examples showing how the internals of pfsspy work.
Note
Click here to download the full example code
pfsspy magnetic field grid¶
A plot of the grid corners, from which the magnetic field values are taken when tracing magnetic field lines.
Notice how the spacing becomes larger at the poles, and closer to the source surface. This is because the grid is equally spaced in \(\cos \theta\) and \(\log r\).
import numpy as np
import matplotlib.pyplot as plt
from pfsspy.grid import Grid
Define the grid spacings
ns = 15
nphi = 360
nr = 10
rss = 2.5
Create the grid
grid = Grid(ns, nphi, nr, rss)
Get the grid edges, and transform to r and theta coordinates
r_edges = np.exp(grid.rg)
theta_edges = np.arccos(grid.sg)
The corners of the grid are where lines of constant (r, theta) intersect, so meshgrid these together to get all the grid corners.
r_grid_points, theta_grid_points = np.meshgrid(r_edges, theta_edges)
Plot the resulting grid corners
fig = plt.figure()
ax = fig.add_subplot(projection='polar')
ax.scatter(theta_grid_points, r_grid_points)
ax.scatter(theta_grid_points + np.pi, r_grid_points, color='C0')
ax.set_ylim(0, 1.1 * rss)
ax.set_theta_zero_location('N')
ax.set_yticks([1, 1.5, 2, 2.5], minor=False)
ax.set_title('$n_{r}$ = ' f'{nr}, ' r'$n_{\theta}$ = ' f'{ns}')
plt.show()
Total running time of the script: ( 0 minutes 0.227 seconds)
Note
Click here to download the full example code
Tracer performance¶
A quick script to compare the performance of the python and fortran tracers.
import timeit
import astropy.units as u
import astropy.coordinates
import numpy as np
import matplotlib.pyplot as plt
import sunpy.map
import pfsspy
Create a dipole map
ntheta = 180
nphi = 360
nr = 50
rss = 2.5
phi = np.linspace(0, 2 * np.pi, nphi)
theta = np.linspace(-np.pi / 2, np.pi / 2, ntheta)
theta, phi = np.meshgrid(theta, phi)
def dipole_Br(r, theta):
return 2 * np.sin(theta) / r**3
br = dipole_Br(1, theta)
br = sunpy.map.Map(br.T, pfsspy.utils.carr_cea_wcs_header('2010-01-01', br.shape))
pfss_input = pfsspy.Input(br, nr, rss)
pfss_output = pfsspy.pfss(pfss_input)
print('Computed PFSS solution')
Trace some field lines
seed0 = np.atleast_2d(np.array([1, 1, 0]))
tracers = [pfsspy.tracing.PythonTracer(),
pfsspy.tracing.FortranTracer()]
nseeds = 2**np.arange(14)
times = [[], []]
for nseed in nseeds:
print(nseed)
seeds = np.repeat(seed0, nseed, axis=0)
r, lat, lon = pfsspy.coords.cart2sph(seeds[:, 0], seeds[:, 1], seeds[:, 2])
r = r * astropy.constants.R_sun
lat = (lat - np.pi / 2) * u.rad
lon = lon * u.rad
seeds = astropy.coordinates.SkyCoord(lon, lat, r, frame=pfss_output.coordinate_frame)
for i, tracer in enumerate(tracers):
if nseed > 64 and i == 0:
continue
t = timeit.timeit(lambda: tracer.trace(seeds, pfss_output), number=1)
times[i].append(t)
Plot the results
fig, ax = plt.subplots()
ax.scatter(nseeds[1:len(times[0])], times[0][1:], label='python')
ax.scatter(nseeds[1:], times[1][1:], label='fortran')
pydt = (times[0][4] - times[0][3]) / (nseeds[4] - nseeds[3])
ax.plot([1, 1e5], [pydt, 1e5 * pydt])
fort0 = times[1][1]
fordt = (times[1][-1] - times[1][-2]) / (nseeds[-1] - nseeds[-2])
ax.plot(np.logspace(0, 5, 100), fort0 + fordt * np.logspace(0, 5, 100))
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('Number of seeds')
ax.set_ylabel('Seconds')
ax.axvline(180 * 360, color='k', linestyle='--', label='180x360 seed points')
ax.legend()
plt.show()
This shows the results of the above script, run on a 2014 MacBook pro with a 2.6 GHz Dual-Core Intel Core i5:
Total running time of the script: ( 0 minutes 0.000 seconds)
API reference¶
pfsspy Package¶
Functions¶
|
Load a saved output file. |
|
Compute PFSS model. |
Classes¶
|
Input to PFSS modelling. |
|
Output of PFSS modelling. |
Class Inheritance Diagram¶
pfsspy.fieldline Module¶
Classes¶
|
A set of closed field lines. |
|
A single magnetic field line. |
|
A collection of |
|
A set of open field lines. |
Class Inheritance Diagram¶
pfsspy.tracing Module¶
Classes¶
|
Tracer using Fortran code. |
|
Tracer using native python code. |
|
Abstract base class for a streamline tracer. |
Class Inheritance Diagram¶
pfsspy.utils Module¶
Functions¶
|
Reproject a plate-carée map in to a cylindrical-equal-area map. |
|
Create a Carrington WCS header for a Cylindrical Equal Area (CEA) projection. |
|
Fix non-compliant FITS metadata in HMI maps. |
|
Returns |
|
Returns |
|
Returns |
|
Parse adapt .fts file as a |
Improving performance¶
numba¶
pfsspy automatically detects an installation of numba, which compiles some of the numerical code to speed up pfss calculations. To enable this simply install numba and use pfsspy as normal.
Streamline tracing¶
pfsspy has two streamline tracers: a pure python pfsspy.tracing.PythonTracer
and a FORTRAN pfsspy.tracing.FortranTracer
. The FORTRAN version is
significantly faster, using the streamtracer package.
Changelog¶
0.6.6¶
Two bugs have been fixed in pfsspy.utils.carr_cea_wcs_header
:
The reference pixel was previously one pixel too large in both longitude and latitude.
The longitude coordinate was previously erroneously translated by one degree.
Both of these are now fixed.
0.6.5¶
This release improves documentation and handling of HMI maps. In particular:
The HMI map downloading example has been updated to use the polar filled data product, which does not have any data missing at the poles.
pfsspy.utils.fix_hmi_meta()
has been added to fix metadata issues in HMI maps. This modifies the metadata of a HMI map to make it FITS compliant, allowing it to be used with pfsspy.
0.6.4¶
This release adds citation information to the documentation.
0.6.3¶
This release contains the source for the accepted JOSS paper describing pfsspy.
0.6.2¶
This release includes several small fixes in response to a review of pfsspy for the Journal of Open Source Software. Thanks to Matthieu Ancellin and Simon Birrer for their helpful feedback!
A permanent code of conduct file has been added to the repository.
Information on how to contribute to pfsspy has been added to the docs.
The example showing the performance of different magnetic field tracers has been fixed.
The docs are now clearer about optional dependencies that can increase performance.
The GONG example data has been updated due to updated data on the remote GONG server.
0.6.1¶
Bug fixes¶
Fixed some messages in errors raised by functions in
pfsspy.utils
.
0.6.0¶
New features¶
The
pfsspy.utils
module has been added, and contains various tools for loading and working with synoptic maps.pfsspy.Output
has a newbunit
property, which returns theUnit
of the input map.Added
pfsspy.Output.get_bvec()
, to sample the magnetic field solution at arbitrary coordinates.Added the
pfsspy.fieldline.FieldLine.b_along_fline
property, to sample the magnetic field along a traced field line.Added a guide to the numerical methods used by pfsspy.
Breaking changes¶
The
.al
property ofpfsspy.Output
is now private, as it is not intended for user access. If you really want to access it, use._al
(but this is now private API and there is no guarantee it will stay or return the same thing in the future).A
ValueError
is now raised if any of the input data topfsspy.Input
is non-finite or NaN. Previously the PFSS computation would run fine, but the output would consist entirely of NaNs.
Behaviour changes¶
The monopole term is now ignored in the PFSS calculation. Previously a non-zero (but small) monopole term would cause floating point precision issues, leading to a very noisy result. Now the monopole term is explicitly removed from the calculation. If your input has a non-zero mean value, pfsspy will issue a warning about this.
The data downloaded by the examples is now automatically downloaded and cached with
sunpy.data.manager
. This means the files used for running the examples will be downloaded and stored in yoursunpy
data directory if they are required.The observer coordinate information in GONG maps is now automatically set to the location of Earth at the time in the map header.
Bug fixes¶
The
date-obs
FITS keyword in GONG maps is now correctly populated.
0.5.3¶
Improved descriptions in the AIA overplotting example.
Fixed the ‘date-obs’ keyword in GONG metadata. Previously this just stored the date and not the time; now both the date and time are properly stored.
Drastically sped up the calculation of source surface and solar surface magnetic field footpoints.
0.5.2¶
Fixed a bug in the GONG synoptic map source where a map failed to load once it had already been loaded once.
0.5.1¶
Fixed some calculations in
pfsspy.carr_cea_wcs_header
, and clarified in the docstring that the input shape must be in[nlon, nlat]
order.Added validation to
pfsspy.Input
to check that the inputted map covers the whole solar surface.Removed ghost cells from
pfsspy.Output.bc
. This changes the shape of the returned arrays by one along some axes.Corrected the shape of
pfsspy.Output.bg
in the docstring.Added an example showing how to load ADAPT ensemble maps into a
CompositeMap
Sped up field line expansion factor calculations.
0.5.0¶
Changes to outputted maps¶
This release largely sees a transition to leveraging Sunpy Map objects. As such, the following changes have been made:
pfsspy.Input
now must take a sunpy.map.GenericMap
as an
input boundary condition (as opposed to a numpy array). To convert a numpy array
to a GenericMap
, the helper function
pfsspy.carr_cea_wcs_header
can be used:
map_date = datetime(...)
br = np.array(...)
header = pfsspy.carr_cea_wcs_header(map_date, br.shape)
m = sunpy.map.Map((br, header))
pfss_input = pfsspy.Input(m, ...)
pfsspy.Output.source_surface_br
now returns a GenericMap
instead of an array. To get the data array use source_surface_br.data
.
The new pfsspy.Output.source_surface_pils
returns the coordinates of
the polarity inversion lines on the source surface.
In favour of directly using the plotting functionality built into SunPy, the following plotting functionality has been removed:
pfsspy.Input.plot_input
. InsteadInput
has a newmap
property, which returns a SunPy map, which can easily be plotted usingsunpy.map.GenericMap.plot
.pfsspy.Output.plot_source_surface
. A map of \(B_{r}\) on the source surface can now be obtained usingpfsspy.Output.source_surface_br
, which again returns a SunPy map.pfsspy.Output.plot_pil
. The coordinates of the polarity inversion lines on the source surface can now be obtained usingpfsspy.Output.source_surface_pils
, which can then be plotted usingax.plot_coord(pil[0])
etc. See the examples section for an example.
Specifying tracing seeds¶
In order to make specifying seeds easier, they must now be a
SkyCoord
object. The coordinates are internally
transformed to the Carrington frame of the PFSS solution, and then traced.
This should make specifying coordinates easier, as lon/lat/r coordinates can be created using:
seeds = astropy.coordinates.SkyCoord(lon, lat, r, frame=output.coordinate_frame)
To convert from the old x, y, z array used for seeds, do:
r, lat, lon = pfsspy.coords.cart2sph
r = r * astropy.constants.R_sun
lat = (lat - np.pi / 2) * u.rad
lon = lon * u.rad
seeds = astropy.coordinates.SkyCoord(lon, lat, r, frame=output.coordinate_frame)
Note that the latitude must be in the range \([-\pi/2, \pi/2]\).
GONG and ADAPT map sources¶
pfsspy now comes with built in sunpy
map sources for GONG and ADAPT synoptic
maps, which automatically fix some non-compliant FITS header values. To use
these, just import pfsspy
and load the .FITS files as normal with sunpy.
Tracing seeds¶
pfsspy.tracing.Tracer
no longer has a transform_seeds
helper method, which
has been replaced by coords_to_xyz
and
pfsspy.tracing.Tracer.xyz_to_coords
. These new methods convert
between SkyCoord
objects, and Cartesian xyz coordinates
of the internal magnetic field grid.
0.4.3¶
Improved the error thrown when trying to use :class`pfsspy.tracing.FotranTracer` without the
streamtracer
module installed.Fixed some layout issues in the documentation.
0.4.2¶
Fix a bug where :class`pfsspy.tracing.FotranTracer` would overwrite the magnetic field values in an
Output
each time it was used.
0.4.1¶
Reduced the default step size for the
FortranTracer
from 0.1 to 0.01 to give more resolved field lines by default.
0.4.0¶
New fortran field line tracer¶
pfsspy.tracing
contains a new tracer,
FortranTracer
. This requires and uses the
streamtracer package
which does streamline tracing rapidly in python-wrapped
fortran code. For large numbers of field lines this results in an ~50x
speedup compared to the PythonTracer
.
Changing existing code to use the new tracer is as easy as swapping out
tracer = pfsspy.tracer.PythonTracer()
for
tracer = pfsspy.tracer.FortranTracer()
. If you notice any issues with the
new tracer, please report them at https://github.com/dstansby/pfsspy/issues.
Changes to field line objects¶
pfsspy.FieldLines
andpfsspy.FieldLine
have moved topfsspy.fieldline.FieldLines
andpfsspy.fieldline.FieldLine
.FieldLines
no longer hassource_surface_feet
andsolar_feet
properties. Instead these have moved to the newpfsspy.fieldline.OpenFieldLines
class. All the open field lines can be accessed from aFieldLines
instance using the newopen_field_lines
property.
Changes to Output
¶
pfsspy.Output.bg
is now returned as a 4D array instead of three 3D arrays. The final index now indexes the vector components; see the docstring for more information.
0.3.2¶
Fixed a bug in
pfsspy.FieldLine.is_open
, where some open field lines were incorrectly calculated to be closed.
0.3.1¶
Fixed a bug that incorrectly set closed line field polarities to -1 or 1 (instead of the correct value of zero).
FieldLine.footpoints
has been removed in favour of the newpfsspy.FieldLine.solar_footpoint
andpfsspy.FieldLine.source_surface_footpoint
. These each return a single footpoint. For a closed field line, see the API docs for further details on this.pfsspy.FieldLines
has been added, as a convenience class to store a collection of field lines. This means convenience attributes such aspfsspy.FieldLines.source_surface_feet
can be used, and their values are cached greatly speeding up repeated use.
0.3.0¶
The API for doing magnetic field tracing has changed. The new
pfsspy.tracing
module containsTracer
classes that are used to perform the tracing. Code needs to be changed from:fline = output.trace(x0)
to:
tracer = pfsspy.tracing.PythonTracer() tracer.trace(x0, output) flines = tracer.xs
Additionally
x0
can be a 2D array that contains multiple seed points to trace, taking advantage of the parallelism of some solvers.The
pfsspy.FieldLine
class no longer inherits fromSkyCoord
, but theSkyCoord
coordinates are now stored inpfsspy.FieldLine.coords
attribute.pfsspy.FieldLine.expansion_factor
now returnsnp.nan
instead ofNone
if the field line is closed.pfsspy.FieldLine
now has a~pfsspy.FieldLine.footpoints
attribute that returns the footpoint(s) of the field line.
0.2.0¶
pfsspy.Input
andpfsspy.Output
now take the optional keyword argument dtime, which stores the datetime on which the magnetic field measurements were made. This is then propagated to the obstime attribute of computed field lines, allowing them to be transformed in to coordinate systems other than Carrington frames.pfsspy.FieldLine
no longer overrrides the SkyCoord__init__
; this should not matter to users, as FieldLine objects are constructed internally by callingpfsspy.Output.trace
0.1.5¶
Output.plot_source_surface
now accepts keyword arguments that are given to Matplotlib to control the plotting of the source surface.
0.1.4¶
Added more explanatory comments to the examples
Corrected the dipole solution calculation
Added
pfsspy.coords.sph2cart
to transform from spherical to cartesian coordinates.
0.1.3¶
pfsspy.Output.plot_pil
now accepts keyword arguments that are given to Matplotlib to control the style of the contour.pfsspy.FieldLine.expansion_factor
is now cached, and is only calculated once if accessed multiple times.
Contributing to pfsspy¶
pfsspy is a community package, and contributions are welcome from anyone. This includes reporting bugs, requesting new features, along with writing code and improving documentation.
Reporting Issues¶
If you use pfsspy and stumble upon a problem, the best way to report it is by opening an issue on the GitHub issue tracker. This way we can help you work around the problem and hopefully fix the problem!
When reporting an issue, please try to provide a short description of the issue with a small code sample, this way we can attempt to reproduce the error. Also provide any error output generated when you encountered the issue, we can use this information to debug the issue.
Requesting new features¶
If there is functionality that is not currently available in pfsspy you can make a feature request on the issue page. Please add as much information as possible regarding the new feature you would like to see.
Improving documentaion¶
If something doesn’t make sense in the online documentation, please open an issue. If you’re comfortable fixing it, please open a pull request to the pfsspy repository.
Contributing code¶
All code contributions to pfsspy are welcome! If you would like to submit a contribution, please open a pull request at the pfsspy repository. It will then be reviewed, and if it is a useful addition merged into the main code branch.
If you are new to open source development, the instructions in these links might be useful for getting started:
Numerical methods¶
For more information on the numerical methods used in the PFSS calculation see
this document
.
Synoptic map FITS conventions¶
FITS is the most common filetype used for the storing of solar images. On this page the FITS metadata conventions for synoptic maps are collected. All of this information can be found in, and is taken from, “Coordinate systems for solar image data (Thompson, 2005)”.
Keyword |
Output |
CRPIXn |
Reference pixel to subtract along axis n. Counts from 1 to N. Integer values refer to the centre of the pixel. |
CRVALn |
Coordinate value of the reference pixel along axis n. |
CDELTn |
Pixel spacing along axis n. |
CTYPEn |
Coordinate axis label for axis n. |
PVi_m |
Additional parameters needed for some coordinate systems. |
Note that CROTAn is ignored in this short guide.
Cylindrical equal area projection¶
In this projection, the latitude pixels are equally spaced in sin(latitude). The reference pixel has to be on the equator, to facilitate alignment with the solar rotation axis.
CDELT2 is set to 180/pi times the pixel spacing in sin(latitude).
CTYPE1 is either ‘HGLN-CEA’ or ‘CRLN-CEA’.
CTYPE2 is either ‘HGLT-CEA’ or ‘CRLT-CEA’.
PVi_1 is set to 1.
LONPOLE is 0.
The abbreviations are “Heliographic Longitude - Cylindrical Equal Area” etc. If the system is heliographic the observer must also be defined in the metadata.