""" Code to query catalogues
Context : SRP.GW
Module  : SRPGWQuery
Author  : Stefano Covino
Date    : 27/03/2019
E-mail  : stefano.covino@brera.inaf.it
URL:    : http://www.merate.mi.astro.it/utenti/covino
Purpose : Query catalogues

usage: SRPGWQuery [-h] [-B] [-c RA DEC] [-d] [-g gaiarad] -i inputfile -o
                  outfile [-s simbadrad] [-v] [-w wsubpath] [--version]

optional arguments:
  -h, --help            show this help message and exit
  -B, --BoMySQL         IGAIA query to MySQL server in Bologna
  -c RA DEC, --coordcols RA DEC
                        Coordinate column labels
  -d, --debug           Debug information
  -g gaiarad, --gaia gaiarad
                        GAIA catalogue query radius (arcsec)
  -i inputfile, --inputfile inputfile
                        Input table
  -o outfile, --outfile outfile
                        Output table
  -s simbadrad, --simbad simbadrad
                        Simbad catalogue query radius (arcsec)
  -v, --verbose         Fully describe operations
  -w wsubpath, --weight wsubpath
                        Get weight frame value (substitute path)
  --version             show program's version number and exit


History : (12/11/2015) First version.
		: (03/01/2016) Work for single tables too.
        : (16/01/2016) Better management of multiple entries.
        : (21/01/2016) Area weights added.
        : (22/01/2016) Get pixel
        : (02/02/2016) Better path management.
        : (03/02/2016) FWHM from FITS headers.
        : (05/02/2016) Neighbors identification.
        : (06/02/2016) Separation value.
        : (10/02/2016) More rationale reading from FITS header.
        : (11/02/2016) Quicker access to FITS data.
        : (15/02/2015) Minor planet search.
        : (17/02/2016) xMatch IGAIA and Simbad search.
        : (04/03/2016) Delta mag in neighbor.
        : (06/03/2016) Score and quicker minor planet search.
        : (16/05/2016) Queries for more telescopes.
        : (31/05/2016) LEDA catalogue.
        : (07/06/2016) Improved score.
        : (30/11/2016) GAIA DR1 catalogue.
        : (07/12/2016) More flexible input coordinate format, better 
            management of non standard tables and test statistics.
        : (15/12/2016) Meta data for FWHM
        : (18/12/2016) FWHM computation.
        : (24/12/2016) Chi2 in score.
        : (05/01/2017) PanSTARRS catalogue.
        : (25/02/2017) USNO catalogue.
        : (02/03/2017) Update.
        : (16/05/2017) Minor update.
        : (05/03/2019) GLADE galaxy catalogue.
        : (06/03/2019) More flexible date reading.
        : (27/03/2019) MINMAG in score.
"""

__version__ = '2.17.0'


import argparse, os
import numpy
from astropy.table import Column, Table
import pyprind
import SRPGW as GW
from SRPGW.GetSimbad import GetSimbad
from SRPGW.GetSimbadxMatch import GetSimbadxMatch
from SRPGW.GetGAIAxMatch import GetGAIAxMatch
from SRPGW.GetGLADExMatch import GetGLADExMatch
from SRPGW.GetLEDAxMatch import GetLEDAxMatch
from SRPGW.GetPixVal import GetPixVal
from SRPGW.GetHeadVal import GetHeadVal
from SRPGW.GetTest import GetTest
from SRPGW.GetWeights import GetWeights
from SRPGW.GetPixXY import GetPixXY
from SRPGW.GetPanSTARRS import GetPanSTARRS
from SRPGW.GetSkyBotMatch import GetSkyBotMatch
from SRPGW.GetUSNOxMatch import GetUSNOxMatch
from SRPGW.Score import Score
from SRPGW.ReadDate import ReadDate
from astropy.coordinates import SkyCoord
from astropy import units as u
from astropy.coordinates import match_coordinates_sky
from SRP.SRPTime.UT2MJD import UT2MJD
from SRPFITS.GetFWHM import GetFWHM






