Tutorials
For the entirety of this section, it is expected that the user has the latest version of the chime-frb-api installed with proper CHIME/FRB credentials active.
How do you query information about FRBs (including repeaters) from the CHIME/FRB database?¶
The event ID, data path, and signal-to-noise ratio (S/N) of bursts associated with a particular repeater can be queried from chime-frb-api using the following python code (credit: Paul Scholz)
Python
The python code provided above can be used to return the event_id, data_path, and snr of each burst in the CHIME/FRB database associated with a particular TNS. There are other properties that can be returned as well, which can be performed analogously by querying the desired parameters from the event_id object.
from chime_frb_api import frb_master
def get_event_datapath_and_snr(event):
for par in event["measured_parameters"]:
if par["pipeline"]["name"] == "realtime":
event_date = par["datetime"].split(" ")[0].split("-")
event_snr = par["snr"]
data_path = (
"/data/frb-archiver/"
+ event_date[0]
+ "/"
+ event_date[1]
+ "/"
+ event_date[2]
+ "/astro_"
+ str(event["id"])
+ "/"
+ "intensity/processed/"
)
return data_path, event_snr
master = (
frb_master.FRBMaster()
)
repeater_name = "<TNS name of FRB>"
repeater_source = master.sources.get_source(repeater_name)
candidates = master.candidates.get_all_candidates()
event_ids = []
for candidate in candidates:
if candidate["id"] in repeater_source["candidate_id"]:
event_ids += candidate["events"]
for event_id in event_ids:
event = master.events.get_event(event_id)
data_path, snr = get_event_datapath_and_snr(event)
print(event_id, data_path, snr)
How do you get access to the raw (msgpack) and processed CHIME/FRB data?¶
- Data archiving locations on site:
- Raw intensity data (secured using read-only everywhere):
- Processed intensity data:
- Raw baseband data (secured using read-only everywhere):
-
Processed baseband data:
-
Scratch space for users doing other analysis:
- The data may also be available at CANFAR:
Additional details and information can be found at: https://github.com/CHIMEFRB/frb-ops/issues/126
How do you get events that have baseband data?¶
Python
import chime_frb_api, requests
master = chime_frb_api.frb_master.FRBMaster()
master.API.authorize()
auth = {"authorization": master.API.access_token}
r = requests.get(
"https://frb.chimenet.ca/frb-master/v1/verification/get-verifications/NEW SOURCE",
headers=auth,
)
events = r.json()
events_with_bb = []
for event in events:
if event["data_status"]["baseband_data"]:
events_with_bb.append(event["event_id"])
How do you get events classified by tsar as FRBs/RFI/REPEATERS/KNOWN?¶
Replace the "XXX" in the following script to:
KNOWN%20CANDIDATE: to get REPEATERsFAINT: to get Faint sourcesKNOWN%20SOURCE: to get Galactic source(pulsar)New%20CANDIDATE: to get FRBRFI: to get RFI (Note that the XXX is in a URL, which cannot accept spaces. Spaces are replaced with %20).
Python
import chime_frb_api, requests
master = chime_frb_api.frb_master.FRBMaster()
master.API.authorize()
auth = {"authorization": master.API.access_token}
r = requests.get(
"https://frb.chimenet.ca/frb-master/v1/verification/get-verifications/XXX",
headers=auth,
)
# "XXX" is the user classification, which is in the content of: event['user_verification'][-1]['classification']"
events = r.json()
events_selected = []
for event in events:
events_selected.append(event["event_id"])
How to provide publication ready waterfall¶
- Quest waterfall data:
Script: get_wfall.py
Python
import click
import numpy as np
from iautils import cascade, spectra_utils
@click.command()
@click.option("--eventid", type=click.INT)
@click.option("--dm", type=click.FLOAT)
@click.option("--nsub", type=click.INT)
def wfall_maker(eventid, dm, nsub):
event_max = cascade.make_cascade_from_eventid(eventid, beams="max")
event_max.max_beam.intensity = spectra_utils.scale_chime_data(
event_max.max_beam.intensity
)
event_max.dm = dm
event_max.max_beam.subband(nsub)
fname = "{}_wfall.npz".format(eventid)
print("save to:", fname)
np.savez(
fname,
intensity=event_max.max_beam.intensity,
dt=event_max.max_beam.dt,
df=event_max.max_beam.df,
dm=event_max.max_beam.dm,
fbottom=event_max.max_beam.fbottom,
)
if __name__ == "__main__":
wfall_maker()
Example command:
nsub is the number of channels, which could also be specified in the plotting script.The waterfall would by default be saved as EVENTID_wfall.npz
- Generate single waterfall:
Script:plot_intensity_wfalls.py
Python
import matplotlib.pyplot as plt
import numpy as np
from frb_common import common_utils
from iautils import spectra
import argparse
parser=argparse.ArgumentParser()
parser.add_argument("--fname",dest='fname',required=True,metavar="FILE",help="the path and name of the npz file containing waterfall, acquired from get_wfall.py")
parser.add_argument("--dm",dest='dm',type=float,help="the targeted DM, default set to the value in npz file")
parser.add_argument("--figname",dest='figname',metavar="FILE",help="the path and name of the figure file")
args = parser.parse_args()
def plot_wfall(fname,ax,dm=None):
'''
fname: the numpy file contains waterfall
ax: 2 element, 0 for the time profile, 1 for the waterfall
dm: the targeted dm, default is None, it will use the DM from the waterfall file.
'''
data = np.load(fname)
intensity = data["intensity"]
dt = data["dt"]
df = data["df"]
fbottom = data["fbottom"]
dmtrial = data["dm"]
print('DM',dm)
weights = np.ones_like(intensity, dtype=float)
weights[intensity == 0.] = 0.
spec = spectra.Spectra(intensity, common_utils.freq_bottom_mhz,
df, 0.0, dt, weights=weights, dm=dmtrial)
if dm is not None:
spec.dedisperse(dm)
#################################################
###ADD MASK here!!
################################################
#mask = np.argwhere(np.abs(np.nansum(spec.intensity, axis=1)) \
# > 2)
#spec.intensity[mask,...] = np.nan
#mask = np.argwhere(np.abs(np.nanvar(spec.intensity, axis=1)) \
# > 0.05)
#spec.intensity[mask,...] = np.nan
spec.subband(128)
spec.intensity[spec.intensity == 0.] = np.nan
# additional masking
#mask = np.argwhere(np.nanstd(spec.intensity, axis=1) > 0.2)
#spec.intensity[mask,...] = np.nan
ts = np.nansum(spec.intensity, axis=0)
pulse_peak_idx = np.argmax(ts)
window = int(np.ceil(0.025 / dt))
###
if pulse_peak_idx<window:
pulse_peak_idx=window-1
vmin = np.nanpercentile(
spec.intensity[...,pulse_peak_idx-window+1:pulse_peak_idx+window],
#100.-95.)
#100.-99.73)
2.5)
vmax = np.nanpercentile(
spec.intensity[...,pulse_peak_idx-window-40:pulse_peak_idx-window+1],
#95.)
#99.73)
100-2.5)
plot_times = (spec.times[pulse_peak_idx-window+1:pulse_peak_idx+window] \
- spec.center_times[pulse_peak_idx]) * 1000.0 # ms
plot_times = np.append(plot_times, plot_times[-1] + spec.dt)
ax[0].plot(plot_times,
np.append(ts[pulse_peak_idx-window+1:pulse_peak_idx+window],
ts[pulse_peak_idx-window+1:pulse_peak_idx+window][-1]),
lw=1, drawstyle="steps-post", color="tab:blue")
ax[1].imshow(
spec.intensity[...,pulse_peak_idx-window+1:pulse_peak_idx+window],
origin="lower", aspect="auto", interpolation="nearest",
extent=(plot_times.min(), plot_times.max(), fbottom, fbottom+400.),
vmin=vmin, vmax=vmax, cmap="viridis")
ax[1].set_xlabel("Time [ms]", fontsize=9)
xpos = (plot_times.max() - plot_times.min()) * 0.03 + plot_times.min()
ylim = ax[0].get_ylim()
ypos = (ylim[1] - ylim[0]) * 0.93 + ylim[0]
ax[0].set_xlim(plot_times.min(), plot_times.max())
ax[1].set_xlim(plot_times.min(), plot_times.max())
ax[0].set_yticklabels([], visible=False)
ax[0].set_yticks([])
ax[0].set_xticklabels([], visible=False)
ax[1].set_xticks([-20, 0, 20])
ax[1].set_xticklabels([-20, 0, 20], fontsize=9)
ax[1].set_yticks([400, 500, 600, 700, 800])
ax[1].set_yticklabels([400, 500, 600, 700, 800], fontsize=9)
ax[1].set_ylabel("Frequency [MHz]", fontsize=9)
if __name__ == "__main__":
fig, ax = plt.subplots(2, figsize=(4.26, 3),
gridspec_kw={"height_ratios": [1, 3],
"hspace": 0, "wspace": 0.1})
plot_wfall(args.fname,ax,dm=args.dm)
#plot_wfall()
#ax[0].set_title("{}".format(tns_name, fontsize=9))
if args.figname is None:
figname=args.fname.replace('.npz','.pdf')
plt.savefig(figname)
print('waterfall save to:',figname)
Example command:
If additional masks are needed, please modify the corresponding section in the code (marked in the code, starting from Line 42)- Generate figures containing multiple waterfall:
Script:plot_intensity_wfalls_multi.py
Python
import matplotlib.pyplot as plt
import numpy as np
from frb_common import common_utils
from iautils import spectra
if __name__ == "__main__":
fnames = ["69541452_wfall.npz", "101903385_wfall.npz", "145207799_wfall.npz"]
tns_names = ["20200120E", "20200718A", "20201129A"]
dm_goals = [87.76, 87.86, 87.81]
outname = "wfall.pdf"
fig, ax = plt.subplots(
2,
3,
figsize=(4.26, 3),
gridspec_kw={
"height_ratios": [1, 3],
"width_ratios": [1, 1, 1],
"hspace": 0,
"wspace": 0.1,
},
)
for i, fname in enumerate(fnames):
data = np.load(fname)
intensity = data["intensity"]
dt = data["dt"]
df = data["df"]
fbottom = data["fbottom"]
dm = data["dm"]
weights = np.ones_like(intensity, dtype=float)
weights[intensity == 0.0] = 0.0
spec = spectra.Spectra(
intensity, common_utils.freq_bottom_mhz, df, 0.0, dt, weights=weights, dm=dm
)
spec.dedisperse(dm_goals[i])
mask = np.argwhere(np.abs(np.nansum(spec.intensity, axis=1)) > 2)
spec.intensity[mask, ...] = np.nan
mask = np.argwhere(np.abs(np.nanvar(spec.intensity, axis=1)) > 0.05)
spec.intensity[mask, ...] = np.nan
spec.subband(128)
spec.intensity[spec.intensity == 0.0] = np.nan
# additional masking
# mask = np.argwhere(np.nanstd(spec.intensity, axis=1) > 0.2)
# spec.intensity[mask,...] = np.nan
ts = np.nansum(spec.intensity, axis=0)
pulse_peak_idx = np.argmax(ts)
window = int(np.ceil(0.025 / dt))
vmin = np.nanpercentile(
spec.intensity[..., pulse_peak_idx - window + 1 : pulse_peak_idx + window],
# 100.-95.)
# 100.-99.73)
2.5,
)
vmax = np.nanpercentile(
spec.intensity[
..., pulse_peak_idx - window - 40 : pulse_peak_idx - window + 1
],
# 95.)
# 99.73)
100 - 2.5,
)
plot_times = (
spec.times[pulse_peak_idx - window + 1 : pulse_peak_idx + window]
- spec.center_times[pulse_peak_idx]
) * 1000.0 # ms
plot_times = np.append(plot_times, plot_times[-1] + spec.dt)
ax[0, i].plot(
plot_times,
np.append(
ts[pulse_peak_idx - window + 1 : pulse_peak_idx + window],
ts[pulse_peak_idx - window + 1 : pulse_peak_idx + window][-1],
),
lw=1,
drawstyle="steps-post",
color="tab:blue",
)
ax[1, i].imshow(
spec.intensity[..., pulse_peak_idx - window + 1 : pulse_peak_idx + window],
origin="lower",
aspect="auto",
interpolation="nearest",
extent=(plot_times.min(), plot_times.max(), fbottom, fbottom + 400.0),
vmin=vmin,
vmax=vmax,
cmap="viridis",
)
ax[1, i].set_xlabel("Time [ms]", fontsize=9)
xpos = (plot_times.max() - plot_times.min()) * 0.03 + plot_times.min()
ylim = ax[0, i].get_ylim()
ypos = (ylim[1] - ylim[0]) * 0.93 + ylim[0]
ax[0, i].set_title("{}".format(tns_names[i]), fontsize=9)
ax[0, i].set_xlim(plot_times.min(), plot_times.max())
ax[1, i].set_xlim(plot_times.min(), plot_times.max())
ax[0, i].set_yticklabels([], visible=False)
ax[0, i].set_yticks([])
ax[0, i].set_xticklabels([], visible=False)
if i > 0:
plt.setp(ax[1, i].get_yticklabels(), visible=False)
ax[1, i].set_xticks([-20, 0, 20])
ax[1, i].set_xticklabels([-20, 0, 20], fontsize=9)
ax[1, 0].set_yticks([400, 500, 600, 700, 800])
ax[1, 0].set_yticklabels([400, 500, 600, 700, 800], fontsize=9)
ax[1, 0].set_ylabel("Frequency [MHz]", fontsize=9)
fig.suptitle("(a)", fontsize=9, y=1)
plt.savefig(outname)
Modify line 10,11,14,15 for the waterfall file name, dm and output image name.