laurence
c6299f41ab
Preparing to merge this repository into a combined astrophotography repo. session-scripts/ becomes pipeline/ because the scripts import layout.py from their own directory and must stay together, and because 'pipeline' says what it is rather than how it came about. observing/ stays at the top level: observing plans are not processing code.
75 lines
2.7 KiB
Python
75 lines
2.7 KiB
Python
"""Choose the background mesh size EMPIRICALLY, using the 589 SIMBAD-catalogued
|
|
Cen A globular clusters that fall in the field as a truth set. For each mesh
|
|
size we measure (a) how many known GCs are recovered, split by projected radius,
|
|
and (b) how many total detections there are (a proxy for spurious detections)."""
|
|
import os, numpy as np, sep, warnings
|
|
from astropy.io import fits
|
|
from astropy.wcs import WCS
|
|
|
|
import layout
|
|
warnings.filterwarnings("ignore")
|
|
S = layout.SESSION
|
|
NUC_RA, NUC_DEC = 201.365063, -43.019113
|
|
PIXSCALE = 0.5376
|
|
|
|
|
|
def gk(fwhm=5.0):
|
|
sig = fwhm / 2.3548
|
|
n = int(2 * round(3 * sig) + 1)
|
|
r = np.arange(n) - n // 2
|
|
xx, yy = np.meshgrid(r, r)
|
|
k = np.exp(-(xx ** 2 + yy ** 2) / (2 * sig ** 2))
|
|
return (k / k.sum()).astype(np.float32)
|
|
|
|
|
|
def match(ra1, d1, ra2, d2, tol_as):
|
|
"""brute-force nearest match, small catalogues"""
|
|
cd = np.cos(np.radians(d1.mean()))
|
|
idx = np.full(len(ra1), -1)
|
|
sep_as = np.full(len(ra1), 9e9)
|
|
for i in range(len(ra1)):
|
|
dd = np.hypot((ra2 - ra1[i]) * cd, d2 - d1[i]) * 3600.0
|
|
j = np.argmin(dd)
|
|
if dd[j] < tol_as:
|
|
idx[i] = j
|
|
sep_as[i] = dd[j]
|
|
return idx, sep_as
|
|
|
|
|
|
z = np.load(layout.path("_simbad.npz"), allow_pickle=True)
|
|
m = z["otype"] == "GlC"
|
|
gra, gdec = z["ra"][m].astype(float), z["dec"][m].astype(float)
|
|
|
|
hdr = fits.getheader(layout.path("master-Luminance.fit"))
|
|
w = WCS(hdr)
|
|
ny, nx = hdr["NAXIS2"], hdr["NAXIS1"]
|
|
gx, gy = w.all_world2pix(gra, gdec, 0)
|
|
inf = (gx > 20) & (gx < nx - 20) & (gy > 20) & (gy < ny - 20)
|
|
gra, gdec, gx, gy = gra[inf], gdec[inf], gx[inf], gy[inf]
|
|
cx, cy = [float(v) for v in w.all_world2pix(NUC_RA, NUC_DEC, 0)]
|
|
grad = np.hypot(gx - cx, gy - cy) * PIXSCALE / 60.0
|
|
print("known GCs inside the frame: %d (r range %.2f - %.2f arcmin)"
|
|
% (len(gra), grad.min(), grad.max()))
|
|
|
|
d = fits.getdata(layout.path("master-Luminance.fit")).astype(np.float32)
|
|
d = np.ascontiguousarray(d)
|
|
sep.set_extract_pixstack(3000000)
|
|
K = gk()
|
|
|
|
bins = [(0, 2), (2, 4), (4, 8), (8, 14), (14, 30)]
|
|
print("\n%6s %7s | %s" % ("mesh", "Ndet", " ".join("%4.0f-%-4.0f'" % b for b in bins)))
|
|
for bw in (12, 16, 24, 32, 48, 64, 128):
|
|
b = sep.Background(d, bw=bw, bh=bw, fw=3, fh=3)
|
|
ds = d - b.back()
|
|
rm = b.rms()
|
|
o = sep.extract(ds, 3.0, err=rm, minarea=5, filter_kernel=K,
|
|
filter_type="matched", deblend_nthresh=32,
|
|
deblend_cont=0.005, clean=True)
|
|
ora, odec = w.all_pix2world(o["x"], o["y"], 0)
|
|
idx, sp = match(gra, gdec, ora, odec, 2.0)
|
|
row = []
|
|
for lo, hi in bins:
|
|
s = (grad >= lo) & (grad < hi)
|
|
row.append("%3d/%-3d" % ((idx[s] >= 0).sum(), s.sum()))
|
|
print("%6d %7d | %s" % (bw, len(o), " ".join(row)))
|
|
del b, ds, rm, o
|