PSRCHIVE from Python — a guided tour¶

Everything the command-line tools (vap, pav, pam, paz, psradd, pat) do is built on the PSRCHIVE C++ library, and that same library is exposed in Python as the psrchive package. Same Archive / Integration / Profile objects, same methods — but now you can mix them with numpy and matplotlib and script whatever you like.

This notebook tours the bindings on the real J1903-7051 data from Part 1 of the tutorial, and finishes on one of the hidden-picture archives so you can do a reveal entirely in Python.

Run this notebook with the PSRCHIVE environment loaded (the same one in which import psrchive works at a terminal). If the first cell fails to import, ask an instructor how to load the environment on this machine.

In [ ]:
import psrchive
import numpy as np
import matplotlib.pyplot as plt

# Paths, relative to the tutorial folder. Adjust if you run from elsewhere.
REAL = "real_data/J1903-7051"
ARCHIVES = "archives"

print("psrchive loaded:", psrchive.__file__)

1. Load an archive¶

psrchive.Archive_load(filename) returns an Archive object for any PSRFITS file (.ar, .std, .sm, …).

In [ ]:
import glob
# Pick one real observation (the first epoch).
obs_file = sorted(glob.glob(f"{REAL}/data/*.ar"))[0]
arch = psrchive.Archive_load(obs_file)
print("loaded:", obs_file.split('/')[-1])
print(type(arch))

2. Read the header¶

Every field vap/psredit shows has a get_…() method. Here are the common ones — this is the Python equivalent of vap -c name,nbin,nchan,nsubint,npol,freq,bw,dm,state file.ar.

In [ ]:
print("source        :", arch.get_source())
print("telescope     :", arch.get_telescope())
print("nbin          :", arch.get_nbin())
print("nchan         :", arch.get_nchan())
print("nsubint       :", arch.get_nsubint())
print("npol          :", arch.get_npol())
print("centre freq   :", arch.get_centre_frequency(), "MHz")
print("bandwidth     :", arch.get_bandwidth(), "MHz")
print("DM            :", arch.get_dispersion_measure(), "pc/cm^3")
print("state         :", arch.get_state())
print("dedispersed?  :", arch.get_dedispersed())
print("length (s)    :", arch.integration_length())

3. The data cube and the weights¶

arch.get_data() returns one numpy array of shape (nsubint, npol, nchan, nbin) — the whole archive at once. arch.get_weights() returns the per-(subint, channel) weights; a zero weight means that channel was zapped (RFI-flagged).

In [ ]:
data = arch.get_data()
print("data shape :", data.shape, " dtype:", data.dtype)

w = arch.get_weights()
print("weights shape    :", w.shape)
print("fraction zapped  :", float((w == 0).mean()))

4. Plot the pulse profile¶

These real archives are already fully scrunched (nsubint = nchan = 1), so the profile is just data[0, 0] (subint 0, Stokes I). This is the Python version of pav -D.

In [ ]:
prof = data[0, 0, 0]                   # Stokes I profile (subint 0, pol 0, chan 0)
phase = np.linspace(0, 1, prof.size, endpoint=False)

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(phase, prof, lw=1)
ax.set_xlabel("pulse phase"); ax.set_ylabel("flux (arb.)")
ax.set_title(f"{arch.get_source()} — total intensity")
plt.show()

5. Polarisation profile (Stokes I, L, V)¶

The real data is in the Stokes state with 4 polarisations, so data[0, :, 0, :] holds [I, Q, U, V] (subint 0, channel 0). Linear polarisation is L = sqrt(Q^2 + U^2). This reproduces what pav -S draws.

In [ ]:
I, Q, U, V = data[0, :, 0, :]           # each is (nbin,)
L = np.hypot(Q, U)

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(phase, I, 'k', lw=1.2, label="I (total)")
ax.plot(phase, L, 'r', lw=1.0, label="L (linear)")
ax.plot(phase, V, 'b', lw=1.0, label="V (circular)")
ax.set_xlabel("pulse phase"); ax.set_ylabel("flux (arb.)")
ax.set_title(f"{arch.get_source()} — polarisation")
ax.legend(); plt.show()

6. Scrunching and dedispersion (in place)¶

The Archive methods mutate the object, mirroring the pam flags:

Python pam
arch.dedisperse() -D
arch.tscrunch() -T
arch.fscrunch() -F
arch.pscrunch() -p
arch.bscrunch(n) combine n bins
arch.convert_state('Stokes') -S
arch.remove_baseline() (baseline subtraction)

Because they change the object, reload a fresh copy before each new experiment. Let's use a 16-channel file so there's a frequency axis to see, and build a frequency-vs-phase image (the Python version of pav -dG).

In [ ]:
f16 = sorted(glob.glob(f"{REAL}/data_16ch/*.ar"))[0]
a16 = psrchive.Archive_load(f16)
a16.remove_baseline()
a16.dedisperse()        # align channels in phase
a16.pscrunch()          # total intensity
a16.tscrunch()          # one subint
img = a16.get_data()[0, 0]   # (nchan, nbin)

fig, ax = plt.subplots(figsize=(8, 4))
ax.imshow(img, aspect="auto", origin="lower", cmap="magma",
          extent=[0, 1, 0, a16.get_nchan()])
