def _init_as_astrodata(self):
        """
           Initialize parameters to be used by as_astrodata.

           Creates a WCS object (pywcs) from the SCI header and
           form the output AD object with the PHU and MDF from
           the input AD. We are adding the TRACEFP extension as well 
           for later use on the spectral reduction process. 

           Input:
              self.ad:  AD object.
           Output:
              adout:  Output AD object with AD phu and MDF 
        """

        ad = self.ad
        # Start output AD with the original phu and the MDF extension.
        adout = AstroData(phu=ad.phu)
        adout.append(ad['MDF'])
        adout.append(ad['TRACEFP'])

        # Get wcs information. It is in the PHU
        try:
            self.wcs = pywcs.WCS(ad.phu.header)
            if not hasattr(self.wcs.wcs, 'cd'):
                self.wcs = None
        except:   # Something wrong with WCS, set it to None
            self.wcs = None

        return adout

def test_method_info_4():
    old_stdout = sys.stdout
    sys.stdout = mystdout = StringIO()
    ad = AstroData(TESTFILE)
    ad.info(oid=True, table=True)
    sys.stdout = old_stdout
    assert isinstance(mystdout.getvalue(), str)

 def recover(self):
     log.debug("OufileETIFile recover()")
     ad = AstroData(self.tmp_name, mode="update")
     ad.filename = self.ad_name
     ad = gemini_tools.obsmode_del(ad)
     log.fullinfo(self.tmp_name + " was loaded into memory")
     return ad

def example():
    """
    This is just an example.  Cut and paste that on the python prompt.
    It can also be run as specplot.example().
    """
    import numpy as np
    import matplotlib.pyplot as plt
    from astropy import wcs
    from astrodata import AstroData
    
    ad = AstroData('JHK.fits')
    x_values = np.arange(ad.get_key_value('NAXIS1'))
    
    wcs_ad = wcs.WCS(ad.header.tostring())
    wlen = wcs_ad.wcs_pix2world(zip(x_values), 0)
    
    plt.plot(wlen, ad.data)
    plt.xlabel('Wavelength [Angstrom]')
    plt.ylabel('Counts')
    plt.axis('tight')
    plt.ylim(-100, 800)
    plt.show()
    
    ad.close()

    #plt.axis[[-100,1000,ymin,ymax]]

def atd6():
    """
    Verify that a MosaicAD method can create a block from a given extension name.

    The test creates a mosaic ndarray from the MosaicAD method mosaic_image_data 
    with the block parameter value (0,0), indicating to output the lower left block.

    NOTE: Having one amp per block, the actual extension data is not the same 
          as the block since it would be trim by the DATASEC image section.

    gmos_file='../data/gS20120420S0033.fits'
    gsaoi_file='../data/guS20120413S0048.fits'
    """
        
    from astrodata import AstroData
    from gempy.adlibrary.mosaicAD import MosaicAD
    #    This is the default Mosaic function
    from gempy.mosaic.gemMosaicFunction import gemini_mosaic_function
    print '\n atd6 REQUIREMENT.......'
    print ('***** From a given AstroData object, the system shall create a block from '
          'a given extension name')

    gmos_file='../data/gS20120420S0033.fits'
    gsaoi_file='../data/guS20120413S0048.fits'

    for file in [gmos_file,gsaoi_file]:
        ad = AstroData(file)
        print 'Instrument: ',ad.instrument()
        mo = MosaicAD(ad, gemini_mosaic_function)
        #    Now use the mosaic_image_data method to generate
        #    an output block by using the parameter block and
        #    value as a tuple (col,row) (0-based) of the block
        #    you want returned. For GMOS the block values are
        #    (0,0), (1,0), (2,0).

        block=(1,0)
        block_data = mo.mosaic_image_data(block=block,tile=True)
         
        #    Get the shape: (height, width) in pixels.
        print 'block_data.shape:',block_data.shape

        extn = block[0] + block[1]*mo.geometry.mosaic_grid[0] + 1
        print 'Input shape for 1-amp per detector:',ad['SCI',extn].data.shape

        #    Check values of the lower 2x2 pixels
        print 'Output block [0:2,0:2] pixels:\n',block_data[:2,:2]
        if ad.instrument() == 'GSAOI':
            # GSAOI FITS extension 1 correspond to block (1,0)
            # and extension 2 to block (0,0). DETSEC image section
            # indicates this.
            extn = [2,1,3,4][extn-1]

        # To get the correct segment we need to look at
        # DATASEC, in case the data is trimm -as it appears in data_list.
        x1,x2,y1,y2 = ad.data_section().as_dict()['SCI',extn]
        print 'Input amp DATASEC[x1:x1+2 pixels:\n',\
              ad['SCI',extn].data[x1:x1+2,y1:y1+2]
        print '\n'

