from pathlib import Path src_dir = str(Path(__file__).parent) + '/src' lc_dir = str(Path(__file__).parent) + '/data/lightcurve' lc_ver = {'DR1': ''} DEFAULT_LC_VERSION = 'DR1' cat_dir = str(Path(__file__).parent) + '/data/catalog' from src.HDF5DataSet import HDF5DataSet from src.CatalogDatabase import CatalogDatabase import numpy as np import warnings from astropy.wcs import WCS, FITSFixedWarning warnings.simplefilter('ignore', category=(FITSFixedWarning, RuntimeWarning)) warnings.filterwarnings( "ignore", message=r"Input line \d+ contained no data and will not be counted towards `max_rows=\d+`" ) def resolve_lc_dir(version=None): """ resolve the light curve data directory based on the version name. """ if version is None: version = DEFAULT_LC_VERSION if version not in lc_ver: warnings.warn(f"Unknown version '{version}', available versions: {list(lc_ver.keys())}. Using default version '{DEFAULT_LC_VERSION}' instead.") version = DEFAULT_LC_VERSION return lc_dir + lc_ver[version] def read_lc(star_full_id, datasets=None, # if not provided, will read from data path defined by star_info catalogdb=None, # if star_info is not provided, need catalog to get star metadata star_info=None, # dict with star metadata, if not provided, will read from catalog bands=['r', 'z'], # bands to read read_cmd=False, read_image=False, version=None, # version name (e.g., 'DR1', 'DR1.5'), defaults to DEFAULT_LC_VERSION verbose=False): # if star_info is None, need to read from catalog if star_info is None: if catalogdb is None: # load catalog database try: catalogdb = CatalogDatabase(cat_dir + '/DREAMS_star_catalog_v1.h5') except Exception as e: raise ValueError("Read catalog database failed.") from e if verbose: print(f"Catalog database loaded from default path: {cat_dir + '/DREAMS_star_catalog_v1.h5'}") # read star_info from catalog try: star_info = catalogdb.search_id(global_ids=star_full_id)[0] except IndexError as e: raise ValueError(f"Star ID {star_full_id} not found in catalog.") from e if verbose: print(f"Star ID {star_full_id} info loaded from catalog database.") # if datasets is None, need to determine datasets based on catalog and load them if datasets is None: if verbose: print("No datasets provided. Loading light curve datasets based on star_info ...") stamp1, stamp2 = star_info['stamp1'], star_info['stamp2'] datasets = {} for stamp in [stamp1, stamp2]: if stamp in ['', '-']: continue datasets[stamp] = {} for band in bands: # result = {'id': star_full_id, 'band': band, 'stamp': stamp} data_dir = resolve_lc_dir(version) datasets[stamp][band] = HDF5DataSet(data_dir + f'/lc_{stamp}_{band}.h5') if verbose: print(f" Loaded dataset for stamp {stamp}, band {band} from {data_dir + f'/lc_{stamp}_{band}.h5'}") if verbose: print("All datasets loaded.") # if datasets is not none, need to make sure it is consistent with star_info and bands else: if verbose: print("Datasets provided. Validating datasets consistency with star_info and bands ...") stamp1, stamp2 = star_info['stamp1'].item(), star_info['stamp2'].item() expected_stamps = [s for s in [stamp1, stamp2] if s not in ['', '-']] if list(datasets.keys()) != expected_stamps: raise ValueError(f"datasets stamps {list(datasets.keys())} != expected {expected_stamps}") for stamp in expected_stamps: if list(datasets[stamp].keys()) != bands: raise ValueError(f"datasets[{stamp}] bands {list(datasets[stamp].keys())} != expected {bands}") if verbose: print("Datasets validation passed.") # read light curves from datasets and organize results light_curves = [] if verbose: print(f"Reading light curves from datasets, star ID {star_full_id} ...") stamp1, stamp2 = star_info['stamp1'], star_info['stamp2'] stamp_to_id = {} if stamp1 not in ['', '-']: stamp_to_id[stamp1] = star_info['internal_id1'] if stamp2 not in ['', '-']: stamp_to_id[stamp2] = star_info['internal_id2'] for stamp in datasets: star_id = stamp_to_id[stamp] for band in datasets[stamp]: dataset = datasets[stamp][band] result = {'id': star_full_id, 'band': band, 'stamp': stamp} result['has_data'] = True # read light curve try: lcdata = dataset.read_lc(star_id) except Exception as e: print(f"Failed to read light curve for star {star_id} in {stamp}_{band}: {e}") result['has_data'] = False light_curves.append(result) continue # read metadata metadata = dataset.read_metadata() mag_zero = metadata['mag_zero'] result['mag_zero'] = mag_zero # read catalog info cat_info = dataset.read_catalog_data() if star_id not in cat_info: print(f"Star ID {star_id} not found in {stamp}_{band}.") result['has_data'] = False light_curves.append(result) continue cat_info = cat_info[star_id] result['ra'] = cat_info['ra'] result['dec'] = cat_info['dec'] result['x'] = cat_info['x'] result['y'] = cat_info['y'] result['cat_mag'] = cat_info['mag'] result['cat_merr'] = cat_info['merr'] # read reference info ref_info = dataset.read_reference_data() ref_wcs = WCS(ref_info['header']) ref_bkg, ref_bkgerr = float(ref_info['header']['PBKGRND']), float(ref_info['header']['PBKGERR']) ref_gain = float(ref_info['header']['GAIN']) ref_nstack = int(ref_info['header']['NSTACK']) result['gain'] = ref_gain f_limit = 3.0*np.sqrt(ref_bkg/ref_gain/ref_nstack) # 3-sigma Poisson flux limit in DN unit if read_image: result['ref_image'] = ref_info['image'] result['ref_mask'] = ref_info['mask'] # read baseline if read_cmd: baseall = dataset.read_baseline() base_flux_all = np.array([b['flux'] for b in baseall]) base_ferr_all = np.array([b['eflux'] for b in baseall]) # if baseline flux is negative, set to the limiting flux base_ferr_all[base_flux_all0: problematic\n') def write_cat_header(handle, light_curve): # cat 文件 header 记录 target 在当前 stamp 中的位置和星等 handle.write( f"# Target star: x={light_curve['x']:.6f} y={light_curve['y']:.6f} ra={light_curve['ra']:.8f} dec={light_curve['dec']:.8f} mag={light_curve['mag_base']:.4f} merr={light_curve['merr_base']:.4f}\n" ) def write_cmd_header(handle, light_curve_r, light_curve_z): # cmd 文件 header 分别记录 target 的 r / z 波段信息 handle.write( f"# Target star r: x={light_curve_r['x']:.6f} y={light_curve_r['y']:.6f} ra={light_curve_r['ra']:.8f} dec={light_curve_r['dec']:.8f} mag={light_curve_r['mag_base']:.4f} merr={light_curve_r['merr_base']:.4f}\n" ) handle.write( f"# Target star z: x={light_curve_z['x']:.6f} y={light_curve_z['y']:.6f} ra={light_curve_z['ra']:.8f} dec={light_curve_z['dec']:.8f} mag={light_curve_z['mag_base']:.4f} merr={light_curve_z['merr_base']:.4f}\n" ) os.makedirs(savedir, exist_ok=True) valid_light_curves = [lc for lc in light_curves if lc['has_data']] if not valid_light_curves: raise ValueError('No valid light curves to save') for light_curve in valid_light_curves: # 每个有效 band 各导出一份 pysis 和 catalog 文件 field = light_curve['stamp'].split('.')[0] band = light_curve['band'] lc_filename = f"lc_{field}_{band}.dat" lc_filepath = os.path.join(savedir, lc_filename) with open(lc_filepath, 'w') as lc_file: write_pysis_header(lc_file, light_curve) lc_file.write('# image_id mjd flux ferr fwhm sky qirr_percent stddev airmass sigma_subt sigma_phot mag merr flag\n') for i in range(light_curve['ndata']): mag = light_curve['mag'][i] merr = light_curve['merr'][i] mjd = light_curve['mjd'][i] lc_file.write( f"{int(light_curve['image_id'][i]):4d} {mjd:.6f} {light_curve['flux'][i]:12.3f} {light_curve['ferr'][i]:10.3f} {light_curve['fwhm'][i]:5.3f} {light_curve['sky'][i]:9.3f} {100*light_curve['qirr'][i]:6.2f} {light_curve['stddev'][i]:6.1f} {light_curve['airmass'][i]:4.2f} {light_curve['sigma_subt'][i]:8.2E} {light_curve['sigma_phot'][i]:8.2E} {mag:7.4f} {merr:6.4f} {light_curve['flag'][i]:1d}\n" ) if all(key in light_curve for key in ['base_x_all', 'base_y_all', 'base_ra_all', 'base_dec_all', 'base_mag_all', 'base_merr_all']): cat_filename = f"catalog_{field}_{band}.dat" cat_filepath = os.path.join(savedir, cat_filename) with open(cat_filepath, 'w') as cat_file: write_cat_header(cat_file, light_curve) cat_file.write('# x y ra dec mag merr\n') for x, y, ra, dec, mag_base, merr_base in zip( light_curve['base_x_all'], light_curve['base_y_all'], light_curve['base_ra_all'], light_curve['base_dec_all'], light_curve['base_mag_all'], light_curve['base_merr_all'] ): cat_file.write(f"{x:8.3f} {y:8.3f} {ra:.8f} {dec:.8f} {mag_base:.4f} {merr_base:.4f}\n") try: # CMD 数据、CMD 图和 finder chart 都使用同一个 selected stamp selected_stamp, light_curve_r, light_curve_z = _select_cmd_light_curves(valid_light_curves, priority_field=priority_field) except ValueError: print(f"Light curve files saved to {savedir}") return field = selected_stamp.split('.')[0] cmd_filename = f"cmd_{field}_r-z.dat" cmd_filepath = os.path.join(savedir, cmd_filename) with open(cmd_filepath, 'w') as cmd_file: write_cmd_header(cmd_file, light_curve_r, light_curve_z) cmd_file.write('# r_x r_y r_ra r_dec r_mag r_merr z_x z_y z_ra z_dec z_mag z_merr\n') nstars = min(len(light_curve_r['base_mag_all']), len(light_curve_z['base_mag_all'])) for i in range(nstars): cmd_file.write( f"{light_curve_r['base_x_all'][i]:8.3f} {light_curve_r['base_y_all'][i]:8.3f} {light_curve_r['base_ra_all'][i]:.8f} {light_curve_r['base_dec_all'][i]:.8f} {light_curve_r['base_mag_all'][i]:.4f} {light_curve_r['base_merr_all'][i]:.4f} " f"{light_curve_z['base_x_all'][i]:8.3f} {light_curve_z['base_y_all'][i]:8.3f} {light_curve_z['base_ra_all'][i]:.8f} {light_curve_z['base_dec_all'][i]:.8f} {light_curve_z['base_mag_all'][i]:.4f} {light_curve_z['base_merr_all'][i]:.4f}\n" ) cmd_plot_path = os.path.join(savedir, f"cmd_{field}_r-z.png") plot_cmd(valid_light_curves, priority_field=priority_field, saveto=cmd_plot_path) finder_chart_path = os.path.join(savedir, 'finder_chart.png') if 'ref_image' in light_curve_z: plot_finder_chart(light_curve_z, saveto=finder_chart_path) print(f"Light curve files saved to {savedir}") if __name__ == '__main__': import argparse parser = argparse.ArgumentParser(description='Read light curves for a given star ID') parser.add_argument('star_id', type=int, nargs='?', default=19178854, help='Global star ID (default: 19178854)') parser.add_argument('--version', type=str, default=None, help='Light curve version (e.g., DR1, DR1.5). Defaults to latest.') args = parser.parse_args() star_full_id = args.star_id lc_version = args.version stamp1, stamp2 = '-', 'D02.N3.0702' catalogdb = CatalogDatabase(cat_dir + '/DREAMS_star_catalog_v1.h5') all_star_info = np.loadtxt(cat_dir + f'/group/stargroup_{stamp1}_{stamp2}.cat', dtype=str) star_info = all_star_info[54] star_info_dict = { 'global_id': int(star_info[0]), 'ra': float(star_info[1]), 'dec': float(star_info[2]), 'mag': float(star_info[3]), 'merr': float(star_info[4]), 'stamp1': star_info[5], 'internal_id1': int(star_info[6]), 'stamp2': star_info[11], 'internal_id2': int(star_info[12]), } datasets = {} data_dir = resolve_lc_dir(lc_version) for stamp in [stamp1, stamp2]: if stamp in ['', '-']: continue datasets[stamp] = {} for band in ['r', 'z']: datasets[stamp][band] = HDF5DataSet(data_dir + f'/lc_{stamp}_{band}.h5') # no other inputs, read light curves internally print(f"{'#'*50:s}\nTest read_lc using only star ID {star_full_id} (version: {lc_version or DEFAULT_LC_VERSION}) ...") light_curves = read_lc(star_full_id, read_cmd=True, read_image=True, version=lc_version, verbose=True) print(f"Successfully read {len(light_curves)} light curves for star ID {star_full_id}.\n{'#'*50:s}\n") # # input id + catalog database only # print(f"{'#'*50:s}\nTest read_lc using star ID {star_full_id} with input catalog database ...") # light_curves = read_lc(star_full_id, catalogdb=catalogdb, read_cmd=True, read_image=True, verbose=True) # print(f"Successfully read {len(light_curves)} light curves for star ID {star_full_id}.\n{'#'*50:s}\n") # # input id + star_info only # print(f"{'#'*50:s}\nTest read_lc using star ID {star_full_id} with input star_info ...") # light_curves = read_lc(star_full_id, star_info=star_info_dict, read_cmd=True, read_image=True, verbose=True) # print(f"Successfully read {len(light_curves)} light curves for star ID {star_full_id}.\n{'#'*50:s}\n") # del all_star_info, star_info # # input id + datasets only # print(f"{'#'*50:s}\nTest read_lc using star ID {star_full_id} with input datasets ...") # light_curves = read_lc(star_full_id, datasets=datasets, read_cmd=True, read_image=True, verbose=True) # print(f"Successfully read {len(light_curves)} light curves for star ID {star_full_id}.\n{'#'*50:s}\n") # del stamp1, stamp2 # # input datasets + star_info # print(f"{'#'*50:s}\nTest read_lc using star ID {star_full_id} with input datasets + star_info ...") # light_curves = read_lc(star_full_id, datasets=datasets, star_info=star_info_dict, read_cmd=True, read_image=True, verbose=True) # print(f"Successfully read {len(light_curves)} light curves for star ID {star_full_id}.\n{'#'*50:s}\n") # # input datasets + catalogdb # print(f"{'#'*50:s}\nTest read_lc using star ID {star_full_id} with input datasets + catalogdb ...") # light_curves = read_lc(star_full_id, datasets=datasets, catalogdb=catalogdb, read_cmd=True, read_image=True, verbose=True) # print(f"Successfully read {len(light_curves)} light curves for star ID {star_full_id}.\n{'#'*50:s}\n") # # input all three: datasets + catalogdb + star_info # print(f"{'#'*50:s}\nTest read_lc using star ID {star_full_id} with input datasets + catalogdb + star_info ...") # light_curves = read_lc(star_full_id, datasets=datasets, catalogdb=catalogdb, star_info=star_info_dict, read_cmd=True, read_image=True, verbose=True) # print(f"Successfully read {len(light_curves)} light curves for star ID {star_full_id}.\n{'#'*50:s}\n") ### 返回的 light_curves 是一个列表，每个元素是一个字典，包含一个波段的光变信息（元信息和光变数据） # 可索引的元信息包括这些，前面标记*的是需要展示的这颗恒星的元信息，只需筛选任意一个 has_data=True 的字典内容展示即可： # *'id': 星的完整编号 # *'band': 波段 # 'stamp': 所处子区域的编号 # 'has_data': 是否有数据文件 # 'mag_zero': 星等零点 # *'ra': 赤经 # *'dec': 赤纬 # *'cat_mag': 星表星等 # *'cat_merr': 星表星等误差 # *'flux_base': 基线流量 # *'ferr_base': 基线流量误差 # *'mag_base': 基线星等 # *'merr_base': 基线星等误差 # *'ndata': 光变曲线测量数据点数量 # 画天区星图相关的信息，参考 plot_finder_chart 示例： # 'ref_image': 仅当read_image=True时返回, 参考图像, numpy二维数组(600x600) # 'ref_mask': 仅当read_image=True时返回, 参考图像掩膜, numpy二维数组(600x600) # 'x': 在参考图像中的X像素坐标 # 'y': 在参考图像中的Y像素坐标 # 'base_mag_all': 仅当时read_cmd=True时返回, 图像中所有恒星的基线星等, numpy一维数组(ndata,) # 'base_merr_all': 仅当时read_cmd=True时返回, 图像中所有恒星的基线星等误差, numpy一维数组(ndata,) # 可索引的光变信息包括这些（每个都是长度为 ndata 的 numpy 数组），前面标记*的是需要展示的信息： # *'mjd': 观测时间（简化儒略日 Modified Julian Day） # *'flux': 流量 # *'ferr': 流量误差 # *'mag': 星等 # *'merr': 星等误差 # 'fwhm': 图像半峰全宽 # 'sky': 天空背景 # 'qirr': 像质不对称度 # 'airmass': 大气质量 # 'stddev': 原始图像标准差 # 'sigma_subt': 图像相减残差标准差 # 'sigma_phot': 测光残差标准差 # *'flag': 数据质量标记，0为良好，1为图像相减异常，2为测光异常 # for l in light_curves: # if l['has_data']: # plot_finder_chart(l) # # print(l.keys()) # plot_cmd(light_curves) save_lc_files(light_curves, savedir=str(Path(__file__).parent) + f'/output/{star_full_id:08d}')