ax.set_xlabel("pulse phase"); ax.set_ylabel("channel")
ax.set_title("frequency vs phase (dedispersed)")
plt.show()

7. Per-subint / per-channel access¶

get_data() is convenient but copies the whole cube. For big archives you can walk the structure instead: an Archive holds Integrations (subints), each holding Profiles (one per pol/channel).

In [ ]:
a = psrchive.Archive_load(f16)
subint0 = a.get_Integration(0)
print("subint duration (s):", subint0.get_duration())
print("centre freq of chan 8:", subint0.get_centre_frequency(8), "MHz")

prof_obj = subint0.get_Profile(0, 8)   # (pol=0, chan=8)
amps = prof_obj.get_amps()             # numpy VIEW into the buffer
print("profile shape:", amps.shape, " peak:", float(amps.max()))

get_amps() returns a view, so you can modify the data in place (amps[:] = new_values) — this is exactly how the hidden-picture archives were painted. After editing, arch.unload('new.ar') writes it back out.

A common idiom — the per-channel frequency array:

In [ ]:
freqs = np.array([subint0.get_centre_frequency(c)
                  for c in range(a.get_nchan())])
print("channel freqs (MHz):", np.round(freqs, 1))

8. State conversion + a reveal in Python¶

Now to one of the hidden-picture archives. Some pictures are stored in a particular polarisation, in the instrument's native (coherency) state. Convert to Stokes, then look at the polarisation you want — the same logic as pam -DTS + psrplot -p freq -c "pol=3".

Try it on archives/delta.ar (read its header comment first for the hint). The reveal lives in Stokes V (index 3).

In [ ]:
rev = psrchive.Archive_load(f"{ARCHIVES}/delta.ar")
print("hint:", rev.get_source())

rev.remove_baseline()
rev.dedisperse()
rev.tscrunch()
rev.convert_state("Stokes")          # native -> Stokes
v_img = rev.get_data()[0, 3]         # Stokes V, (nchan, nbin)

fig, ax = plt.subplots(figsize=(8, 6))
ax.imshow(v_img, aspect="auto", origin="lower", cmap="gray")
ax.set_xlabel("pulse phase bin"); ax.set_ylabel("channel")
ax.set_title("Stokes V — what develops?")
plt.show()

Your turn: repeat the trick on archives/alpha.ar — but that one is a total-intensity image, so you only need dedisperse() + tscrunch() and then imshow of get_data()[0, 0]. Compare what you get to running pav -GTd archives/alpha.ar on the command line: the Python and CLI results are identical.

9. Generating TOAs from Python¶

pat is the one major tool whose Python equivalent is awkward to drive, so in practice people shell out to it from a notebook and parse the resulting .tim. Here we reuse the smoothed template you can make on the command line (Part 1, §1.5); if you haven't made one yet, this cell will tell you.

In [ ]:
import subprocess, os

template = f"{REAL}/data/grand.average.Tp.sm"   # made in Part 1 §1.6
if not os.path.exists(template):
    print(f"No template at {template!r} yet — make one in Part 1 §1.6, "
          "or point this at any .sm/.std file.")
else:
    obs = sorted(glob.glob(f"{REAL}/data/*.ar"))
    out = subprocess.run(
        ["pat", "-f", "tempo2", "-s", template, *obs],
        capture_output=True, text=True)
    tim_lines = [l for l in out.stdout.splitlines()
                 if l and not l.startswith("FORMAT")]
    print(f"{len(tim_lines)} TOAs generated. First few:")
    print("\n".join(tim_lines[:3]))

Parse the TOA uncertainties (4th column) and plot them — a quick look at how well each epoch was timed.

In [ ]:
if os.path.exists(template):
    mjd = np.array([float(l.split()[2]) for l in tim_lines])
    err = np.array([float(l.split()[3]) for l in tim_lines])  # microseconds
    fig, ax = plt.subplots(figsize=(8, 4))
    ax.plot(mjd, err, 'o', ms=4)
    ax.set_xlabel("MJD"); ax.set_ylabel("TOA uncertainty (us)")
    ax.set_title("Per-epoch timing precision")
    plt.show()
    print("median TOA uncertainty: %.3f us" % np.median(err))

10. Saving your work¶

Any Archive you've modified can be written back out; PSRCHIVE picks the format from the extension.

arch.unload("my_result.ar")

Where to read more¶

  • The HTML manuals document every class/method: https://psrchive.sourceforge.net/manuals/
  • In a notebook, tab-completion on an Archive object (arch.<TAB>) lists the available methods; help(arch.tscrunch) shows docstrings.

Exercises¶

  1. Open archives/beta.ar. It's a time vs phase image — dedisperse, fscrunch() (not tscrunch!), and imshow get_data()[0, 0] transposed appropriately. What develops?
  2. Open archives/zeta.ar. Without scrunching frequency, dedisperse + tscrunch + pscrunch and imshow the (nchan, nbin) plane. Then separately fscrunch it and plot the Stokes I and V profiles on the same axes. You should find two different things.
  3. Loop over every get_… getter on an Archive and print what it returns:
    for name in dir(arch):
    if name.startswith("get_") and "Integration" not in name:
        try: print(name, "=", getattr(arch, name)())
        except TypeError: pass