def test_method_insert_7():
    ad = AstroData(TESTFILE)
    xname = "TEST"
    xver  = 99
    del header1['EXTNAME']
    ad.insert(1, header=header1, data=data1, extname=xname,
              extver=xver, auto_number=True)
    assert ad[1].header.get("EXTNAME") == xname
    assert ad[1].header.get("EXTVER") == xver

def test_finding_an_easier_center():
    ad = AstroData(get_data_file_name("N20060131S0012.fits"))
    selection = get_selection_peak(ad.data, (489.07, 478.99), float(ad.pixel_scale()))

    predicted_center = selection.get_center()

    assert_tolerance(predicted_center,
                     (489.11025833697016, 478.68088198208636),
                     tolerance=0.02)

def test_finding_the_center_of_titan():
    ad = AstroData(get_data_file_name("N20060131S0011.fits"))
    selection = get_selection_peak(ad.data, (391.0, 539.0), float(ad.pixel_scale()))

    predicted_center = selection.get_center()

    # the object is titan, doesn't have a very clearly defined center
    assert_tolerance(predicted_center,
                     (384.87060478881705, 542.73266988305159),
                     tolerance=0.07)

def test_out_of_bound_aperture():
    ad = AstroData(get_data_file_name("N20060131S0012.fits"))
    selection = get_selection_peak(ad.data, (10.0, 10.0), float(ad.pixel_scale()))

    predicted_center = selection.get_center()

    # should be rather non-sense that is returned
    assert_tolerance(predicted_center,
                     (50.0, 50.0),
                     tolerance=50.0)

 def recover(self):
     log.debug("OutAtList recover()")
     adlist = []
     for i, tmpname in enumerate(self.diskoutlist):
         ad = AstroData(tmpname, mode="update")
         ad.filename = self.ad_name[i]
         ad = gemini_tools.obsmode_del(ad)
         adlist.append(ad)
         log.fullinfo(tmpname + " was loaded into memory")
     return adlist

    def as_astrodata(self):
        """
            
          With each cut object in the cut_list having the SCI,DQ,VAR set,
          form an hdu and append it to adout.  Update keywords EXTNAME= 'SCI', 
          EXTVER=<footprint#>, CCDSEC, DISPAXIS, CUTSECT, CUTORDER in the header
          and reset WCS information if there was a WCS in the input AD header.

          ::

           Input:
              self.cut_list: List of Cut objects.
              self.adout:    Output AD object with MDF and
                             TRACEFP extensions.
           Output:
              adout: contains the appended HDUs.
        """

        adout = self._init_as_astrodata()

        ad = self.ad
        scihdr =        ad['SCI',1].header.copy()
        if self.has_dq:
            dqheader =  ad['DQ', 1].header.copy()
        if self.has_var:
            varheader = ad['VAR',1].header.copy()

        # Update NSCIEXT keyword to represent the current number of cuts.
        if new_pyfits_version:
            adout.phu.header.update = adout.phu.header.set
        adout.phu.header.update('NSCIEXT',len(self.cut_list)) 

        # This is a function renaming when using Pyfits 3.1
        if new_pyfits_version:
            scihdr.update = scihdr.set
        extver = 1

        # Generate the cuts using the region's sci_cut,var_cut and
        # dq_cut
        for region,sci_cut,var_cut,dq_cut in self.cut_list: 
            rx1,rx2,ry1,ry2 = np.asarray(region) + 1   # To 1-based
            csec = '[%d:%d,%d:%d]'%(rx1,rx2,ry1,ry2)
            scihdr.update('NSCUTSEC',csec,
                          comment="Region extracted by 'cut_footprints'")
            scihdr.update('NSCUTSPC',extver,comment="Spectral order")
            form_extn_wcs(scihdr, self.wcs, region)
            new_sci_ext = AstroData(data=sci_cut,header=scihdr)
            new_sci_ext.rename_ext(name='SCI',ver=extver)
            adout.append(new_sci_ext)
            if self.has_dq:
                new_dq_ext = AstroData(data=dq_cut, header=dqheader)
                new_dq_ext.rename_ext(name='DQ',ver=extver)
                adout.append(new_dq_ext)
            if self.has_var:
                new_var_ext = AstroData(data=var_cut, header=varheader)
                new_var_ext.rename_ext(name='VAR',ver=extver)
                adout.append(new_var_ext)
            extver += 1

        return adout

 def write_new_table(self, fname):
     cols = list(self.get_columns())
     cols.extend(self.get_fiber_positions_columns())
     
     # Create the table HDU
     tablehdu = pf.new_table(cols)
     
     # Create an AstroData object to contain the table
     # and write to disk.
     new_ad = AstroData(tablehdu)
     new_ad.rename_ext('SCI', 1)
     new_ad.write(fname, clobber=True)