parser = argparse.ArgumentParser()
parser.add_argument("-c", "--coordcols", action="store", nargs=2, default=(GW.RACOL,GW.DECCOL), help="Coordinate column labels", metavar=('RA','DEC'))
parser.add_argument("-d", "--debug", action="store_true", help="Debug information")
parser.add_argument("-D", "--date", action="store_true", help="Get frame date")
parser.add_argument("-F", "--FWHM", action="store", type=float, help="Compute frame FWHM (size of computation area in pixel", metavar='fwhmarea')
parser.add_argument("-g", "--gaia", action="store", type=float, help="GAIA catalogue query radius (arcsec)", metavar='gaiarad')
parser.add_argument("-G", "--GLADE", action="store", type=float, help="GLADE catalogue query radius (arcsec)", metavar='gladerad')
parser.add_argument("-H", "--headerinfo", action="store", help="Header for DATE", default=GW.DATEHEAD, metavar='datehead')
parser.add_argument("-i", "--inputfile", action="store", help="Input table", required=True, metavar='inputfile')
parser.add_argument("-l", "--leda", action="store", type=float, help="LEDA catalogue query radius (arcsec)", metavar='ledarad')
parser.add_argument("-m", "--minorplanet", action="store", type=float, help="Minor planet search (radius, arcsec)", metavar='minorplanet')
parser.add_argument("-n", "--neighbor", action="store", help="Neighbor search (table)", metavar='neightable')
parser.add_argument("-o", "--outfile", action="store", help="Output table", required=True, metavar='outfile')
parser.add_argument("-p", "--panstarrs", action="store", type=float, help="PanSTARRS catalogue query radius (arcsec)", metavar='spanstarrsrad')
parser.add_argument("-s", "--simbad", action="store", type=float, help="Simbad catalogue query radius (arcsec)", metavar='simbadrad')
parser.add_argument("-S", "--score", action="store_true", help="Compute source score")
parser.add_argument("-t", "--testgooddata", action="store", type=int, help="Compute statistics ibn object area (halfsize (px))", metavar='halfsize')
parser.add_argument("-u", "--usno", action="store", type=float, help="USNO catalogue query radius (arcsec)", metavar='usnorad')
parser.add_argument("-v", "--verbose", action="store_true", help="Fully describe operations")
parser.add_argument("-w", "--weight", action="store_true", help="Get weight frame value")
parser.add_argument("-W", "--weightarea", action="store", nargs=2, type=float, help="Check weight frame in an area around the source (halfsize (px), limvalue)", metavar=('halfsize','limvalue'))
parser.add_argument("-x", "--xypix", action="store_true", help="Get X,Y values")
parser.add_argument("--version", action="version", version=__version__)
options = parser.parse_args()


