How to use the astropy.wcs.WCS function in astropy

def _load_single_6dfgs_hdu(hdu):
    """
    Helper function to handle loading spectra from a single 6dFGS HDU
    """

    header = hdu.header
    w = WCS(naxis=1)
    w.wcs.crpix[0] = header["CRPIX1"]
    w.wcs.crval[0] = header["CRVAL1"]
    w.wcs.cdelt[0] = header["CDELT1"]
    w.wcs.cunit[0] = header["CUNIT1"]
    if w.wcs.cunit[0] == ANGSTROMS:
        w.wcs.cunit[0] = Unit("Angstrom")
    meta = {"header": header}
    flux = hdu.data[0] * Unit("count") / w.wcs.cunit[0]
    uncertainty = VarianceUncertainty(hdu.data[1])

    return Spectrum1D(flux=flux, wcs=w, meta=meta, uncertainty=uncertainty)

----------
        centers: list of centers for each camera/chip FITS file
        pos: [RA,Dec] position of the source 
    Returns
    ----------
        centers[i][0], centers[i][1]: the camera and chip numbers. Returns 0, 0 if source 
        is not found in any file
    """
    for i in np.arange(1,5,1):
        for j in np.arange(1,5,1):
            dir  = './2019/2019_1_{}-{}/ffis/'.format(i, j)
            file = 'tess2019132000826-{}-{}-0016-s_ffic.fits'.format(i, j)
            if i == 3 and (j == 2 or j == 3):
                file = 'tess2019130000826-{}-{}-0016-s_ffic.fits'.format(i, j)
            mast, mheader = fits.getdata(dir+file, header=True)
            xy = WCS(mheader).all_world2pix(pos[0], pos[1], 1, quiet=True)
            if xy[0] >= 0. and xy[0] <= len(mast) and xy[1] >= 0. and xy[1] <= len(mast[0]):
                return dir, xy, i, j

def _define_wcs(self):
        """
        Given what we know about the scale of the image,
        define a nearly-correct world coordinate system to use with it.
        """
        w = WCS(naxis=2)
        w.wcs.crpix = self.size_pix / 2
        w.wcs.crval = self.center
        w.wcs.cd = np.array([[-1, 0], [0, 1]]) * self.scale / 3600.
        w.wcs.ctype = ['RA---TAN', 'DEC--TAN']
        w.wcs.cunit = ['deg', 'deg']
        w.wcs.radesys = 'ICRS'
        w.wcs.equinox = 2000.0
        self.wcs = w

"""
    import astropy.wcs as pywcs

    new_header = header.copy()

    if 'TELESCOP' in new_header:
        if new_header['TELESCOP'] != 'JWST':
            keys = ['TELESCOP', 'FILTER', 'DETECTOR', 'INSTRUME']
            for key in keys:
                if key in header:
                    new_header.remove(key)

    if simplify_wcs:
        # Make simple WCS header
        orig_wcs = pywcs.WCS(new_header)
        new_header = orig_wcs.to_header()

        new_header['EXTNAME'] = 'SCI'
        new_header['RADESYS'] = 'ICRS'
        new_header['CDELT1'] = -new_header['PC1_1']
        new_header['CDELT2'] = new_header['PC2_2']
        new_header['PC1_1'] = -1
        new_header['PC2_2'] = 1

    return new_header

self.cd=np.array([[hdu[0].header['CD1_1'],hdu[0].header['CD1_2']],
                              [hdu[0].header['CD2_1'],hdu[0].header['CD2_2']]])
            self.ctype=[hdu[0].header['CTYPE1'],hdu[0].header['CTYPE2']]
            self.cunit=[hdu[0].header['CUNIT1'],hdu[0].header['CUNIT2']]
            self.crres=np.array([hdu[0].header['CRRES1'],hdu[0].header['CRRES2']])

        # Compute image properties
        self.sx=np.sqrt(self.cd[0,0]**2+self.cd[1,0]**2)
        self.sy=np.sqrt(self.cd[0,1]**2+self.cd[1,1]**2)
        self.wx=self.nx*self.sx
        self.wy=self.ny*self.sy
        self.vmin=np.mean(self.zmax)-2.0*np.std(self.zmax)
        self.vmax=np.mean(self.zmax)+6.0*np.std(self.zmax)

        # Setup WCS
        self.w=wcs.WCS(naxis=2)
        self.w.wcs.crpix=self.crpix
        self.w.wcs.crval=self.crval
        self.w.wcs.cd=self.cd
        self.w.wcs.ctype=self.ctype
        self.w.wcs.set_pv([(2,1,45.0)])

pixlim_coms = ["Start of valid flux calibration",
                   "End of valid flux calibration",
                   'Start of region with at least one fiber',
                   'End of region with at least one fiber',
                   'Start of region with all fibers',
                   'End of region with all fibers',
                   ]

    if ref not in [0, 1]:
        raise ValueError("ref must be 0 or 1")

    off = (ref + 1) % 2

    if wcs is None:
        wcs = astropy.wcs.WCS(header)

    for key, com in zip(PIXLIM_KEYS, pixlim_coms):
        header[key] = (pixlims[key] + off, com)

    r1 = numpy.array([pixlims[key] for key in PIXLIM_KEYS])
    r2 = numpy.zeros_like(r1)
    lm = numpy.array([r1, r2])
    # Values are 0-based
    wavelen_ = wcs.all_pix2world(lm.T, ref)
    if wcs.wcs.cunit[0] == u.dimensionless_unscaled:
        # CUNIT is empty, assume Angstroms
        wavelen = wavelen_[:, 0] * u.AA
    else:
        wavelen = wavelen_[:, 0] * wcs.wcs.cunit[0]

    for idx, (key, com) in enumerate(zip(WAVLIM_KEYS, pixlim_coms)):

