from astropy.io import fits
import sys
from astropy.wcs import WCS
from astropy import wcs
from astropy.table import QTable, Table, Column

from astropy.coordinates import SkyCoord # High-level coordinates
from astropy.coordinates import ICRS, Galactic, FK4, FK5 # Low-level frames
from astropy.coordinates import Angle, Latitude, Longitude # Angles
import astropy.units as u





filename='Gaia_unWISE_Coma.fits.catalog'
fout=filename.replace(".fits.catalog", ".footprint.reg")


hdul = fits.open(filename)
#hdul.info()

tbdata = hdul[1].data
hdul.close()

detmask='../../products/mosa_parts_tm0_DetectorMask_en00.fits'
print("reading sensmap {}".format(detmask))
hdul = fits.open(detmask)
mask = hdul[0].data
hdr = hdul[0].header
hdul.close()
wcs = WCS(hdr)


tab = Table.read(filename, format='fits')
tbr=[]
for idx in range(len(tab['RA'])):
    ra=tab['RA'][idx]
    dec=tab['DEC'][idx]
    crd = SkyCoord(ra, dec, frame=FK5(), unit="deg")
    pix = wcs.wcs_world2pix([(ra, dec),], 1)
    px=round(pix[0][0])-1
    py=round(pix[0][1])-1
    if not (px >=0 and py >= 0 and px <= (10000-1) and py < (10000-1)):
        tbr.append(idx)
        continue
    if mask[py,px] == 0:
        tbr.append(idx)
tab.remove_rows(tbr)
tab.write("Gaia_unWISE_Coma.footprint.fits.catalog", format='fits', overwrite='True')



with open("./{}".format(fout), 'w') as writer:
    for rec in tbdata:
        ra=rec['RA']
        dec=rec['DEC']
        crd = SkyCoord(ra, dec, frame=FK5(), unit="deg")
        pix = wcs.wcs_world2pix([(ra, dec),], 1)
        px=round(pix[0][0])-1
        py=round(pix[0][1])-1
        if not (px >=0 and py >= 0 and px <= (10000-1) and py < (10000-1)):
            continue
        if mask[py,px] == 1:
            print("fk5;circle({}, {}, 0.008)".format(rec['RA'],rec['DEC']))
            writer.write("fk5;circle({}, {}, {})\n".format(rec['RA'],rec['DEC'],0.0080000))