I/O for NIKA2 time ordered data.
Currently supports the calibrated (Jy/beam) IDL .xdr save files produced by
the NIKA2 pipeline (stru_scan<date>s<scan>_<band>_<kind>.xdr).
read_xdr(fname, max_nan_frac=0.1)
Read a calibrated NIKA2 TOD from an IDL .xdr save file.
Parameters:
| Name |
Type |
Description |
Default |
fname
|
str
|
Path to the stru_scan*.xdr file.
|
required
|
max_nan_frac
|
float
|
KIDs with a larger fraction of NaN samples than this are cut.
Remaining NaNs are linearly gapfilled.
|
0.1
|
Returns:
| Name | Type |
Description |
dat |
dict | None
|
Dict that can be passed to minkasi.tods.Tod.
dat_calib is in Jy/beam, dx and dy are RA/Dec in radians.
Returns None if no KIDs survive the cuts.
|
Source code in witch/external/nika2/io.py
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115 | def read_xdr(fname: str, max_nan_frac: float = 0.1) -> Optional[dict[str, Any]]:
"""
Read a calibrated NIKA2 TOD from an IDL ``.xdr`` save file.
Parameters
----------
fname : str
Path to the ``stru_scan*.xdr`` file.
max_nan_frac : float, default: 0.1
KIDs with a larger fraction of NaN samples than this are cut.
Remaining NaNs are linearly gapfilled.
Returns
-------
dat : dict | None
Dict that can be passed to ``minkasi.tods.Tod``.
``dat_calib`` is in Jy/beam, ``dx`` and ``dy`` are RA/Dec in radians.
Returns ``None`` if no KIDs survive the cuts.
"""
m = XDR_RE.search(os.path.basename(fname))
if m is None:
raise ValueError(f"{fname} does not look like a NIKA2 stru_scan xdr file")
date, scan_num, band, kind = m.groups()
raw = readsav(fname)
if "scan" not in raw:
raise ValueError(f"{fname} has no 'scan' structure")
data = raw["scan"]
names = data.dtype.names
units = [u.decode() if isinstance(u, bytes) else str(u) for u in data["UNITS"][0]]
if not any("Jy/beam" in u for u in units):
raise ValueError(f"{fname} is not in Jy/beam, units are {units}")
signal = np.asarray(data["SIGNAL_PLW"][0], dtype="float64")
nan_frac = np.mean(~np.isfinite(signal), axis=-1)
good_dets = np.where(nan_frac <= max_nan_frac)[0]
if len(good_dets) == 0:
print(f"No good detectors in tod {fname}")
return None
ndet = len(good_dets)
ndata = signal.shape[1]
dat: dict[str, Any] = {}
dat["pixid"] = good_dets
for key in ("ARRAY", "ACQBOX", "ACQSUBBAND"):
if key in names:
dat[key.lower()] = np.asarray(data[key][0])[good_dets]
signal = signal[good_dets]
signal[~np.isfinite(signal)] = np.nan
signal, cuts = gapfill_lin(signal)
dat["dat_calib"] = np.ascontiguousarray(signal, dtype="float64")
dat["cuts"] = cuts
ra = np.asarray(data["XRA_BOLOS_PLW"][0], dtype="float64")[good_dets]
dec = np.asarray(data["YDEC_BOLOS_PLW"][0], dtype="float64")[good_dets]
dat["dx"] = np.ascontiguousarray(np.deg2rad(_unwrap_deg(ra)))
dat["dy"] = np.ascontiguousarray(np.deg2rad(dec))
subscans = np.asarray(data["IND_SUBSCANS"][0])
dat["subscans"] = subscans
dat["turnarounds"] = np.where(np.diff(subscans))[0]
dat["dt"] = float(np.median(np.diff(data["TIME_S"][0])))
if "ELEV" in names:
dat["elev"] = float(data["ELEV"][0])
if "TAU_EFF" in names:
dat["tau_eff"] = float(data["TAU_EFF"][0])
dat["fname"] = fname
dat["calinfo"] = {
"calinfo": True,
"date": date,
"scan": scan_num,
"band": band,
"kind": kind,
"units": "Jy/beam",
}
ff = 180 / np.pi
print(
f"ndet {ndet} ndata {ndata} on {fname} with lims "
f"{dat['dx'].min() * ff:.4f} {dat['dx'].max() * ff:.4f} "
f"{dat['dy'].min() * ff:.4f} {dat['dy'].max() * ff:.4f}"
)
return dat
|