object_extract_rate.py

import os

import numpy as np
import matplotlib.pyplot as plt
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.io import fits, votable
from astropy.wcs import WCS
from astropy.table import Table
from astropy.coordinates import Angle
import pandas as pd
'''
midify time:2021-08-10
verion:V2.0
Note: modify the input file path and name(137rows)
'''
def qurey (acatfile2,wcsfitsfile,refc,newxypanstar,regionpanstar):
    print('----------------------------------------')
    print('obs file:  ',acatfile2)
    # 提取观测数据*3acat文件中恒星数据
    hdu=  fits.open(acatfile2)
    obsdata = hdu[1].data
    obs_ra = obsdata['RA']
    num_obs_star = len(obsdata)
    obs_dec = obsdata['DEC']
    obs_mag = obsdata['VIS_MAG']
    print('obs star:  ', len(obs_ra))
    obsc = SkyCoord(obs_ra, obs_dec, unit=(u.deg, u.deg), frame='fk5')


    # 读取星表数据（ra,dec,mag）
    f = fits.open(wcsfitsfile) #读取wcs信息，将re dec 转换为x y,并写入newxypanstar***.csv中，
    w = WCS(f[0].header)    # 每张图像的WCS不相同，因此有几张图，就生成几个带有x,y坐标的星表
    refdata = Table.read(refc)
    refdata = np.array(refdata)
    print('ref star:  ', len(refdata))
    ref_ra0 = refdata['raMean']
    ref_dec0 = refdata['decMean']
    ref_gmag0 = refdata['gMeanPSFMag']
    ref_rmag0 = refdata['rMeanPSFMag']
    ref_vmag0 = refdata['vMeanPSFMag']
    ref_coord = list(zip(ref_ra0, ref_dec0))
    ref_coord = SkyCoord(ref_coord, frame='fk5', unit='degree')
    c = w.world_to_pixel(ref_coord)
    x0 = c[0]
    y0 = c[1]
    df = pd.DataFrame(
        {'raMean': ref_ra0, 'decMean': ref_dec0, 'gMeanPSFMag': ref_gmag0, 'rMeanPSFMag': ref_rmag0,
         'vMeanPSFMag': ref_vmag0, 'x': x0, 'y': y0})
    df.to_csv(newxypanstar, index=False)

    '''
    读取带有xy坐标的星表，根据xy坐标提取与观测图像区域一致的参考恒星，并将其写入4681***.region.csv中
    '''

    refdata1 = Table.read(newxypanstar)
    x1 = refdata1['x']
    y1 = refdata1['y']
    mask = (x1 < 4076) & (x1 > 20) & (y1 < 4076) & (y1 > 20)
    newrefdata = refdata1[mask]
    print('region ref star',len(newrefdata))

    # 将参考数据的检索范围内的恒星数据提取到*4acat文件中
    ref_ra1 = newrefdata['raMean']
    ref_dec1 = newrefdata['decMean']
    ref_gmag1 = newrefdata['gMeanPSFMag']
    ref_rmag1 = newrefdata['rMeanPSFMag']
    ref_vmag1 = newrefdata['vMeanPSFMag']
    x1 = newrefdata['x']
    y1 = newrefdata['y']
    # print(r1,d1)
    df = pd.DataFrame(
        {'raMean': ref_ra1, 'decMean': ref_dec1, 'gMeanPSFMag': ref_gmag1, 'rMeanPSFMag': ref_rmag1,
         'vMeanPSFMag': ref_vmag1, 'x': x1, 'y': y1})
    # print(df)
    df.to_csv(regionpanstar, index=False)


    '''
    提取公共区域内大于极限星等的参考恒星，与观测数据进行交叉匹配，匹配成功的数量作为提取概率的分子，
    公共区域内大于极限星等的参考恒星数量作为提取概率的分母
    '''
    ls_mag = []
    ls_extractor_rate = []
    ls_false_rate = []
    for mag in range(10, 20, 1):
        ls_mag.append(mag)
        mask3 = (ref_vmag1 < mag)  # 提取星表中大于极限星等的恒星
        # print(type(newrefdata))
        new_refdata = newrefdata[mask3]
        # print(new_refdata)
        new_refdata = np.array(new_refdata)
        denominator = len(new_refdata)
        print('region ref star(refmag<limitmag)fenmu:  ', denominator)

        mask4 =(obs_mag<mag)
        new_obsdata = obsdata[mask4]
        new_obsdata = np.array(new_obsdata)
        num_new_obsstar = len(new_obsdata)  #xiao yu jixian xingdeng de guance hengxing shuliang


        ref_ra = new_refdata['raMean']
        ref_dec = new_refdata['decMean']
        ref_radec = list(zip(ref_ra, ref_dec))
        newrefc = SkyCoord(ref_radec, frame='fk5', unit='degree')


        odd = open('extract_rate.txt', 'a')
        idx1, d2d1, d3d1 = newrefc.match_to_catalog_sky(obsc, nthneighbor=1)
        idindex = []
        for j in range(len(idx1)):
            idindex.append(j)

        uidlim1 = d2d1[idindex].arcsecond < 3
        # print(uidlim1)
        t = []
        for i in range(len(uidlim1)):
            if uidlim1[i] == True:
                t.append(uidlim1[i])
            else:
                i = i + 1
        fenzi = len(t)
        print('region match star(fenzi)', len(t))

        extractor_odds =fenzi/ denominator  # denominator 所在检索区域内星等小于极限星等的参考恒星的数量
        false_odds = abs(num_new_obsstar - fenzi) / denominator
        print(str('%.2f' % extractor_odds) + ' ' + str('%.2f' % false_odds))
        ls_extractor_rate.append(extractor_odds)
        ls_false_rate.append(false_odds)

        odd.write(acatfile2[num:-5] + ' ' + str('%.2f' % extractor_odds) + '(Extraction rate)' + '  ' + str(
            '%.2f' % false_odds) + '(False rate)' +'  limitmag'+str(mag)+ '\n')


    plt.scatter(ls_mag,ls_extractor_rate)
    # plt.scatter(ls_mag,ls_false_rate)
    plt.grid(linewidth = 0.5)
    plt.plot(ls_mag,ls_extractor_rate)
    # plt.plot(ls_mag,ls_false_rate)
    # 设置图表标题，并给坐标轴添加标签
    # plt.title("limit magnitude", fontsize=20)
    plt.ylim(0,1)
    plt.xlabel("Vmag ", fontsize=12)
    plt.ylabel("Probability", fontsize=12)

    # 设置数字标签
    for a, b in zip(ls_mag, ls_extractor_rate):
        plt.text(a, b, '%.2f'%b, ha='center', va='bottom')
    # for c, d in zip(ls_mag, ls_false_rate):
    #     plt.text(c, d, '%.2f'%d, ha='center', va='bottom')
    plt.legend(['Extraction probability','False probability'],loc='lower left')
    # 设置坐标轴刻度标记的大小
    # plt.tick_params(axis='both', labelsize=10)
    # plt.show()
    extractor_png ='extractor_odd'+ acatfile2[0:-5]+'png'
    plt.savefig(extractor_png)
    # plt.show()
    plt.close()
    print('please check result---extractor_odd*png')

cmd = 'ls *2acat >acatfile2.dat'
os.system(cmd)
alist = [line.rstrip()for line in open('acatfile2.dat')]
num = 0
for f in alist:
    acatfile2 = f
    refc = 'panstarr_ali.csv'
    wcsfitsfile = f[0:-6]+'.fits'
    newxypanstar = 'newxypanstar'+acatfile2[num:-6]+'.csv'
    regionpanstar = 'region_panstar'+acatfile2[num:-6] +'.csv'
    qurey(acatfile2,wcsfitsfile,refc,newxypanstar,regionpanstar)

cmd2 = 'rm newxypanstar* region_panstar*'
os.system(cmd2)