#!/usr/bin/python
# Usage example:
#     python caom_cutout.py cE_candidates.csv -b 30 > result_caom_cE_candidates.csv
# This returns 30" cutouts around positions from cE_candidates.csv, which is a file with "Name RA Dec" lines, possibly comma-separated.
# Then you can load result.csv in TOPCAT and setup activation action to open "cutout_url" column e.g. in Aladin or ds9
#
# Query for David Valls-Gabaud on nearby early type galaxies observed in UV with HST:
#     python caom_cutout.py nearby_ellipticals2.csv -i WFPC2 ACS WFC3 NICMOS -f F225W F330W > result_caom_nearby_ellipticals_hst_f225w.csv
import sys
import urllib
import argparse
URL = "http://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/tap/sync"
params={'REQUEST': 'doQuery',
        'LANG':    'ADQL',
        'FORMAT':  'csv'}
cutoutParams={'URI': '',
              'CUTOUT' : ''}
dataURL = "http://www.cadc.hia.nrc.gc.ca/caom2ops/cutout"
# Beware! a.productType is empty for the HST data!
baseQuery = """
            SELECT 
                '%s' as Target,
                o.collection,
                p.energy_bandpassName,
                p.time_exposure,
                a.uri as accessURL,
                o.instrument_name, 
                a.productType, 
                p.dataProductType,
                p.dataRelease
            FROM 
                caom2.Observation AS o 
                JOIN caom2.Plane AS p 
                    ON o.obsID=p.obsID 
                JOIN caom2.Artifact AS a 
                    ON p.planeID=a.planeID
            WHERE 
                CONTAINS(POINT('ICRS',%f,%f),position_bounds) = 1 
                AND p.calibrationLevel >= 2
                --AND a.productType = 'science'
                AND ((a.productType not in ('auxiliary','preview')) OR (a.productType IS NULL))
	        AND %s
	        AND %s
	        """
cutoutRegion = "circle ICRS %f %f %f"
default_box = 30. # arcsec
default_instr = ['WFPC2', 'ACS', 'WFC3', 'NICMOS', 'MegaPrime', 'WIRCam']
# command line arguments parsing
parser = argparse.ArgumentParser(description='CAOM query script (uses ObsTAP to query CADC archive)')
group = parser.add_mutually_exclusive_group()
parser.add_argument("input_file", help="input ascii file name with 3 columns required (name, ra, dec) either whitespace- or comma-separated")
parser.add_argument("-b", "--box", help="cutout box size in arcsec (default: %.1f)" % default_box, type=float, default=default_box)
group.add_argument("-i", "--instruments", help="instruments to query (default: %s)" % (" ".join(default_instr)), nargs="+", type=str, default=default_instr)
group.add_argument("-ni", "--no-instruments", help="do not filter on instruments, get everything", action="store_true")
parser.add_argument("-f", "--filters", help="telescope filters to check, e.g. g.MP9401, F225W", nargs="+", type=str)
parser.add_argument("-v", "--verbose", help="go verbose", action="store_true")
args = parser.parse_args()
positions = file(args.input_file).readlines()
cutout = float(args.box) / 3600.0 / 2.0
print "object, collection, bandpass, exptime, instrument_name, product_type, data_type, cutout_url"
for position in positions:
    v = position.split()
    if len(v) != 3: # CSV input instead of whitespace-separated
        v = position.split(',')
    name = v[0].replace(',', '').replace(' ', '')
    try:
        ra = float(v[1].replace(',', '').replace(' ', ''))
        dec = float(v[2].replace(',', '').replace(' ', ''))
        if args.verbose:
            print >> sys.stderr, 'Querying %f %f' % (ra, dec)
    except ValueError:
        if args.verbose:
            print >> sys.stderr, 'Skipping header line'
        continue
    if not args.no_instruments:
        instruments = "o.instrument_name in ('%s')" % ("','".join(args.instruments))
    else:
        instruments = "1=1"
    if args.filters:
        filters = "p.energy_bandpassName IN ('%s')" % ("','".join(args.filters))
    else:
        filters = "1=1"
    query = baseQuery % ( name, ra, dec, instruments, filters )
    params['QUERY'] = query
    if args.verbose:
        print >> sys.stderr, query
        print URL, urllib.urlencode(params)
    result = urllib.urlopen(URL, urllib.urlencode(params))
    title = result.readline()
    for res_cur in result:
        res_cur = res_cur.replace('\n', '')
        output_params = res_cur.split(',')
        cutoutParams['URI'] = output_params[4]
        cutoutParams['CUTOUT'] = cutoutRegion % (ra, dec, cutout)
        cutoutURL = dataURL + "?" + urllib.urlencode(cutoutParams)
        print "%s, %s, %s, %s, %s, %s, %s, %s" % (output_params[0], output_params[1], output_params[2], output_params[3],
                                                    output_params[5], output_params[6], output_params[7], cutoutURL)
#        print output_params[0]+"_"+output_params[2]+".fits.gz"+" \""+cutoutURL+"\""