Plot data files from https://lpsc.in2p3.fr/cosmic-rays-db without unpacking the tarball

plot_cosmic_rays.py

import tarfile
import pandas as pd
from io import TextIOWrapper
import matplotlib.pyplot as plt
import os
import argparse


parser = argparse.ArgumentParser()
parser.add_argument('outputfile')
args = parser.parse_args()


def read_data(tar, member):
    f = TextIOWrapper(tar.extractfile(member), encoding='ascii')

    for i in range(2):
        line = next(f)

    names = line.lstrip('#').strip().replace('<', '').replace('>', '').split()

    return pd.read_csv(f, sep=r'\s+', comment='#', names=names, header=None)


f = tarfile.open(f'proton_spectra.tar.gz')

experiments = read_data(f, 'data_exps.txt')
experiments['file_name'] = experiments['file_name'].apply(os.path.basename)

# they fucked up the column names
if 'qty' in experiments.columns:
    experiments[['qty', 'exp_name']] = experiments[['exp_name', 'qty']]

data = {}
for row in experiments.itertuples():
    df = read_data(f, row.file_name)
    experiments.loc[row.Index, 'e_min'] = df.E.min()
    experiments.loc[row.Index, 'e_max'] = df.E.max()
    data[row.exp_name] = df

experiments['file_num'] = experiments.index + 1
experiments['experiment'] = experiments['exp_name'].str.split('(', expand=True)[0]

# select 20 experiments with the highest energy
experiments = experiments.sort_values('e_max').tail(20)
experiments.set_index('experiment', inplace=True)

names = ['AMS02', 'ATIC02', 'CREAM-I', 'PAMELA', 'JACEE', 'RUNJOB']
all_data = pd.DataFrame()

for i, name in enumerate(names):
    exp = experiments.loc[name]

    df = read_data(f, f'data_exp{exp.file_num}.dat')
    df['experiment'] = name
    all_data = all_data.append(df)

    plt.errorbar(
        df.E,
        df.y * df.E**2,
        yerr=(
            df.yerrtot_lo * df.E**2,
            df.yerrtot_up * df.E**2,
        ),
        linestyle='',
        marker='.',
        label=name,
        ms=5,
        mew=0,
    )


plt.grid()
plt.legend(ncol=3, bbox_to_anchor=[0.5, 1.025], loc='lower center')
plt.yscale('log')
plt.xscale('log')
plt.xlabel(r'$E \,\,/\,\, \mathrm{GeV}$')

plt.ylabel(
    r'$E^2 \frac{\mathrm{d}N}{\mathrm{d}E} /'
    r' (\mathrm{GeV}^{-1}\, \mathrm{m}^{-2}\, \mathrm{s}^-1 \mathrm{sr}^{-1})$'
)

plt.tight_layout(pad=0)
plt.savefig(args.outputfile)