def check_gmos_image(filepath, calibrations=False):
    from astrodata import AstroData
    reason = "GMOS image"
    try:
        ad = AstroData(filepath)
    except:
        reason = "can't load file"
        return False, reason
    
    try:
        fp_mask = ad.focal_plane_mask().as_pytype()
    except:
        fp_mask = None

    if "GMOS" not in ad.types:
        reason = "not GMOS"
        return False, reason
    elif "GMOS_DARK" in ad.types:
        reason = "GMOS dark"
        return False, reason
    elif "GMOS_BIAS" in ad.types and not calibrations:
        reason = "GMOS bias"
        return False, reason
    elif "GMOS_IMAGE_FLAT" in ad.types and not calibrations:
        reason = "GMOS flat"
        return False, reason
    elif ("GMOS_IMAGE" in ad.types and
          fp_mask!="Imaging"):
        reason = "GMOS slit image"
        return False, reason
    elif (("GMOS_IMAGE" in ad.types and
           fp_mask == "Imaging" and
           "GMOS_DARK" not in ad.types) or
          "GMOS_BIAS" in ad.types or 
          "GMOS_IMAGE_FLAT" in ad.types):

        # Test for 3-amp mode with e2vDD CCDs
        # This mode has not been commissioned.
        dettype = ad.phu_get_key_value("DETTYPE")
        if dettype == "SDSU II e2v DD CCD42-90":
            namps = ad.phu_get_key_value("NAMPS")
            if namps is not None and int(namps)==1:
                reason = "uncommissioned 3-amp mode"
                return False, reason
            else:
                return True, reason
        else:
            return True, reason
    else:
        reason = "not GMOS image"
        return False, reason

 def recover(self):
     log.debug("OutAtList recover()")
     adlist = []
     for i, tmpname in enumerate(self.diskoutlist):
         ad = AstroData(tmpname, mode="update")
         ad.filename = self.ad_name[i]
         ad = gemini_tools.obsmode_del(ad)
         # Read the database back in, if it exists
         try:
             ad = gemini_tools.read_database(
                 ad, database_name=self.database_name, 
                 input_name=self.tmpin_name[i], 
                 output_name=ad.phu_get_key_value("ORIGNAME"))
         except:
             pass
         adlist.append(ad)
         log.fullinfo(tmpname + " was loaded into memory")
     return adlist

def check_gmos_longslit(filepath):
    from astrodata import AstroData
    reason = "GMOS longslit"

    try:
        ad = AstroData(filepath)
    except:
        reason = "can't load file"
        return False, reason
    
    try:
        fp_mask = ad.focal_plane_mask().as_pytype()
    except:
        fp_mask = None

    if "GMOS" not in ad.types:
        reason = "not GMOS"
        return False, reason
    elif "GMOS_DARK" in ad.types:
        reason = "GMOS dark"
        return False, reason
    elif "GMOS_BIAS" in ad.types:
        reason = "GMOS bias"
        return False, reason
    elif "GMOS_NODANDSHUFFLE" in ad.types:
        reason = "GMOS nod and shuffle"
        return False, reason
    elif "GMOS_LS" in ad.types:

        # Test for 3-amp mode with e2vDD CCDs
        # This mode has not been commissioned.
        dettype = ad.phu_get_key_value("DETTYPE")
        if dettype == "SDSU II e2v DD CCD42-90":
            namps = ad.phu_get_key_value("NAMPS")
            if namps is not None and int(namps)==1:
                reason = "uncommissioned 3-amp mode"
                return False, reason
            else:
                return True, reason
        else:
            return True, reason
    else:
        reason = "not GMOS longslit"
        return False, reason

