from astropy.io import fits as pyfits
import numpy
import pyvo


def stringify(s):
    """returns s utf-8-decoded if it's bytes, s otherwise.

    (that's working around old astropy votable breakage)
    """
    if isinstance(s, bytes):
        return s.decode("utf-8")
    return s


def get_image_metadata(field):
    """returns metadata rows from the ppakm31 service for field.

    Access URL and schema are taken from a TOPCAT exploration of the service.
    """
    svc = pyvo.dal.TAPService("http://dc.zah.uni-heidelberg.de/tap")

    res = []
    for row in svc.run_sync(
            "select accref, field, imagetitle, bandpassid, cube_link"
            "  from ppakm31.maps"
            "  where field={}".format(repr(field))):
        res.append({
            "bandpassid": stringify(row["bandpassid"]),
            "accref": stringify(row["accref"]),})

    return res


def get_line_maps(image_meta):
    """returns a dict band->hdus for image_meta as returned by
    get_image_maps.
    """
    res = {}
    for row in image_meta:
        if row["bandpassid"]=='[NII]6583':
             # we don't need that one later, so skip it
             continue
        res[row["bandpassid"]] = pyfits.open(row["accref"])
    return res


def get_kewley_map(line_maps):
    """does the math to compute the Kewley criterion, 2001ApJ...556..121K.

    Line_maps is a dict line label -> fits hdus.
    """
    return (numpy.log(
            line_maps["[OIII]5007"][0].data
            /line_maps["Hbeta"][0].data)
        /(0.72/(
            numpy.log((
                    line_maps["[SII]6716"][0].data
                    +line_maps["[SII]6730"][0].data)
                /line_maps["Halpha"][0].data-0.32))
            +1.30))


if __name__=="__main__":
    line_maps = get_line_maps(get_image_metadata("F2"))
    kewley_map = get_kewley_map(line_maps)

    # hack: re-use Halpha for our line map, with just a hint of
    # sanitation
    new_hdr = line_maps["Halpha"][0].header
    for outdated_card in ["CRPIX3", "CRVAL3", "CDELT3", "BUNIT",
            "CUNIT3", "DATAMIN", "DATAMAX"]:
        del new_hdr[outdated_card]
    new_hdr["HISTORY"] = "butchered by build_i_map.py"
    pyfits.HDUList([
        pyfits.PrimaryHDU(header=new_hdr, data=kewley_map)]
        ).writeto("kewley_map.fits", overwrite=True)