Repository navigation
Expand file tree
/
Copy patherosita_dr2.py
More file actions
687 lines (570 loc) · 31.8 KB
/
Copy patherosita_dr2.py
File metadata and controls
687 lines (570 loc) · 31.8 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
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
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
"""eROSITA-DE DR2 (eRASS:3) cross-match for the MCSNRcat catalog.
DR2 (released 2026-07-31) is a **catalogue-only** release — no event lists, no
images, no exposure maps — so it changes nothing about the imagery this site
streams. What it adds is a single uniform X-ray source list, twice as deep as
eRASS1, over exactly the hemisphere the Magellanic Clouds sit in: 1 975 540
sources at ``DET_LIKE_0 >= 5`` in 0.2-2.3 keV (Ramos-Ceja et al. 2026).
For this catalog that buys three things the existing X-ray block cannot give:
* **completeness** — every object gets a detection *or* a quantified 3σ flux
upper limit, so an empty cell stops meaning "nobody published one";
* **uniformity** — one pipeline applied to all 124 positions, where the
``xray_rate_ctss`` / ``lx_1e35`` / ``fx_z26_*`` columns are three papers with
three selection functions;
* **independence** — Maggi+16 and Zangrandi+24/26 are this catalog's own
sources; DR2 is a different team, pipeline and data cut.
**What it is not.** A detection is not a confirmation. The classification rule
(README: two of three criteria, one of them *diffuse* X-ray emission) is about
the remnant's own emission, while an eRASS:3 entry may be any point source that
happens to lie in the aperture — and at eRASS:3 depth the LMC field carries
~250 sources/deg², so "there is an X-ray source here" is nearly free. Hence
``erass3_p_chance`` (the Poisson chance of a coincidental match at that
aperture) is written next to every detection, and **no ``erass3_*`` column may
ever change an object's ``klass``**: ``merge`` asserts it and
``tests/test_erosita_dr2.py`` tests it.
Usage
-----
Fetch the cross-match (network; results cached under ``data/raw/erosita_dr2``,
so a re-run is offline apart from the upper limits)::
python erosita_dr2.py fetch
Then fold the ``erass3_*`` columns into a **new** catalog version::
python catalog.py new-version
python erosita_dr2.py merge
Endpoints, verbosity levels and the two-stage fetch are documented inline
below; the same services are described in ``../VLMism/docs/DATA_SOURCING.md``
§1b, whose ``vlmism/archives/erosita_cat.py`` this is a dependency-free port of
(stdlib HTTP + XML, numpy geometry — so it runs in the same environment as
``build.py``, with no ``requests``/``astropy``).
"""
from __future__ import annotations
import argparse
import json
import logging
import subprocess
import urllib.error
import urllib.parse
import urllib.request
from dataclasses import dataclass
from datetime import datetime, timezone
from pathlib import Path
from xml.etree import ElementTree
import numpy as np
import pandas as pd
from build import latest_catalog
from catalog import write_catalog
logger = logging.getLogger(__name__)
DATA_DIR = Path(__file__).parent / "data"
#: Cross-match product: one row per catalog object, committed so ``merge`` (and
#: anyone auditing the numbers) needs no network.
CROSSMATCH_CSV = DATA_DIR / "erosita_dr2_erass3.csv"
#: Raw per-cone service responses (gitignored) — delete to force a re-fetch.
CACHE_DIR = DATA_DIR / "raw" / "erosita_dr2"
ERODAT_BASE = "https://erosita.mpe.mpg.de/erodat"
#: IVOA Simple Cone Search over the released source catalogues.
SCS_URL = f"{ERODAT_BASE}/catalogue/SCS"
#: Upper-flux-limit server; the ``_multi`` endpoint takes many positions per POST.
UPPER_LIMIT_URL = f"{ERODAT_BASE}/upperlimit/service_multi"
#: DR2_Main is the 0.2-2.3 keV detection catalogue (the band this catalog's
#: other eROSITA columns are measured in). DR2_Hard is the 2.3-5.0 keV run.
CATALOGUE = "DR2_Main"
#: Service-side cap on a single cone-search radius, degrees.
MAX_CONE_RADIUS_DEG = 5.0
#: The sample spans ~9 deg on the sky, so the sweep is tiled into legal cones.
#: The pad only has to exceed the largest match aperture (5.05' here).
SWEEP_PAD_ARCMIN = 20.0
#: Pad on each per-object detail cone, so a detection sitting exactly on the
#: aperture edge is still returned by the second-stage query.
DETAIL_PAD_ARCMIN = 1.0
#: Floor on the per-object aperture. Below the ~26" (0.43') eROSITA survey PSF
#: a detection centroid could not land inside the aperture even in principle.
#: (No catalog object is currently this small — the smallest r_arcmin is 0.65 —
#: but the floor keeps the join honest if a compact object is ever added.)
MIN_APERTURE_ARCMIN = 0.5
#: Outer radius of the annulus used to measure the local source density around
#: each object, and from it the chance of a coincidental match. Apertures span
#: 0.65'-5.05' — a ~60x range in area — so a raw detection rate would ask a
#: different question of each object without this.
DENSITY_ANNULUS_ARCMIN = 15.0
#: Upper-limit server settings. 024 = 0.2-2.3 keV, matching CATALOGUE's band.
UL_BAND = "024"
UL_SURVEY = "DR2_eRASSc3"
#: Service cap is 10 000 positions per POST; 124 objects fit in one request.
UL_BATCH = 10_000
#: Network timeout, seconds. A 5 deg cone answers in ~10 s; the ceiling is
#: generous because the archive occasionally queues.
TIMEOUT_S = 300.0
#: NASA ADS bibcodes stamped into the cross-match product.
DR2_CATALOGUE_REF = "2026A&A...712A.171R" # Ramos-Ceja et al. 2026, A&A 712, A171
UPPER_LIMIT_REF = "2024A&A...682A..35T" # Tubín-Arenas et al. 2024, A&A 682, A35
#: Columns VERB=1 returns — the cheap sweep. Measured against the live service
#: (2026-09-05): ~262 bytes/row at VERB=1 against ~6.4 kB/row at VERB=3, so
#: sweeping the whole LMC footprint at full verbosity would move ~225 MB to
#: gain one flag column. Hence the two-stage fetch: sweep wide at VERB=1, then
#: re-query at VERB=3 only the few small apertures that actually matched.
SWEEP_COLUMNS = (
"iauname", "detuid", "ra", "dec", "ext", "ext_like", "det_like_0",
"ml_rate_1", "ml_flux_1",
)
#: The only VERB=3 column this catalog keeps. ``flag_sp_lga`` is a
#: *spuriousness* flag — "source might lie within an overdense region near a
#: local large galaxy" — and the LMC *is* that galaxy, so it fires on ~18% of
#: detections here and means "possible blend or artefact of the galaxy's own
#: diffuse emission", never "belongs to the LMC".
#:
#: Deliberately not kept: ``flag_sp_snr``. Its overdensities are Galactic and
#: it is 0 for every LMC detection — a column of zeros named ``snr`` on a
#: supernova-remnant catalog is a trap, not data.
DETAIL_COLUMNS = ("detuid", "flag_sp_lga")
#: The block written into ``lmc_snrs_extended_v<N>.csv`` by ``merge``, in page
#: order. Nothing outside this list is copied into the catalog.
ERASS3_COLUMNS = (
"erass3_name",
"erass3_detuid",
"erass3_sep_arcmin",
"erass3_n_in_aperture",
"erass3_n_extended",
"erass3_det_like",
"erass3_ext_arcsec",
"erass3_ext_like",
"erass3_flux",
"erass3_flux_ul",
"erass3_flux_ul_src",
"erass3_p_chance",
"erass3_flag_lga",
)
#: Identity/classification columns ``merge`` must leave byte-identical. A DR2
#: detection is not a confirmation (see the module docstring); the census is
#: set by the literature, never by a cross-match.
PROTECTED_COLUMNS = ("snr_key", "id", "name", "klass")
class ServiceError(RuntimeError):
"""The eRODat service answered, but not with a usable catalogue."""
@dataclass(frozen=True)
class Cone:
"""A cone-search footprint, in degrees."""
ra: float
dec: float
radius_deg: float
def key(self) -> str:
"""Filename-safe cache key (also the human-readable footprint)."""
return f"ra{self.ra:+09.4f}_dec{self.dec:+09.4f}_r{self.radius_deg:.4f}"
# --------------------------------------------------------------------------
# geometry (numpy only — the sample is 124 objects, no spatial index needed)
# --------------------------------------------------------------------------
def _unit_vectors(ra_deg: np.ndarray, dec_deg: np.ndarray) -> np.ndarray:
"""Positions as unit vectors, shape ``(n, 3)``."""
ra = np.radians(np.asarray(ra_deg, dtype=float))
dec = np.radians(np.asarray(dec_deg, dtype=float))
return np.stack([np.cos(dec) * np.cos(ra), np.cos(dec) * np.sin(ra), np.sin(dec)], axis=-1)
def separation_deg(ra1: float, dec1: float,
ra2: np.ndarray, dec2: np.ndarray) -> np.ndarray:
"""Angular separation of one position from many, in degrees.
Haversine rather than ``arccos`` of a dot product: the separations that
matter here are arcminutes, where the dot-product form loses precision.
"""
lat1, lat2 = np.radians(dec1), np.radians(np.asarray(dec2, dtype=float))
dlat = lat2 - lat1
dlon = np.radians(np.asarray(ra2, dtype=float) - ra1)
hav = np.sin(dlat / 2.0) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2.0) ** 2
return np.degrees(2.0 * np.arcsin(np.sqrt(np.clip(hav, 0.0, 1.0))))
def cover_positions(ra: np.ndarray, dec: np.ndarray, *,
max_radius_deg: float = MAX_CONE_RADIUS_DEG,
pad_arcmin: float = SWEEP_PAD_ARCMIN) -> list[Cone]:
"""Tile a set of positions into as few legal cone searches as possible.
Deterministic greedy cover: take the first uncovered position, claim
everything within ``max_radius_deg`` of it, re-centre on the claimed set to
shrink the cone, repeat. Ported from VLMism's ``erosita_cat`` so the two
repos ask the archive the same question.
"""
ra = np.atleast_1d(np.asarray(ra, dtype=float))
dec = np.atleast_1d(np.asarray(dec, dtype=float))
if ra.size == 0:
raise ValueError("no positions to cover")
if max_radius_deg > MAX_CONE_RADIUS_DEG:
raise ValueError(f"cone radius {max_radius_deg} deg exceeds the eRODat cap")
pad_deg = pad_arcmin / 60.0
claim_radius = max(max_radius_deg - pad_deg, pad_deg)
vectors = _unit_vectors(ra, dec)
remaining = np.ones(ra.size, dtype=bool)
cones: list[Cone] = []
while remaining.any():
seed = int(np.flatnonzero(remaining)[0])
claimed = remaining & (separation_deg(ra[seed], dec[seed], ra, dec) <= claim_radius)
centroid = vectors[claimed].mean(axis=0)
norm = float(np.linalg.norm(centroid))
if norm == 0.0: # antipodal set; cannot happen for one galaxy, but be explicit
raise ValueError("positions have no meaningful centroid")
centroid /= norm
c_dec = float(np.degrees(np.arcsin(centroid[2])))
c_ra = float(np.degrees(np.arctan2(centroid[1], centroid[0])) % 360.0)
radius = float(separation_deg(c_ra, c_dec, ra[claimed], dec[claimed]).max()) + pad_deg
cones.append(Cone(c_ra, c_dec, min(max(radius, pad_deg), max_radius_deg)))
remaining &= ~claimed
logger.info("covering %d positions with %d cone(s): %s", ra.size, len(cones),
", ".join(f"{c.radius_deg:.2f} deg" for c in cones))
return cones
# --------------------------------------------------------------------------
# eRODat services
# --------------------------------------------------------------------------
def parse_votable(payload: bytes) -> pd.DataFrame:
"""Parse an SCS VOTable response into a DataFrame.
The service emits VOTable 1.0 in TABLEDATA form with no XML namespace, but
tags are matched by local name anyway so a namespaced future response still
parses. Values arrive as text; numeric columns are coerced by the caller's
use (``pd.to_numeric``), except ``ra``/``dec``, which are needed as floats
everywhere and so are converted here.
"""
try:
root = ElementTree.fromstring(payload)
except ElementTree.ParseError as err:
head = payload[:200].decode("utf-8", errors="replace")
raise ServiceError(f"eRODat response is not parseable XML: {err}; body starts {head!r}") from err
def local(elem: ElementTree.Element) -> str:
return elem.tag.rsplit("}", 1)[-1]
fields = [e.get("name") for e in root.iter() if local(e) == "FIELD"]
if not fields:
head = payload[:400].decode("utf-8", errors="replace")
raise ServiceError(f"eRODat response carries no FIELD definitions; body starts {head!r}")
rows = []
for tr in (e for e in root.iter() if local(e) == "TR"):
cells = [(td.text or "").strip() for td in tr if local(td) == "TD"]
if len(cells) != len(fields):
raise ServiceError(f"row has {len(cells)} cells, header declares {len(fields)} fields")
rows.append(cells)
frame = pd.DataFrame(rows, columns=fields)
for column in ("ra", "dec"):
if column in frame.columns:
frame[column] = pd.to_numeric(frame[column], errors="coerce")
return frame
def cone_search(cone: Cone, *, verbosity: int, cache_dir: Path = CACHE_DIR,
refresh: bool = False, timeout: float = TIMEOUT_S) -> pd.DataFrame:
"""One cone query against DR2, cached to ``cache_dir`` as CSV + provenance.
An empty cone is a legitimate answer here (a small detail cone can come
back empty if the archive re-runs source detection), so it is returned as
an empty frame rather than raised.
"""
cache_dir.mkdir(parents=True, exist_ok=True)
stem = f"{CATALOGUE}_verb{verbosity}_{cone.key()}"
csv_path = cache_dir / f"{stem}.csv"
if csv_path.exists() and not refresh:
return pd.read_csv(csv_path, dtype={"iauname": str, "detuid": str})
params = {"CAT": CATALOGUE, "RA": f"{cone.ra:.6f}", "DEC": f"{cone.dec:.6f}",
"SR": f"{cone.radius_deg:.6f}", "VERB": str(verbosity)}
url = f"{SCS_URL}?{urllib.parse.urlencode(params)}"
logger.info("cone search: %s VERB=%d at (%.4f, %.4f) r=%.4f deg",
CATALOGUE, verbosity, cone.ra, cone.dec, cone.radius_deg)
try:
with urllib.request.urlopen(url, timeout=timeout) as response:
payload = response.read()
except urllib.error.URLError as err:
raise ServiceError(f"eRODat cone search failed ({url}): {err}") from err
frame = parse_votable(payload)
logger.info(" -> %d source(s)", len(frame))
frame.to_csv(csv_path, index=False)
(cache_dir / f"{stem}.meta.json").write_text(json.dumps({
"url": url, "catalogue": CATALOGUE, "verbosity": verbosity,
"ra": cone.ra, "dec": cone.dec, "radius_deg": cone.radius_deg,
"n_rows": len(frame), "ref": DR2_CATALOGUE_REF, "fetched_utc": _now_utc(),
}, indent=2), encoding="utf-8")
return frame
def upper_limits(ra: np.ndarray, dec: np.ndarray, *, band: str = UL_BAND,
survey: str = UL_SURVEY, batch_size: int = UL_BATCH,
timeout: float = TIMEOUT_S) -> pd.DataFrame:
"""3σ flux upper limits at arbitrary positions, from the eROSITA UL server.
This is what turns "not in the catalogue" into a number. Limits are
one-sided 3σ (99.87% CL) for an **absorbed power law with Γ = 2.0 and
N_H = 3e20 cm^-2** (Tubín-Arenas et al. 2024) — not a thermal remnant
model, so for an SNR the limit is indicative rather than exact. That
caveat travels with the column onto the object pages.
Returns one row per position, in input order, with ``UL_B`` (upper limit),
``UL_S`` (same, treating any detected source at the position as
background — the right one when the target is *not* the detected source),
``Flag_pos`` and ``Exposure``.
"""
ra = np.atleast_1d(np.asarray(ra, dtype=float))
dec = np.atleast_1d(np.asarray(dec, dtype=float))
if ra.size == 0:
raise ValueError("no positions given")
if ra.shape != dec.shape:
raise ValueError(f"ra/dec shape mismatch: {ra.shape} vs {dec.shape}")
records: list[dict] = []
for start in range(0, ra.size, batch_size):
chunk = [{"ra": float(r), "dec": float(d), "band": band, "dr_survey": survey}
for r, d in zip(ra[start:start + batch_size], dec[start:start + batch_size])]
logger.info("upper limits: %s band %s, %d position(s)", survey, band, len(chunk))
request = urllib.request.Request(
UPPER_LIMIT_URL, data=json.dumps(chunk).encode("utf-8"),
headers={"Content-Type": "application/json"}, method="POST")
try:
with urllib.request.urlopen(request, timeout=timeout) as response:
body = json.loads(response.read())
except (urllib.error.URLError, json.JSONDecodeError) as err:
raise ServiceError(f"upper-limit request failed: {err}") from err
error = str(body.get("error", "None"))
if error not in {"None", "", "null"}:
raise ServiceError(f"upper-limit server error: {error}")
limits = body.get("limits")
if limits is None or len(limits) != len(chunk):
raise ServiceError(f"upper-limit server returned {0 if limits is None else len(limits)} "
f"row(s) for {len(chunk)} position(s)")
records.extend(limits)
return pd.DataFrame.from_records(records)
# --------------------------------------------------------------------------
# cross-match
# --------------------------------------------------------------------------
def _numeric(frame: pd.DataFrame, column: str) -> np.ndarray:
"""Column as a float array, or all-NaN when absent (VERB=1 responses)."""
if column not in frame.columns:
return np.full(len(frame), np.nan)
return pd.to_numeric(frame[column], errors="coerce").to_numpy(dtype=float)
def _round_sig(value: float, digits: int = 4) -> float:
"""Round to significant digits, so the CSV reads as a measurement.
The services return full float64 repr (``2.007534361279991`` arcmin for a
separation good to ~1"); keeping that would imply precision nobody has.
"""
if not np.isfinite(value) or value == 0.0:
return float(value)
return float(np.round(value, digits - 1 - int(np.floor(np.log10(abs(value))))))
def apertures(objects: pd.DataFrame) -> np.ndarray:
"""Per-object match aperture in arcmin: its own radius, floored.
Not a fixed radius: the catalog spans 0.65'-5.05' and a fixed aperture
would ask a different question at each end — an eROSITA centroid outside a
remnant's shell is a different source, not a detection of it.
"""
radius = pd.to_numeric(objects["r_arcmin"], errors="coerce").to_numpy(dtype=float)
if not np.isfinite(radius).all():
missing = objects.loc[~np.isfinite(radius), "id"].tolist()
raise ValueError(f"objects without r_arcmin cannot be matched: {missing}")
return np.maximum(radius, MIN_APERTURE_ARCMIN)
def match_objects(objects: pd.DataFrame, sources: pd.DataFrame) -> pd.DataFrame:
"""Cross-match catalog objects against eRASS:3 detections.
Each object is matched inside its own aperture (see :func:`apertures`);
every detection in it is counted, the nearest one's properties are
reported, and the local source density measured in a surrounding annulus
gives the Poisson chance that a detection landed there by coincidence.
Objects with nothing inside the aperture carry 0/NaN and an empty name —
that is a measurement, not a failure; :func:`upper_limits` then puts a
number on it.
"""
aperture = apertures(objects)
src_ra = sources["ra"].to_numpy(dtype=float)
src_dec = sources["dec"].to_numpy(dtype=float)
ext_like = _numeric(sources, "ext_like")
flag_lga = _numeric(sources, "flag_sp_lga")
annulus = max(DENSITY_ANNULUS_ARCMIN, float(aperture.max()))
rows: list[dict] = []
for i, (_, obj) in enumerate(objects.iterrows()):
sep = separation_deg(float(obj["ra"]), float(obj["dec"]), src_ra, src_dec) * 60.0
inside = sep <= aperture[i]
in_annulus = (sep <= annulus) & ~inside
# Local density from the annulus only, so an object's own detections
# do not inflate the chance rate they are being compared against.
annulus_area = np.pi * (annulus ** 2 - aperture[i] ** 2) / 3600.0 # deg^2
aperture_area = np.pi * (aperture[i] / 60.0) ** 2 # deg^2
density = in_annulus.sum() / annulus_area if annulus_area > 0 else 0.0
p_chance = 1.0 - np.exp(-density * aperture_area)
row = {
"snr_key": obj["snr_key"], "id": obj["id"], "klass": obj["klass"],
"ra": float(obj["ra"]), "dec": float(obj["dec"]),
"aperture_arcmin": _round_sig(float(aperture[i]), 3),
"erass3_n_in_aperture": int(inside.sum()),
"erass3_n_extended": int((ext_like[inside] > 0).sum()),
"erass3_p_chance": _round_sig(float(p_chance), 3),
"erass3_flag_lga": bool((flag_lga[inside] > 0).any()),
"local_density_deg2": _round_sig(float(density), 4),
}
if inside.any():
nearest = int(np.flatnonzero(inside)[np.argmin(sep[inside])])
near = sources.iloc[nearest]
row.update({
"erass3_name": str(near.get("iauname", "")),
"erass3_detuid": str(near.get("detuid", "")),
"erass3_sep_arcmin": _round_sig(float(sep[nearest]), 3),
"erass3_det_like": _round_sig(float(_numeric(sources, "det_like_0")[nearest]), 4),
"erass3_ext_arcsec": _round_sig(float(_numeric(sources, "ext")[nearest]), 3),
"erass3_ext_like": _round_sig(float(_numeric(sources, "ext_like")[nearest]), 3),
"erass3_flux": _round_sig(float(_numeric(sources, "ml_flux_1")[nearest]), 4),
})
else:
row.update({"erass3_name": "", "erass3_detuid": "", "erass3_sep_arcmin": np.nan,
"erass3_det_like": np.nan, "erass3_ext_arcsec": np.nan,
"erass3_ext_like": np.nan, "erass3_flux": np.nan})
rows.append(row)
frame = pd.DataFrame(rows)
logger.info("matched %d object(s): %d with >=1 detection in aperture "
"(%.1f expected by chance), %d with an extended detection",
len(frame), int((frame["erass3_n_in_aperture"] > 0).sum()),
float(frame["erass3_p_chance"].sum()),
int((frame["erass3_n_extended"] > 0).sum()))
return frame
def attach_upper_limits(match: pd.DataFrame, **kwargs) -> pd.DataFrame:
"""Fill ``erass3_flux_ul`` / ``erass3_flux_ul_src`` for undetected objects.
Detected objects are deliberately not queried: the server would answer
``Flag_pos = 1`` and an unusable limit, and the page shows a flux for them
anyway — an object gets a flux *or* a limit, never both.
"""
out = match.copy()
for column in ("erass3_flux_ul", "erass3_flux_ul_src", "ul_exposure_s"):
out[column] = np.nan
undetected = out["erass3_n_in_aperture"] == 0
if not undetected.any():
logger.info("every object has a detection in aperture; no upper limits needed")
return out
logger.info("querying upper limits for %d undetected object(s)", int(undetected.sum()))
limits = upper_limits(out.loc[undetected, "ra"].to_numpy(dtype=float),
out.loc[undetected, "dec"].to_numpy(dtype=float), **kwargs)
out.loc[undetected, "erass3_flux_ul"] = [
_round_sig(v, 4) for v in pd.to_numeric(limits["UL_B"], errors="coerce")]
out.loc[undetected, "erass3_flux_ul_src"] = [
_round_sig(v, 4) for v in pd.to_numeric(limits["UL_S"], errors="coerce")]
out.loc[undetected, "ul_exposure_s"] = [
_round_sig(v, 4) for v in pd.to_numeric(limits["Exposure"], errors="coerce")]
return out
# --------------------------------------------------------------------------
# commands
# --------------------------------------------------------------------------
def _now_utc() -> str:
return datetime.now(timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ")
def _sweep(objects: pd.DataFrame, *, cache_dir: Path, refresh: bool) -> pd.DataFrame:
"""Stage 1: VERB=1 over the whole footprint, de-duplicated on ``detuid``."""
cones = cover_positions(objects["ra"].to_numpy(dtype=float),
objects["dec"].to_numpy(dtype=float))
frames = [cone_search(c, verbosity=1, cache_dir=cache_dir, refresh=refresh) for c in cones]
sources = pd.concat(frames, ignore_index=True).drop_duplicates(subset="detuid")
logger.info("sweep: %d unique detection(s) over %d cone(s)", len(sources), len(cones))
return sources
def _detail(objects: pd.DataFrame, match: pd.DataFrame, sources: pd.DataFrame, *,
cache_dir: Path, refresh: bool) -> pd.DataFrame:
"""Stage 2: VERB=3 in one small cone per matched object, for ``flag_sp_lga``.
Only the objects that matched something are queried, and only their own
aperture (plus :data:`DETAIL_PAD_ARCMIN`) — the flag is the single column
worth the ~25x bigger response, and only for detections already in hand.
"""
aperture = apertures(objects)
matched = match["erass3_n_in_aperture"].to_numpy() > 0
cones: dict[str, Cone] = {} # keyed so two objects at one position query once
for i in np.flatnonzero(matched):
obj = objects.iloc[int(i)]
cone = Cone(float(obj["ra"]), float(obj["dec"]),
(aperture[i] + DETAIL_PAD_ARCMIN) / 60.0)
cones.setdefault(cone.key(), cone)
logger.info("detail pass: %d cone(s) at VERB=3", len(cones))
frames = [cone_search(c, verbosity=3, cache_dir=cache_dir, refresh=refresh)
for c in cones.values()]
if not frames:
return sources
detail = pd.concat(frames, ignore_index=True).drop_duplicates(subset="detuid")
keep = [c for c in DETAIL_COLUMNS if c in detail.columns]
if "flag_sp_lga" not in keep:
raise ServiceError("VERB=3 response carries no flag_sp_lga column — the "
"release schema changed; check the DR2 column list")
merged = sources.merge(detail[keep], on="detuid", how="left")
logger.info("detail pass: flag_sp_lga known for %d of %d detection(s), %d flagged",
int(merged["flag_sp_lga"].notna().sum()), len(merged),
int((_numeric(merged, "flag_sp_lga") > 0).sum()))
return merged
def cmd_fetch(args: argparse.Namespace) -> None:
"""Query eRODat and write the per-object cross-match CSV."""
catalog_path = Path(args.catalog) if args.catalog else latest_catalog(DATA_DIR)
objects = pd.read_csv(catalog_path)
logger.info("cross-matching %d object(s) from %s", len(objects), catalog_path.name)
cache_dir = Path(args.cache)
sources = _sweep(objects, cache_dir=cache_dir, refresh=args.refresh)
# First pass finds *which* objects matched; the detail pass then adds
# flag_sp_lga for exactly those, and the match is redone with it in hand.
match = match_objects(objects, sources)
if not args.no_detail:
sources = _detail(objects, match, sources, cache_dir=cache_dir, refresh=args.refresh)
match = match_objects(objects, sources)
if not args.no_upper_limits:
match = attach_upper_limits(match)
else:
for column in ("erass3_flux_ul", "erass3_flux_ul_src", "ul_exposure_s"):
match[column] = np.nan
match["ref_catalogue"] = DR2_CATALOGUE_REF
match["ref_upper_limit"] = UPPER_LIMIT_REF
# Terse on purpose — repeated on all 124 rows. The full endpoints live in
# this module and in each cached response's .meta.json.
match["provenance"] = f"{CATALOGUE}/eRASS:3 (eRODat SCS + UL server {UL_SURVEY} band {UL_BAND})"
match["source_catalog"] = catalog_path.name
match["run_utc"] = _now_utc()
ordered = (["snr_key", "id", "klass", "ra", "dec", "aperture_arcmin"]
+ list(ERASS3_COLUMNS)
+ ["local_density_deg2", "ul_exposure_s", "ref_catalogue",
"ref_upper_limit", "provenance", "source_catalog", "run_utc"])
out = Path(args.out)
match[ordered].to_csv(out, index=False)
detected = int((match["erass3_n_in_aperture"] > 0).sum())
logger.info("wrote %s: %d object(s), %d detected, %d with an upper limit",
out, len(match), detected, len(match) - detected)
def _is_tracked(path: Path) -> bool:
"""True when git already tracks ``path`` (i.e. the version is published).
Used only to protect the repo's own rule — never edit an already-published
CSV version, cut a new one. A missing git is not an error here.
"""
repo = Path(__file__).resolve().parent
try:
done = subprocess.run(
["git", "ls-files", "--error-unmatch", "--", str(path.resolve())],
cwd=repo, capture_output=True, check=False)
except OSError:
logger.warning("git not available; skipping the published-version check")
return False
return done.returncode == 0
def cmd_merge(args: argparse.Namespace) -> None:
"""Fold the ``erass3_*`` columns into the latest catalog version."""
crossmatch = pd.read_csv(args.crossmatch,
dtype={"snr_key": str, "erass3_name": str, "erass3_detuid": str})
catalog_path = Path(args.catalog) if args.catalog else latest_catalog(DATA_DIR)
if _is_tracked(catalog_path) and not args.in_place:
raise SystemExit(
f"{catalog_path.name} is already committed — this repo never edits a "
"published version. Run `python catalog.py new-version` first, or pass "
"--in-place if you really mean to rewrite this file.")
# dtype=object, as catalog.py's `set` does: every column the merge does not
# touch is written back byte-for-byte, so the version diff is the eRASS:3
# block and nothing else. Round-tripping through float64 would rewrite
# every integer-valued cell ("6" -> "6.0") for no reason.
df = pd.read_csv(catalog_path, dtype=object)
missing = sorted(set(df["snr_key"]) - set(crossmatch["snr_key"]))
if missing:
raise SystemExit(f"{len(missing)} object(s) absent from the cross-match "
f"(re-run `fetch` against this version): {missing}")
before = df[list(PROTECTED_COLUMNS)].copy()
for column in ERASS3_COLUMNS:
df[column] = df["snr_key"].map(crossmatch.set_index("snr_key")[column])
# Trap 1, enforced rather than promised: a cross-match may add columns, and
# may not touch identity or classification. A DR2 detection can be any
# point source in the aperture; the confirmed/candidate census is set by
# the literature. Also tested in tests/test_erosita_dr2.py.
changed = [c for c in PROTECTED_COLUMNS if not df[c].equals(before[c])]
if changed:
raise SystemExit(f"BUG: merge altered protected column(s) {changed} — refusing to write")
write_catalog(df, catalog_path)
detected = int(pd.to_numeric(df["erass3_n_in_aperture"], errors="coerce").fillna(0).gt(0).sum())
logger.info("merged %d eRASS:3 column(s) into %s: %d detection(s), %d upper limit(s)",
len(ERASS3_COLUMNS), catalog_path.name, detected,
int(df["erass3_flux_ul"].notna().sum()))
def main() -> None:
logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")
parser = argparse.ArgumentParser(
description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
sub = parser.add_subparsers(dest="cmd", required=True)
pf = sub.add_parser("fetch", help="query eRODat and write the cross-match CSV")
pf.add_argument("--catalog", default=None, help="catalog CSV (default: latest in data/)")
pf.add_argument("--out", default=str(CROSSMATCH_CSV), help=f"output CSV (default: {CROSSMATCH_CSV})")
pf.add_argument("--cache", default=str(CACHE_DIR), help="raw response cache directory")
pf.add_argument("--refresh", action="store_true", help="re-query even when cached")
pf.add_argument("--no-detail", action="store_true",
help="skip the VERB=3 pass (leaves erass3_flag_lga all-False)")
pf.add_argument("--no-upper-limits", action="store_true",
help="skip the upper-limit server (leaves erass3_flux_ul empty)")
pf.set_defaults(func=cmd_fetch)
pm = sub.add_parser("merge", help="write the erass3_* columns into the latest catalog CSV")
pm.add_argument("--crossmatch", default=str(CROSSMATCH_CSV), help="cross-match CSV from `fetch`")
pm.add_argument("--catalog", default=None, help="catalog CSV (default: latest in data/)")
pm.add_argument("--in-place", action="store_true",
help="allow rewriting an already-committed catalog version")
pm.set_defaults(func=cmd_merge)
args = parser.parse_args()
args.func(args)
if __name__ == "__main__":
main()