#
if options.inputfile and options.outfile:
    #
    if not os.path.isfile(options.inputfile):
        parser.error ("Input file does not exist.")
    #
    if not (options.gaia or options.GLADE or options.leda or options.simbad or options.weight or options.date or options.weightarea or options.xypix or options.FWHM or options.neighbor or options.minorplanet or options.score or options.testgooddata or options.panstarrs or options.usno):
        parser.error( "At least one option among gaia/GLADE/leda/simbad/weight/weightarea/date/xypix/FWHM/neighbor/minorplanet/score/testgooddata/panstarrs/usno should be selected.")
    #
    if options.simbad is not None and options.simbad <= 0.0:
        parser.error ("Simbad query radius must be positive.")
    #
    if options.gaia is not None and options.gaia <= 0.0:
        parser.error ("GAIA query radius must be positive.")
    #
    if options.GLADE is not None and (options.GLADE <= 0.0 and options.GLADE > 180.0):
        parser.error ("GLADE query radius must be positive and lower than 180 arcsec.")
    #
    if options.leda is not None and (options.leda <= 0.0 and options.leda > 180.0):
        parser.error ("LEDA query radius must be positive and lower than 180 arcsec.")
    #
    if options.FWHM is not None and options.FWHM <= 0.0:
        parser.error ("FWHM computation size must be positive.")
    #
    if options.panstarrs is not None and options.panstarrs <= 0.0:
        parser.error ("PanSTARRS query radius must be positive.")
    #
    if options.usno is not None and options.usno <= 0.0:
        parser.error ("USNO query radius must be positive.")
    #
    if options.weightarea:
        try:
            halfsize = int(round(options.weightarea[0]))
        except TypeError:
            parser.error ("Weight area halfsize not correct.")
        weilimit = options.weightarea[1]
    #
    if options.testgooddata:
        thalfsize = options.testgooddata
        if thalfsize <= 0:
                parser.error ("Testgooddata area halfsize not correct.")
    #
    if options.minorplanet:
        if options.minorplanet <= 0:
            parser.error("Search radius for minor planets must be positive.")
    #
    if options.verbose:
        print ("Reading table ", options.inputfile)
    dt = Table.read(options.inputfile, format='ascii.ecsv')
    #
    if (options.coordcols[0]+'_1' in dt.columns.keys()) and (options.coordcols[1]+'_1' in dt.columns.keys()):
        racol = options.coordcols[0]+'_1'
        deccol = options.coordcols[1]+'_1'
    elif (options.coordcols[0] in dt.columns.keys()) and (options.coordcols[1] in dt.columns.keys()):
        racol = options.coordcols[0]
        deccol = options.coordcols[1]
    else:
        parser.error ("Coordinate columns not recognized.")
    #
    if options.neighbor:
        if options.verbose:
            print ("Neighbor search...")
        #
        if not os.path.isfile(options.neighbor):
            parser.error ("Neighborhood table does not exist.")
        #
        if options.verbose:
            print ("Reading table ", options.neighbor)
        nt = Table.read(options.neighbor, format='ascii.ecsv')
        #
        cdt = SkyCoord(ra=dt[racol], dec=dt[deccol])
        cnt = SkyCoord(ra=nt[racol], dec=nt[deccol])
        #
        idx, d2d, d3d = cdt.match_to_catalog_sky(cnt,nthneighbor=2)
        dt[GW.NEICOL] = d2d.arcsec
        if GW.SRPMAG+'_1' in dt.columns.keys():
            dt[GW.DNEICOL] = dt[GW.SRPMAG+'_1'][idx] - dt[GW.SRPMAG+'_1']
    #
    if options.simbad:
        if options.verbose:
            print ("Simbad query...")
        dt = GetSimbadxMatch(dt,options.simbad,options.coordcols)
    #
    if options.gaia:
        if options.verbose:
            print ("GAIA query...")
        dt = GetGAIAxMatch(dt,options.gaia,options.coordcols)
    #
    if options.GLADE:
        if options.verbose:
            print ("GLADE query...")
        dt = GetGLADExMatch(dt,options.GLADE,options.coordcols)
    #
    if options.leda:
        if options.verbose:
            print ("LEDA query...")
        dt = GetLEDAxMatch(dt,options.leda,options.coordcols)
    #
    if options.panstarrs:
        if options.verbose:
            print ("PanSTARRS query...")
        res = GetPanSTARRS(dt[racol],dt[deccol],options.panstarrs)
        if GW.PANSTARRS in dt.columns.keys():
            dt.remove_column(GW.PANSTARRS)
        dt[GW.PANSTARRS] = res
    #
    if options.usno:
        if options.verbose:
            print ("USNO query...")
        dt = GetUSNOxMatch(dt,options.usno,options.coordcols)
    #
    if options.score:
        if options.verbose:
            print ("Score evaluation...")
        #
        if GW.CHI2 in dt.columns.keys():
            schi2 = dt[GW.MINMAG]
        else:
            schi2 = -99
        #
        if GW.LEDACOL in dt.columns.keys():
            sleda = dt[GW.LEDACOL]
        else:
            sleda = -99
        #
        if GW.MINMAG in dt.columns.keys() and GW.VARINDCOL in dt.columns.keys() and GW.NEICOL in dt.columns.keys() and GW.DNEICOL in dt.columns.keys():
            dt[GW.SCOCOL] = Score(dt[GW.MINMAG],dt[GW.VARINDCOL],dt[GW.NEICOL],dt[GW.DNEICOL],sleda,schi2)
        else:
            if options.verbose:
                print ("Some of the required information for score evaluation are missing.")
            dt[GW.SCOCOL] = 0.0
    #
    fnamesf = []
    fnamesw = []
    if options.minorplanet or options.weight or options.weightarea or options.date or options.FWHM or options.xypix or options.testgooddata:
        for e in dt.columns.keys():
            if GW.FNAMECOLF in e:
                fnamesf.append(e)
        if len(fnamesf) == 0:
            parser.error('FITS filename column(s) not present.')
        #
        for e in dt.columns.keys():
            if GW.FNAMECOLW in e:
                fnamesw.append(e)
        if len(fnamesw) == 0 or len(fnamesf) != len(fnamesw):
            parser.error('Weight filename column(s) not present.')
    #
    for fn in range(len(fnamesf)):
        try:
            nfn = fnamesf[fn].split('_')[1]
            nfn = '_'+nfn
        except IndexError:
            nfn = ''
        #
        if options.minorplanet:
            if options.verbose:
                print ("Minor planet queries: epoch %s..." % fnamesf[fn])
            fname = dt[fnamesf[fn]][0]
            if os.path.isfile(fname):
                if GW.DATECOL+nfn in dt.columns.keys():
                    jd = dt[GW.DATECOL+nfn][0] + 2400000.5
                    res = GetSkyBotMatch(dt[racol],dt[deccol],fname,jd,options.minorplanet,options.minorplanet)
                else:
                    parser.error("For minor planet search MJD are required.")
            else:
                res = 'No'
            if GW.MPCOL+nfn in dt.columns.keys():
                dt.remove_column(GW.MPCOL+nfn)
            dt[GW.MPCOL+nfn] = res
        #
        if options.weight:
            if options.verbose:
                print ("Weight queries: epoch %s..." % fnamesw[fn])
            fname = dt[fnamesw[fn]][0]
            if os.path.isfile(fname):
                res = GetPixVal(dt[racol],dt[deccol],fname)
            else:
                res = -99
            if GW.WEIGHTCOL+nfn in dt.columns.keys():
                dt.remove_column(GW.WEIGHTCOL+nfn)
            dt[GW.WEIGHTCOL+nfn] = res
        #
        if options.weightarea:
            if options.verbose:
                print ("Weight area queries: epoch %s..." % fnamesw[fn])
            fname = dt[fnamesw[fn]][0]
            if os.path.isfile(fname):
                res = GetWeights(dt[racol],dt[deccol],fname,weilimit,halfsize)
            else:
                res = -99
            if GW.WEIGHTACOL+nfn in dt.columns.keys():
                dt.remove_column(GW.WEIGHTACOL+nfn)
            dt[GW.WEIGHTACOL+nfn] = res
        #
        if options.testgooddata:
            if options.verbose:
                print ("Testgooddata queries: epoch %s..." % fnamesf[fn])
            fname = dt[fnamesf[fn]][0]
            if os.path.isfile(fname):
                res = GetTest(dt[racol],dt[deccol],fname,thalfsize)
            else:
                res = -99
            if GW.TESTCOL+nfn in dt.columns.keys():
                dt.remove_column(GW.TESTCOL+nfn)
            dt[GW.TESTCOL+nfn] = res
        #
        if options.date:
            if options.verbose:
                print ("Date queries: epoch %s..." % fnamesf[fn])
            fname = dt[fnamesf[fn]][0]
            if os.path.isfile(fname):
                res = GetHeadVal(options.headerinfo,fname)
            else:
                res = -99
            #
            valuedate = ReadDate(res,options.headerinfo)
            #
            if GW.DATECOL+nfn in dt.columns.keys():
                dt.remove_column(GW.DATECOL+nfn)
            dt[GW.DATECOL+nfn] = valuedate
        #
        if options.FWHM:
            if options.verbose:
                print ("FWHM queries: epoch %s..." % fnamesf[fn])
            if GW.XCOL+nfn in dt.columns.keys() and GW.YCOL+nfn in dt.columns.keys():
                fname = dt[fnamesf[fn]][0]
                if os.path.isfile(fname):
                    res = GetFWHM(dt[GW.XCOL+nfn],dt[GW.YCOL+nfn],fname,size=options.FWHM)
                else:
                    res = -99
                if GW.FWCOL+nfn in dt.columns.keys():
                    dt.remove_column(GW.FWCOL+nfn)
                dt[GW.FWCOL+nfn] = res
                dt.meta[GW.FWCOL+nfn] = "%.2f" % numpy.median(dt[GW.FWCOL+nfn])
            else:
                parser.error("Pixel coordinates not present in table.")
        #
        if options.xypix:
            if options.verbose:
                print ("X,Y queries: epoch %s..." % fnamesf[fn])
            fname = dt[fnamesf[fn]][0]
            if os.path.isfile(fname):
                res = GetPixXY(dt[racol],dt[deccol],fname)
            else:
                res = -99,-99
            if GW.XCOL+nfn in dt.columns.keys():
                dt.remove_column(GW.XCOL+nfn)
            if GW.YCOL+nfn in dt.columns.keys():
                dt.remove_column(GW.YCOL+nfn)
            dt[GW.XCOL+nfn] = res[0]
            dt[GW.YCOL+nfn] = res[1]
        #
    dt.write(options.outfile,format='ascii.ecsv', overwrite=True)
    if options.verbose:
        print ("Table %s with %d entries saved." % (options.outfile, len(dt)))
    else:
        print (options.outfile, len(dt))
else:
    parser.print_help()
#