def runinsert(index=None, f1=None, f2=None, auto=False):
    ad = AstroData(f1)
    md = AstroData(f2)
    pstr = "\n\n             >>>>>>>     AD     <<<<<<<<\n"
    pstr += str(ad.infostr())
    pstr += "\n\n             >>>>>>>    AD APPEND   <<<<<<<<\n"
    pstr += str(md.infostr())
    ad.insert(index=index, moredata=md, auto_number=auto)
    pstr +="\n\n             >>>>>>>  NEW AD <<<<<<<<\n"
    pstr += str(ad.infostr())
    print(pstr)
    return ad

def atd1():
    """
    With a GMOS AstroData object, the test instantiates a MosaicAD object
    containing 'coords' as one of the attributes. The test verify that 
    coords['amp_mosaic_coord'] and ad['SCI'].detector_array.as_dict() values
    match.
    """
   
    print '\n atd1 REQUIREMENT.......'
    print ('*****MosaicAD shall instantiate an object from a supported AstroData objec')
 
    gmos_file = '../data/gS20120420S0033.fits'
    gsaoi_file = '../data/guS20110324S0146.fits'
    nongem_file='../data/kp620765.fits'
        
    from astrodata import AstroData
    from gempy.adlibrary.mosaicAD import MosaicAD
    #     This is the default Mosaic function
    from gempy.mosaic.gemMosaicFunction import gemini_mosaic_function

    #    Success Criteria 1. (GMOS data)
    #    The tester should provide her/his own GMOS file.
    for file in [gmos_file,gsaoi_file,nongem_file]:
        ad = AstroData(file)
        print '\n ...........CASE for:',file,ad.instrument()

        #    Create the MosaicAD object
        mo = MosaicAD(ad,gemini_mosaic_function)

        #    print DETECTOR values for all the 'SCI' extensions as
        #    a dictionary.
        print ad['SCI'].detector_section().as_dict()
          
        #    print the 'amp_mosaic_coord' key value from the 'coords'
        #    attribute. This list is in increasing order of extver.
        print mo.coords['amp_mosaic_coord']

    def __init__(self, filename, mosmask=None, mdfdir=None):
        self.ad = AstroData(filename)
        self.mosmask = mosmask
        self.mdfdir = mdfdir

        # Determine extension
        nsci = len(self.ad)
        debug("...nsci = ", nsci)
            
        if nsci > 1:
            l_sci_ext = 1 
        else:
            l_sci_ext = 0

        debug("...using extension [" + str(l_sci_ext) + "]")

        overscan_dv = self.ad[l_sci_ext].overscan_section()

        if self.is_mos_mode():
            self.box_coords = parse_box_coords(self, self.get_mdf_filename())
            self.box_mosaic = BoxMosaic(self, self.box_coords)
            self.scidata = self.box_mosaic.get_science_data()
        elif self.is_new_gmosn_ccd():
            # tile the 2 center parts of the new GMOS image
            self.scidata = gmultiamp(self.ad)
        elif not overscan_dv.is_none():
            # remove the overscan so we don't have to take it into account when guessing the slit location
            self.scidata = subtract_overscan(self.ad[l_sci_ext])

            # it still affects the center of rotation however
            ox1, ox2, oy1, oy2 = overscan_dv.as_list()
            correction = np.array([ox2 - ox1, 0])
            center = self.get_binned_data_center() - correction
            self.fieldcenter = center * self.detector_y_bin()
        else:
            self.scidata = self.ad[l_sci_ext].data

def test_method_phu_get_key_val_2():
    ad = AstroData(TESTFILE2)
    assert isinstance(ad.phu_get_key_value('BITPIX'), int)

# it would be useful for the DA/SOSs.
# It can be run on any file with a sextractor style OBJCAT in it, for example
# the _forStack files output by the QAP. If no REFCAT is present, it bails out as no refmags
#
# Paul Hirst 20120321
import math
import sys
from astrodata import AstroData
import numpy as np
import matplotlib.pyplot as plt

from random import random

filename = sys.argv[1]

ad = AstroData(filename)
objcat = ad['OBJCAT']

if (ad['REFCAT'] is None):
  print "No Reference Catalog in this file, thus no Zeropoints. Sorry"
  sys.exit(0)

mag = objcat.data.field("MAG_AUTO")
magerr = objcat.data.field("MAGERR_AUTO")
refmag = objcat.data.field("REF_MAG")
refmagerr = objcat.data.field("REF_MAG_ERR")
sxflag = objcat.data.field("FLAGS")
dqflag = objcat.data.field("IMAFLAGS_ISO")

# set mag to None where we don't want to use the object
mag = np.where((mag==-999), None, mag)