""" Synopsis Magnetogram Stitcher Usage: py -3 make_combined.py 2026-02-25T00:00:00 py -3 make_combined.py 2024-02-25T00:00:00 out.fits py -3 make_combined.py 2026-02-25T00:00:00 out.png py -3 make_combined.py # interactive """ import os, sys, re import numpy as np from astropy.io import fits from astropy.time import Time import warnings warnings.filterwarnings("ignore", category=Warning) FITS_DIR = os.path.dirname(__file__) or "." DAT_FILE = os.path.join(FITS_DIR, "carr_long.dat") OUT_DIR = FITS_DIR # ---- 1. 加载 carr_long.dat（必需，不扫描 FITS） ---- def load_cr_table(): rows = [] for line in open(DAT_FILE): p = line.strip().split() if len(p) < 2: continue try: cr = int(p[0]) t = Time(p[1].replace("/", "-"), format="iso", scale="utc") rows.append((cr, t)) except Exception: continue rows.sort(key=lambda x: x[0]) return rows # ---- 2. 二分查找：目标时间落在哪个 CR ---- def find_cr(target_t, table): n = len(table) lo, hi = 0, n while lo < hi: mid = (lo + hi) // 2 if target_t < table[mid][1]: hi = mid else: lo = mid + 1 idx = lo - 1 if idx < 0: raise ValueError(f"Date earlier than CR{table[0][0]} in {DAT_FILE}") if idx + 1 >= n: raise ValueError(f"Date within last CR in {DAT_FILE}: CR{table[idx][0]}, no next CR") return table[idx][0] # ---- 3. FITS 路径（文件名直接构造，不扫描） ---- def fits_path(cr): return os.path.join(FITS_DIR, f"fmg.Synoptic_Ml.{cr}.fits.gz") # ---- 4. 拼图 ---- def stitch(target_str, out_path=None): # 4.1 解析时间 t = Time(target_str.replace("_", "-").replace(" ", "T")) # 4.2 carr_long.dat → 找 CR table = load_cr_table() cr1 = find_cr(t, table) # 找下一个 CR（table 中紧邻的） idx1 = next(i for i, (c, _) in enumerate(table) if c == cr1) cr2 = table[idx1 + 1][0] # 4.3 文件路径（直接构造） f1 = fits_path(cr1) f2 = fits_path(cr2) if not os.path.isfile(f1): raise FileNotFoundError(f"Missing: {os.path.basename(f1)}") if not os.path.isfile(f2): raise FileNotFoundError(f"Missing: {os.path.basename(f2)}") # 4.4 读数据 with fits.open(f1, memmap=False) as hdul: d1 = np.array(hdul[0].data, dtype=np.float64) h1 = dict(hdul[0].header) with fits.open(f2, memmap=False) as hdul: d2 = np.array(hdul[0].data, dtype=np.float64) nx = int(h1["NAXIS1"]) ny = int(h1["NAXIS2"]) # 4.5 col_cut + B0 重算（基于 CR1 FITS 头秒级时间） t_start_hdr = Time(str(h1["T_START"]).strip()) t_stop_hdr = Time(str(h1["T_STOP"]).strip()) total = (t_stop_hdr - t_start_hdr).to_value("day") elapsed = (t - t_start_hdr).to_value("day") frac = max(0.0, min(1.0, elapsed / total)) col_cut = int(round(frac * (nx - 1))) b0_a = float(h1["B0_FRST"]) b0_b = float(h1["B0_LAST"]) b0_slope = (b0_b - b0_a) / total b0_t = b0_a + frac * (b0_b - b0_a) cr_period = t_stop_hdr - t_start_hdr t_stop_out = t + cr_period b0_frst_out = b0_t - b0_slope * elapsed b0_last_out = b0_t + b0_slope * cr_period.to_value("day") # 4.6 横向拼接 # Col: 0 ---- col_cut ---------- nx ----- 2*nx-1 # [丢弃][ CR1后半+CR2前半 ][丢弃] combined = np.concatenate([d1, d2], axis=1) out_data = combined[:, col_cut: col_cut + nx].copy() # # 平滑过渡 ±50 列 # bw = 50 # L = max(0, col_cut - bw) # R = min(nx, col_cut + bw) # if R > L: # w = np.linspace(0.0, 1.0, R - L) # out_data[:, L:R] = out_data[:, L:R] * (1.0 - w) + d2[:, L:R] * w # 4.7 输出路径 if out_path is None: # 格式: fmg.Synoptic_Ml.CR1_CR2.YYYYMMDD_HHMMSS.fits dt_part = t.datetime.strftime("%Y%m%d_%H%M%S") out_path = os.path.join(OUT_DIR, f"fmg.Synoptic_Ml.{cr1}_{cr2}.{dt_part}.fits") ext = os.path.splitext(out_path)[1].lower() # 4.8 写文件 if ext == ".fits": hdu = fits.PrimaryHDU(out_data) hdu.header["SIMPLE"] = True hdu.header["BITPIX"] = -64 hdu.header["NAXIS"] = 2 hdu.header["NAXIS1"] = nx hdu.header["NAXIS2"] = ny hdu.header["BSCALE"] = float(h1.get("BSCALE", 1.0)) hdu.header["BZERO"] = float(h1.get("BZERO", 0.0)) hdu.header["CTYPE1"] = str(h1.get("CTYPE1", "Carrington longitude")) hdu.header["CTYPE2"] = str(h1.get("CTYPE2", "Latitude")) hdu.header["CRPIX1"] = float(h1.get("CRPIX1", nx / 2.0)) hdu.header["CRPIX2"] = float(h1.get("CRPIX2", ny / 2.0)) hdu.header["CRVAL1"] = float(h1.get("CRVAL1", 0.0)) hdu.header["CRVAL2"] = float(h1.get("CRVAL2", 0.0)) hdu.header["CDELT1"] = float(h1.get("CDELT1", 1.0)) hdu.header["CDELT2"] = float(h1.get("CDELT2", 1.0)) hdu.header["BUNIT"] = str(h1.get("BUNIT", "Gauss")) hdu.header["R_SUN"] = float(h1.get("R_SUN", 0.0)) hdu.header["TELESCOP"] = str(h1.get("TELESCOP", "ASO-S/FMG")) hdu.header["ORIGIN"] = "SODC" hdu.header["CAR_ROT"] = f"{cr1}+{cr2}" hdu.header["CR1"] = int(cr1) hdu.header["CR2"] = int(cr2) hdu.header["T_START"] = t.iso hdu.header["T_STOP"] = t_stop_out.iso hdu.header["B0_FRST"] = float(b0_frst_out) hdu.header["B0_LAST"] = float(b0_last_out) hdu.header["B0_T"] = float(b0_t) hdu.header["COL_CUT"] = col_cut hdu.header["FRAC"] = float(frac) hdu.writeto(out_path, overwrite=True, output_verify="fix") elif ext in (".png", ".pdf"): import matplotlib.pyplot as plt fig, ax = plt.subplots(figsize=(18, 9), dpi=120) vmax = np.nanpercentile(np.abs(out_data), 99.5) ax.imshow(out_data, cmap="RdBu_r", vmin=-vmax, vmax=vmax, origin="lower", extent=[0, 360, -90, 90], aspect="auto") plt.colorbar(ax.images[0], ax=ax, label="Br [Gauss]") ax.set_title( f"360-deg Synoptic Magnetogram | {t.iso[:10]} | " f"CR{cr1}+CR{cr2} | B0={b0_t:.2f} deg | col_cut={col_cut}", fontsize=11) ax.set_xlabel("Carrington Longitude [deg]") ax.set_ylabel("Latitude [deg]") plt.tight_layout() plt.savefig(out_path, dpi=150) plt.close() elif ext == ".npy": np.save(out_path, out_data) # 4.9 摘要 total_sec = cr_period.to_value("day") * 86400 print(f"\n{'='*60}") print(f" T_START : {t.iso}") print(f" T_STOP : {t_stop_out.iso}") print(f" 周期 : {cr_period.to_value('day'):.6f} d ({total_sec:.0f} s)") print(f" 每度 : {total_sec/360:.2f} s") print(f" CR : CR{cr1} + CR{cr2}") print(f" 文件 : {os.path.basename(f1)} + {os.path.basename(f2)}") print(f" col_cut : {col_cut} (frac={frac:.5f})") print(f" B0_T : {b0_t:.4f} deg") print(f" B0_FRST : {b0_frst_out:.4f} deg (T_START)") print(f" B0_LAST : {b0_last_out:.4f} deg (T_STOP)") print(f" shape : {out_data.shape}") print(f" 数值范围 : [{np.nanmin(out_data):.1f}, {np.nanmax(out_data):.1f}] Gauss") print(f" 输出 : {os.path.basename(out_path)}") print(f"{'='*60}\n") return out_data # ---- main ---- def _fmt(t_str): """标准化用户输入时间字符串：分隔符统一为 T""" return t_str.replace(" ", "T").replace("-_-", "T") if __name__ == "__main__": if len(sys.argv) == 1: inp = input("起始时间 (默认 2024-02-01T00:00:00): ").strip() target = inp or "2024-02-01T00:00:00" out = None elif len(sys.argv) == 2: target = _fmt(sys.argv[1]) out = None else: target = _fmt(sys.argv[1]) out = sys.argv[2] stitch(target, out)