##^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009
#! /usr/bin/env python

##This script will use the TABLES routine TPROJECT to dump the RA and DEC
##from a catalog entry into a table. This table is used by the by the WCSTRAN
##routine to convert the RA/DEC coordinates to X and Y coordinates in the image
##[FitsFileName].
##
##The ID, Magnitude, RA, DEC, etc. for the sources that have X and Y that 
##fall in the image are then added to the original catalog and these entries
##are put in a format that can be used by the MKCATALOG routine in 
##digiphot/photcal.
##
##The user must specify the columns in the input catalog that contain the 
##RA and DEC and the format of the RA (hours or degrees) in the program.
##
#CALLING SEQUENCE:
#run ~/python/pyraf/SDSSCatalogToTable.py '/nfs/slac/g/ki/ki09/lances/Reduced/07Apr26/M13/Photometry/M13SDSSPhotometricTargets.csv' '/nfs/slac/g/ki/ki09/lances/Reduced/07Apr26/M13/M13_G_SlopeFlatFieldedCenteredDitheredMeanWCS.fits' 

import sys,os,re,string,getopt
from pyraf import iraf
import pyfits
import numpy
from numpy import *
from pylab import *
import TableIO

#Now parse the command line keyword arguments 
keywords = ['TvFlag=']
TvFlag      = 0

opts, extraparams = getopt.getopt(sys.argv[3:],'',keywords)
for o,p in opts:
  if o in ['--TvFlag']:
    TvFlag=int(p)

#The filename will be specified as the first argument. These files have
#the form [Object Name]_[Filter Used]_SlopeFlatFielded...fits
TableName    = sys.argv[1]
FitsFileName = sys.argv[2]
FitsHDU      = pyfits.open(FitsFileName)
FitsHeader   = FitsHDU[0].header
Naxis1	     = FitsHeader['NAXIS1']
Naxis2       = FitsHeader['NAXIS2']

#Extract various keys from the full path
BaseName     = os.path.basename(FitsFileName)
DirName      = os.path.dirname(FitsFileName)
BaseParts    = BaseName.split('_')
ObjectName   = BaseParts[0]
ThisFilter   = BaseParts[1]
BaseParts    = BaseName.split('.')
BaseKey      = BaseParts[0]

BaseName     = os.path.basename(TableName)
DirName      = os.path.dirname(TableName)
BaseParts    = BaseName.split('.')
BaseKey      = BaseParts[0]

#Create a directory for Astrometry output
OutDir = DirName+'/'
if not os.path.exists(OutDir): os.mkdir(OutDir)

#The table to which the output will be written
OutWCSTableName = OutDir+BaseKey+'RADEC.tab'
OutPhyTableName = OutDir+BaseKey+'_'+ThisFilter+'_'+'XY.tab'
OutStandCatalog = OutDir+BaseKey+'Catalog.dat'

#Import the good packages from IRAF
iraf.digiphot(_doprint=0)
iraf.daophot(_doprint=0)
iraf.immatchx(_doprint=0)

#print out the description of the catalog
iraf.tinfo.setParam('table', TableName)
iraf.tinfo()

#Now extract the columns that contain the RA and DEC
#First, write it to file and then to an array
iraf.tproject.setParam('intable',TableName)
iraf.tproject.setParam('outtable',OutWCSTableName)
iraf.tproject.setParam('columns','c8,c9')
iraf.tproject(mode='h')
RADEC=iraf.tproject(Stdout=1, mode ='h')

#Now get all of the columns to do a full-fledge cataloging
iraf.tproject.setParam('columns','c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12')
AllCols = iraf.tproject(Stdout=1, mode ='h')
#Then to convert the RA/DEC to x,y coordinates, we could use the command line
#
#$> sky2xy M13_G_SlopeFlatFieldedWCSTest.fits @M13UBVPhotometryAndProperMotionsKadla1966RADEC.tab > Test.txt
#
#But this outputs the coordinates with a lot of funky extra text. So use 
#iraf.imcoords.wcsctran
iraf.wcsctran.setParam('input',OutWCSTableName)
iraf.wcsctran.setParam('output',OutPhyTableName)
iraf.wcsctran.setParam('image',FitsFileName)
iraf.wcsctran.setParam('inwcs','world')
iraf.wcsctran.setParam('outwcs','tv')
iraf.wcsctran.setParam('formats','%10.4f %10.4f')
iraf.wcsctran.setParam('columns','1,2')
iraf.wcsctran(mode='h')
#Set up vectors to hold X and Y coordinates that fall in image boundaries
iraf.wcsctran.setParam('output',"STDOUT")
XY           = iraf.wcsctran(Stdout=1,mode='h')
XYLen        = len(XY)
XCoords      = []
YCoords      = []
XYCoordList  = []
MatchNum = 0

#Open up a file to write these values to 
F=open(OutStandCatalog,'w')

for i, line in zip(range(XYLen), XY):
    if re.search('\d\.\d\d\d\d',line): 
      XYRow = line.split()
      X     = float(XYRow[0])
      Y     = float(XYRow[1])
      print 'X :'+str(X) + ' Y: '+str(Y)
      if X> 0 and X < Naxis1 and Y > 0 and Y < Naxis2 : 
        print 'Found star in image at X:' + str(X) + ' Y: ' +str(Y)
        XCoords = append(XCoords, X)
        YCoords = append(YCoords, Y)
        MatchNum+=1
        F.write(str(X).ljust(10))
        F.write(str(Y).ljust(10)+'\n')
	XYCoordList.extend([str(X)+'   '+str(Y)])

#Close the file
F.close()

#Write out the X and Y coordinates in a more usable form
XYList = [XCoords,YCoords]
XYArr  = numpy.array(XYList).transpose()
TableIO.writeArray(OutPhyTableName, XYArr)

#Now display them on an image if so desired
if TvFlag > 0:
  Im = pyfits.getdata(FitsFileName)
  MaxDisp = Im.mean()+3*Im.std()
  iraf.set(stdimage='imt4096')
  iraf.display(FitsFileName,1,z1=1,z2=MaxDisp,zrange='no')
  iraf.tvmark.setParam('coords',"STDIN")
  iraf.tvmark(1,Stdin=XYCoordList,col=205)