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.


This release adds citation information to the documentation.


This release contains the source for the accepted JOSS paper describing pfsspy.


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.


Bug fixes

  • Fixed some messages in errors raised by functions in pfsspy.utils.


New features

Breaking changes

  • The .al property of pfsspy.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 to pfsspy.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 This means the files used for running the examples will be downloaded and stored in your sunpy 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.


  • 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.


  • Fixed a bug in the GONG synoptic map source where a map failed to load once it had already been loaded once.


  • 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 in the docstring.

  • Added an example showing how to load ADAPT ensemble maps into a CompositeMap

  • Sped up field line expansion factor calculations.


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 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 =, 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

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. Instead Input has a new map property, which returns a SunPy map, which can easily be plotted using

  • pfsspy.Output.plot_source_surface. A map of \(B_{r}\) on the source surface can now be obtained using pfsspy.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 using pfsspy.Output.source_surface_pils, which can then be plotted using ax.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.


  • Improved the error thrown when trying to use :class`pfsspy.tracing.FotranTracer` without the streamtracer module installed.

  • Fixed some layout issues in the documentation.


  • Fix a bug where :class`pfsspy.tracing.FotranTracer` would overwrite the magnetic field values in an Output each time it was used.


  • Reduced the default step size for the FortranTracer from 0.1 to 0.01 to give more resolved field lines by default.


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

Changes to field line objects

Changes to Output

  • 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.


  • Fixed a bug in pfsspy.FieldLine.is_open, where some open field lines were incorrectly calculated to be closed.


  • 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 new pfsspy.FieldLine.solar_footpoint and pfsspy.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 as pfsspy.FieldLines.source_surface_feet can be used, and their values are cached greatly speeding up repeated use.


  • The API for doing magnetic field tracing has changed. The new pfsspy.tracing module contains Tracer classes that are used to perform the tracing. Code needs to be changed from:

    fline = output.trace(x0)


    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 from SkyCoord, but the SkyCoord coordinates are now stored in pfsspy.FieldLine.coords attribute.

  • pfsspy.FieldLine.expansion_factor now returns np.nan instead of None if the field line is closed.

  • pfsspy.FieldLine now has a ~pfsspy.FieldLine.footpoints attribute that returns the footpoint(s) of the field line.


  • pfsspy.Input and pfsspy.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 calling pfsspy.Output.trace


  • Output.plot_source_surface now accepts keyword arguments that are given to Matplotlib to control the plotting of the source surface.


  • Added more explanatory comments to the examples

  • Corrected the dipole solution calculation

  • Added pfsspy.coords.sph2cart to transform from spherical to cartesian coordinates.


  • 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.