from astropy.io import fits
import numpy as np
from astropy.nddata import Cutout2D
from astropy.wcs import WCS
from astropy.coordinates import SkyCoord
from astropy import units as u
from astropy.table import Table
import glob
import os

#filenames = ['../SUPM139A914F02F38600.fits', '../SUPM139B30920E507F00.fits', '../SUPM139B3B860362D800.fits', '../SUPM139B6C0D02B1F800.fits', '../SUPM139BA4B50383E300.fits', '../SUPM139C01B30E611200.fits']
filenames = []
filenames.extend(glob.glob('PSG/G*.fits'))
filenames.extend(glob.glob('PSG/Mega*.fits'))
filenames.extend(glob.glob('u_adami/*mean.fits'))
filenames.extend(glob.glob('WIRCAM/*.fits'))
print filenames
psg_list = Table.read('/u/kirg/PSG_list.fits')
psg_list.add_row(['587741602572861450', 194.9839, 27.7464])
psg_list.add_row(['3016', 195.0043, 28.0821])

print psg_list
for filename in filenames:
	frame = fits.open(filename)[0]
	wcs_inp = WCS(frame.header)
	base = os.path.basename(filename)
	base_name = os.path.splitext(base)[0]
	for obj in psg_list:
		inp_coords = "%.5fd %.5fd" % (obj['RA'], obj['DEC'])
		coords = SkyCoord(inp_coords)
		pixcrd = wcs_inp.wcs_world2pix(obj['RA'], obj['DEC'], 0)
		try:
			if ((pixcrd[0] > 0) & (pixcrd[1] > 0) & (pixcrd[0] < frame.data.shape[1]) & (pixcrd[1] < frame.data.shape[0])): 
				print pixcrd, obj['NAME']
				size = u.Quantity(300, u.arcsec)
				cutout = Cutout2D(frame.data,coords, size, wcs=WCS(frame.header))
				hdul=fits.HDUList([fits.PrimaryHDU(data=cutout.data, header=cutout.wcs.to_header())])
				hdul.writeto('PSG_cutouts/'+str(obj['NAME'])+'_'+frame.header['FILTER']+'_'+base_name+'_extended.fits', overwrite=True)
		except ValueError:
			print filename
		except AttributeError:
                        print filename