Saving and loading coordinates with asdf

In this example we are going to look at saving and loading collections of coordinates with asdf.

asdf is a modern file format designed to meet the needs of the astronomy community. It has deep integration with Python and SunPy and Astropy as well as implementations in other languages. It can be used to store known Python objects in a portable, well defined file format. It is primarily useful for storing complex Astropy and SunPy objects in a way that can be loaded back into the same form as they were saved.


This example requires Astropy 3.1 and asdf 2.3.0

import numpy as np
import scipy.optimize
import matplotlib.pyplot as plt

import asdf
import astropy.units as u
import astropy.constants as const
from astropy.coordinates import SkyCoord

from sunpy.coordinates import frames
from import AIA_171_IMAGE


/home/docs/checkouts/ DeprecationWarning: Using or importing the ABCs from 'collections' instead of from '' is deprecated, and in 3.8 it will stop working
  from collections import MutableMapping, Sequence  # noqa

To get started let’s use a function to get the coordinates of a semi-circular loop from this blog post by Will Barnes to generate ourselves some coordinates.

def semi_circular_loop(length: u.m, theta0: u.deg=0*u.deg):
    Return a Heliographic Stonyhurst coordinate object with points of a semi circular loop in it.
    r_sun = const.R_sun

    def r_2_func(x):
        return np.arccos(0.5 * x / - np.pi + / 2. / x

    r_2 = scipy.optimize.bisect(r_2_func,
                       / (2 * np.pi),
                       / np.pi) *
    alpha = np.arccos(0.5 * (r_2 / r_sun).decompose())
    phi = np.linspace(-np.pi * u.rad + alpha, np.pi * u.rad - alpha, 2000)

    # Quadratic formula to find r
    a = 1.
    b = -2 * ( * np.cos(
    c =**2 -**2
    r = (-b + np.sqrt(b**2 - 4 * a * c)) / 2 / a
    # Choose only points above the surface
    i_r = np.where(r > r_sun)
    r = r[i_r]
    phi = phi[i_r]
    hcc_frame = frames.Heliocentric(
        observer=SkyCoord(lon=0 * u.deg, lat=theta0, radius=r_sun, frame='heliographic_stonyhurst'))

    return SkyCoord( * np.sin(,
        y=u.Quantity(r.shape[0] * [0 *]), * np.cos(,

Use this function to generate a SkyCoord object.

loop_coords = semi_circular_loop(500*u.Mm, 30*u.deg)
# print the first and last coordinate point
print(loop_coords[[0, -1]])


<SkyCoord (HeliographicStonyhurst: obstime=None): (lon, lat, radius) in (deg, deg, cm)
    [(-14.09138695, 29.2478234, 6.98633187e+10),
     ( 14.09138695, 29.2478234, 6.98633187e+10)]>

This is a regular coordinate object that can be transformed to other frames or overplotted on images. For instance we could overplot it on an AIA image

aiamap =

ax = plt.subplot(projection=aiamap)


We can now save these loop points to an asdf file to use later. The advantage of saving them to asdf is that all the metadata about the coordinates will be preserved, and when we load the asdf, we will get back an identical HeliographicStonyhurst object.

asdf files save a dictionary to a file, so to save the loop coordinates we need to put them into a dictionary. This becomes what asdf calls a tree.

The asdf file can not save the SkyCoord object, but it can save the underlying frame. Therefore we construct a tree with the frame.

tree = {'loop_points': loop_coords.frame}

with asdf.AsdfFile(tree) as asdf_file:

This asdf file is a portable file and can be safely loaded by anyone with Astropy and SunPy installed. We can reload the file like so:

with"loop_coords.asdf") as input_asdf:
    coords = input_asdf['loop_points']
    new_coords = SkyCoord(coords)

# print the first and last coordinate point
print(new_coords[[0, -1]])


<SkyCoord (HeliographicStonyhurst: obstime=None): (lon, lat, radius) in (deg, deg, cm)
    [(-14.09138695, 29.2478234, 6.98633187e+10),
     ( 14.09138695, 29.2478234, 6.98633187e+10)]>

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

Gallery generated by Sphinx-Gallery