dapdb.Template,
                                    dapdb.Structure.template_kin_pk == dapdb.Template.pk).filter(
                                        dapdb.BinType.name == self.bintype.name,
                                        dapdb.Template.name == self.template.name).use_cache(self.cache_region).all()

                if len(db_modelcube) > 1:
                    raise MarvinError('more than one ModelCube found for '
                                      'this combination of parameters.')

                elif len(db_modelcube) == 0:
                    raise MarvinError('no ModelCube found for this combination of parameters.')

                self.data = db_modelcube[0]

            self.header = self.data.file.primary_header
            self.wcs = WCS(self.data.file.cube.wcs.makeHeader())
            self._wavelength = np.array(self.data.file.cube.wavelength.wavelength, dtype=np.float)
            self._redcorr = np.array(self.data.redcorr[0].value, dtype=np.float)
            self._shape = self.data.file.cube.shape.shape

            self.plateifu = str(self.header['PLATEIFU'].strip())
            self.mangaid = str(self.header['MANGAID'].strip())

def make_source_catalog():
    """Make source catalog from images.

    TODO: use other images to do measurements for the sources,
    e.g. excess, npred, flux.
    """
    significance_threshold = 7

    hdu = fits.open(TS_IMAGES)['sqrt_ts']
    header = fits.getheader(REF_IMAGE)
    wcs = WCS(header)

    print('Running find_peaks ...')
    table = find_peaks(data=hdu.data, threshold=significance_threshold, wcs=wcs)
    print('Number of sources detected: {}'.format(len(table)))

    # Add some useful columns
    icrs = SkyCoord(table['icrs_ra_peak'], table['icrs_dec_peak'], unit='deg')
    galactic = icrs.galactic
    table['Source_Name'] = coordinate_iau_format(icrs, ra_digits=5, prefix='J')
    table['GLON'] = galactic.l.deg
    table['GLAT'] = galactic.b.deg

    # table.show_in_browser(jsviewer=True)
    print('Writing {}'.format(SOURCE_CATALOG))
    table.write(SOURCE_CATALOG, overwrite=True)

if self.pos == None and self.tic != None:
            id, pos, tmag = data_products.tic_pos_by_ID(self)
            self.pos = pos

        in_file=[None]
        # Searches through rows of the table
        for i in range(len(t)):
            data=[]
            # Creates a list of the data in a row
            for j in range(146):
                data.append(t[i][j])
            d = dict(zip(t.colnames[0:146], data))
            hdr = fits.Header(cards=d)

            xy = WCS(hdr).all_world2pix(self.pos[0], self.pos[1], 1, quiet=True)
            x_cen, y_cen, l, w = t['POST_CEN_X'][i], t['POST_CEN_Y'][i], t['POST_SIZE1'][i]/2., t['POST_SIZE2'][i]/2.
            # Checks to see if xy coordinates of source falls within postcard
            if (xy[0] >= x_cen-l) & (xy[0] <= x_cen+l) & (xy[1] >= y_cen-w) & (xy[1] <= y_cen+w):
                if in_file[0]==None:
                    in_file[0]=i
                else:
                    in_file.append(i)
                # If more than one postcard is found for a single source, choose the postcard where the
                # source is closer to the center
                if len(in_file) > 1:
                    dist1 = np.sqrt( (xy[0]-t['POST_CENX'][in_files[0]])**2 + (xy[1]-t['POST_CENY'][in_files[0]])**2  )
                    dist2 = np.sqrt( (xy[0]-t['POST_CENX'][in_files[1]])**2 + (xy[1]-t['POST_CENY'][in_files[1]])**2  )
                    if dist1 >= dist2:
                        in_file[0]=in_file[1]
                    else:
                        in_file[0]=in_file[0]

xmin = x1 if x1 < x2 else x2
        xmax = x2 if x2 > x1 else x1

        ymin = y1 if y1 < y2 else y2
        ymax = y2 if y2 > y1 else y1

        # create a new WCS with the same transform but new center coordinates
        whdr = w.to_header()
        whdr['CRPIX1'] = (xmax - xmin)/2
        whdr['CRPIX2'] = (ymax - ymin)/2
        whdr['CRVAL1'] = cntra
        whdr['CRVAL2'] = cntdecl
        whdr['NAXIS1'] = xmax - xmin
        whdr['NAXIS2'] = ymax - ymin
        w = WCS(whdr)

    else:
        xmin, xmax, ymin, ymax = 0, hdr['NAXIS2'], 0, hdr['NAXIS1']

    # add the axes with the WCS projection
    # this should automatically handle subimages because we fix the WCS
    # appropriately above for these
    fig.add_subplot(111,projection=w)

    if scale is not None and stretch is not None:

        norm = ImageNormalize(img,
                              interval=scale,
                              stretch=stretch)

        plt.imshow(img[ymin:ymax,xmin:xmax],