How to use drizzlepac - 10 common examples

pixel variances.
    
    ds9: `pyds9.DS9`
        Show intermediate steps of the drizzling
    
    Returns
    -------
    hdu : `~astropy.io.fits.HDUList`
        FITS HDUList with the drizzled 2D spectrum and weight arrays
        
    """
    from drizzlepac.astrodrizzle import adrizzle
    from astropy import log
    log.setLevel('ERROR')
    #log.disable_warnings_logging()
    adrizzle.log.setLevel('ERROR')
    
    NX = int(np.round(np.diff(wlimit)[0]*1.e4/dlam)) // 2
    center = np.mean(wlimit[:2])*1.e4
    out_header, output_wcs = utils.make_spectrum_wcsheader(center_wave=center,
                                 dlam=dlam, NX=NX, 
                                 spatial_scale=spatial_scale, NY=NY)
                                 
    sh = (out_header['NAXIS2'], out_header['NAXIS1'])
    
    outsci = np.zeros(sh, dtype=np.float32)
    outwht = np.zeros(sh, dtype=np.float32)
    outctx = np.zeros(sh, dtype=np.int32)
    
    outvar = np.zeros(sh, dtype=np.float32)
    outwv = np.zeros(sh, dtype=np.float32)
    outcv = np.zeros(sh, dtype=np.int32)

def apply_tweak_shifts(wcs_ref, shift_dict, grism_matches={}, verbose=True):
    """
    
    """
    from drizzlepac import updatehdr

    hdu = wcs_ref.to_fits(relax=True)
    file0 = list(shift_dict.keys())[0].split('.fits')[0]
    tweak_file = '{0}_tweak_wcs.fits'.format(file0)
    hdu.writeto(tweak_file, clobber=True)
    for file in shift_dict:
        updatehdr.updatewcs_with_shift(file, tweak_file,
                                        xsh=shift_dict[file][0],
                                        ysh=shift_dict[file][1],
                                        rot=0., scale=1.,
                                        wcsname='SHIFT', force=True,
                                        reusename=True, verbose=verbose,
                                        sciext='SCI')
        
        ### Bug in astrodrizzle? Dies if the FLT files don't have MJD-OBS
        ### keywords
        im = pyfits.open(file, mode='update')
        im[0].header['MJD-OBS'] = im[0].header['EXPSTART']
        im.flush()
        
        # Update paired grism exposures
        if file in grism_matches:
            for grism_file in grism_matches[file]:

for drz_file in drz_files[1:]:
        root = drz_file.split('_drz')[0]
        result = align_drizzled_image(root=root, mag_limits=mag_limits,
                                      radec=rd_ref,
                                      NITER=5, clip=20)

        orig_wcs, drz_wcs, out_shift, out_rot, out_scale = result
        
        im = pyfits.open(drz_file)
        files = []
        for i in range(im[0].header['NDRIZIM']):
          files.append(im[0].header['D{0:03d}DATA'.format(i+1)].split('[')[0])
        
        
        for file in files:
            updatehdr.updatewcs_with_shift(file, drz_files[0],
                                      xsh=out_shift[0], ysh=out_shift[1],
                                      rot=out_rot, scale=out_scale,
                                      wcsname=ref_catalog, force=True,
                                      reusename=True, verbose=True,
                                      sciext='SCI')

            im = pyfits.open(file, mode='update')
            im[0].header['MJD-OBS'] = im[0].header['EXPSTART']
            im.flush()

        ### Second drizzle
        if len(files) > 1:
            AstroDrizzle(files, output=root, clean=True, context=False, preserve=False, skysub=True, driz_separate=False, driz_sep_wcs=False, median=False, blot=False, driz_cr=False, driz_cr_corr=False, driz_combine=True, final_bits=576, coeffs=True, resetbits=0, build=False, final_wht_type='IVM')        
        else:
            AstroDrizzle(files, output=root, clean=True, final_scale=None, final_pixfrac=1, context=False, final_bits=576, preserve=False, driz_separate=False, driz_sep_wcs=False, median=False, blot=False, driz_cr=False, driz_cr_corr=False, driz_combine=True, build=False, final_wht_type='IVM')

### Contamination-cleaned
        adrizzle.do_driz(beam_data, beam_wcs, wht, output_wcs, 
                         outsci, outwht, outctx, 1., 'cps', 1,
                         wcslin_pscale=beam.grism.wcs.pscale, uniqid=1, 
                         pixfrac=pixfrac, kernel=kernel, fillval=0, 
                         stepsize=10, wcsmap=None)
        
        ### Continuum
        adrizzle.do_driz(beam_continuum, beam_wcs, wht, output_wcs, 
                         coutsci, coutwht, coutctx, 1., 'cps', 1, 
                         wcslin_pscale=beam.grism.wcs.pscale, uniqid=1, 
                         pixfrac=pixfrac, kernel=kernel, fillval=0, 
                         stepsize=10, wcsmap=None)
        
        ### Contamination
        adrizzle.do_driz(beam.contam, beam_wcs, wht, output_wcs, 
                         xoutsci, xoutwht, xoutctx, 1., 'cps', 1, 
                         wcslin_pscale=beam.grism.wcs.pscale, uniqid=1, 
                         pixfrac=pixfrac, kernel=kernel, fillval=0, 
                         stepsize=10, wcsmap=None)
        
        ### Direct thumbnail
        if direct_extension == 'REF':
            if beam.direct['REF'] is None:
                direct_extension = 'SCI'
                
        if direct_extension == 'REF':
            thumb = beam.direct['REF']
            thumb_wht = np.cast[np.float32]((thumb != 0)*1)
        else:
            thumb = beam.direct[direct_extension]#/beam.direct.photflam
            thumb_wht = 1./(beam.direct.data['ERR']/beam.direct.photflam)**2

wht = 1./im['ERR'].data**2
        wht[(im['ERR'].data == 0) | (dq > 0) | (~np.isfinite(wht))] = 0
        wht[im['SCI'].data < -3*im['ERR'].data] = 0

        wht = np.cast[np.float32](wht)

        exp_wcs = pywcs.WCS(im[1].header, relax=True)
        exp_wcs.pscale = utils.get_wcs_pscale(exp_wcs)

        #pf = 0.5
        # import drizzlepac.wcs_functions as dwcs
        # xx = out_wcs.deepcopy()
        # #xx.all_pix2world = xx.wcs_world2pix
        # map = dwcs.WCSMap(exp_wcs, xx)

        astrodrizzle.adrizzle.do_driz(im['SCI'].data, exp_wcs, wht,
                                      inter_wcs, outsci, outwht, outctx,
                                      1., 'cps', 1,
                                      wcslin_pscale=exp_wcs.pscale,
                                      uniqid=1,
                                      pixfrac=pixfrac, kernel=kernel,
                                      fillval=0, stepsize=10,
                                      wcsmap=SIP_WCSMap)

        if ds9 is not None:
            ds9.view(outsci, header=outh)

    #outsci /= out_wcs.pscale**2
    rms = 1/np.sqrt(outwht)
    mask = (outwht == 0) | (rms > 100)
    rms[mask] = 0
    outsci[mask] = 0.

line_wcs = pywcs.WCS(h, relax=True)
            line_wcs.pscale = utils.get_wcs_pscale(line_wcs)

            # Science and wht arrays
            sci = flt.grism['SCI'] - flt.model
            wht = 1/(flt.grism['ERR']**2)
            scl = np.exp(-(fcontam*np.abs(flt.model)/flt.grism['ERR']))
            wht *= scl
            
            wht[~np.isfinite(wht)] = 0
            
            # Drizzle it
            if verbose:
                print('Drizzle {0} to wavelength {1:.2f}'.format(flt.grism.parent_file, wave))
                                                        
            adrizzle.do_driz(sci, line_wcs, wht, out_wcs, 
                             outsci, outwht, outctx, 1., 'cps', 1,
                             wcslin_pscale=line_wcs.pscale, uniqid=1, 
                             pixfrac=pixfrac, kernel=kernel, fillval=0, 
                             stepsize=10, wcsmap=None)
        
        # Done!
        return outsci, outwht

wht[~np.isfinite(data_i+scl)] = 0
            contam_weight[~np.isfinite(data_i+scl)] = 0
            data_i[~np.isfinite(data_i+scl)] = 0
        
        ###### Go drizzle
        
        ### Contamination-cleaned
        adrizzle.do_driz(data_i, beam_wcs, wht, output_wcs, 
                         outsci, outwht, outctx, 1., 'cps', 1,
                         wcslin_pscale=1.0, uniqid=1, 
                         pixfrac=pixfrac, kernel=kernel, fillval=0, 
                         stepsize=10, wcsmap=None)
        
        # For variance
        adrizzle.do_driz(contam_weight, beam_wcs, wht, output_wcs, 
                         outvar, outwv, outcv, 1., 'cps', 1,
                         wcslin_pscale=1.0, uniqid=1, 
                         pixfrac=pixfrac, kernel=kernel, fillval=0, 
                         stepsize=10, wcsmap=None)
        
        if ds9 is not None:
            ds9.view(outsci/output_wcs.pscale**2, header=out_header)
    
    ### Correct for drizzle scaling
    outsci /= output_wcs.pscale**2
    
    # variance
    outvar /= outwv*output_wcs.pscale**2
    outwht = 1/outvar
    outwht[(outvar == 0) | (~np.isfinite(outwht))] = 0

else:
            wht *= sens**2
            beam_data /= sens
            beam_continuum /= sens
        
        ###### Go drizzle
        
        ### Contamination-cleaned
        adrizzle.do_driz(beam_data, beam_wcs, wht, output_wcs, 
                         outsci, outwht, outctx, 1., 'cps', 1,
                         wcslin_pscale=beam.grism.wcs.pscale, uniqid=1, 
                         pixfrac=pixfrac, kernel=kernel, fillval=0, 
                         stepsize=10, wcsmap=None)
        
        ### Continuum
        adrizzle.do_driz(beam_continuum, beam_wcs, wht, output_wcs, 
                         coutsci, coutwht, coutctx, 1., 'cps', 1, 
                         wcslin_pscale=beam.grism.wcs.pscale, uniqid=1, 
                         pixfrac=pixfrac, kernel=kernel, fillval=0, 
                         stepsize=10, wcsmap=None)
        
        ### Contamination
        adrizzle.do_driz(beam.contam, beam_wcs, wht, output_wcs, 
                         xoutsci, xoutwht, xoutctx, 1., 'cps', 1, 
                         wcslin_pscale=beam.grism.wcs.pscale, uniqid=1, 
                         pixfrac=pixfrac, kernel=kernel, fillval=0, 
                         stepsize=10, wcsmap=None)
        
        ### Direct thumbnail
        if direct_extension == 'REF':
            if beam.direct['REF'] is None:
                direct_extension = 'SCI'

dlam = np.interp(wave, beam.beam.lam[1:], np.diff(beam.beam.lam))
        # 1e-17 erg/s/cm2 #, scaling closer to e-/s
        sens *= 1./dlam
        
        if sens == 0:
            continue
        else:
            wht *= sens**2
            beam_data /= sens
            beam_continuum /= sens
        
        ###### Go drizzle
        
        ### Contamination-cleaned
        adrizzle.do_driz(beam_data, beam_wcs, wht, output_wcs, 
                         outsci, outwht, outctx, 1., 'cps', 1,
                         wcslin_pscale=beam.grism.wcs.pscale, uniqid=1, 
                         pixfrac=pixfrac, kernel=kernel, fillval=0, 
                         stepsize=10, wcsmap=None)
        
        ### Continuum
        adrizzle.do_driz(beam_continuum, beam_wcs, wht, output_wcs, 
                         coutsci, coutwht, coutctx, 1., 'cps', 1, 
                         wcslin_pscale=beam.grism.wcs.pscale, uniqid=1, 
                         pixfrac=pixfrac, kernel=kernel, fillval=0, 
                         stepsize=10, wcsmap=None)
        
        ### Contamination
        adrizzle.do_driz(beam.contam, beam_wcs, wht, output_wcs, 
                         xoutsci, xoutwht, xoutctx, 1., 'cps', 1, 
                         wcslin_pscale=beam.grism.wcs.pscale, uniqid=1,

scl = convert_to_flambda*beam.beam.total_flux/1.e-17
            scl *= 1./beam.flat_flam.reshape(beam.beam.sh_beam).sum(axis=0)
            #scl = convert_to_flambda/beam.beam.sensitivity
            
            data_i *= scl
            wht *= (1/scl)**2            
            #contam_weight *= scl
            
            wht[~np.isfinite(data_i+scl)] = 0
            contam_weight[~np.isfinite(data_i+scl)] = 0
            data_i[~np.isfinite(data_i+scl)] = 0
        
        ###### Go drizzle
        
        ### Contamination-cleaned
        adrizzle.do_driz(data_i, beam_wcs, wht, output_wcs, 
                         outsci, outwht, outctx, 1., 'cps', 1,
                         wcslin_pscale=1.0, uniqid=1, 
                         pixfrac=pixfrac, kernel=kernel, fillval=0, 
                         stepsize=10, wcsmap=None)
        
        # For variance
        adrizzle.do_driz(contam_weight, beam_wcs, wht, output_wcs, 
                         outvar, outwv, outcv, 1., 'cps', 1,
                         wcslin_pscale=1.0, uniqid=1, 
                         pixfrac=pixfrac, kernel=kernel, fillval=0, 
                         stepsize=10, wcsmap=None)
        
        if ds9 is not None:
            ds9.view(outsci/output_wcs.pscale**2, header=out_header)
    
    ### Correct for drizzle scaling