Python Angle.wrap_at方法代码示例

# 需要导入模块: from astropy.coordinates import Angle [as 别名]
# 或者: from astropy.coordinates.Angle import wrap_at [as 别名]
def plot_vhs_tiles(table=None,ra=None, dec=None,
    wrap_ra24hr=False, PA=90.0, overplot=True, savefig=True):
    """

    table is astropy table
    ra, dec in degrees

    """

    plot_polygon=True

    if wrap_ra24hr:
        angle=Angle(xdata * u.deg)
        angle.wrap_at('180d', inplace=True)
        xdata=angle.degree

    # assumes filename is in table.meta
    print('plotting: ', table.meta['filename'])
    filename = table.meta['filename']

    if not plot_polygon:
      plt.plot(xdata, ydata, 's', ms=7.0,
       color='yellow', markeredgecolor='yellow',
       alpha=0.1,
       label='OB Progress (submitted)\n' + filename)

    if plot_polygon:
      plt.plot(xdata, ydata, '.', ms=7.0,
       color='yellow', markeredgecolor='yellow',
       alpha=0.1,
       label='OB Progress (submitted):' + filename)

      # filled polygon
      ra=xdata
      dec=ydata
      for i in range(len(ra)):
        if i==0 or i == len(ra): print('i: ', i)
        ra_poly, dec_poly=plt_tile.mk_polytile(ra_cen=ra[i], dec_cen=dec[i],
         coverage='twice', PA=PA)

        xypolygon=np.column_stack((ra_poly, dec_poly))
        if i==0 or i == len(ra):
          print('xypolygon.shape: ',   xypolygon.shape)
          print('xypolygon: ',xypolygon)

        polygon = Polygon(xypolygon, True, color='green', alpha=0.1)
        plt.gca().add_patch(polygon)

        #print(ob_table.meta)

      plt.suptitle(dqcfile_vhs)
      plt.legend(fontsize='small')

      figname='vhs_des_check_progress_vhs_obprogress_'+ datestamp + '.png'
      print('Saving: '+figname)
      plt.savefig(plotdir+figname)
      plt.suptitle('')

# 需要导入模块: from astropy.coordinates import Angle [as 别名]
# 或者: from astropy.coordinates.Angle import wrap_at [as 别名]
def test_3rd_body_Curtis(test_params):
    # based on example 12.11 from Howard Curtis
    body = test_params['body']
    with solar_system_ephemeris.set('builtin'):
        j_date = 2454283.0 * u.day
        tof = (test_params['tof']).to(u.s).value
        body_r = build_ephem_interpolant(body, test_params['period'], (j_date, j_date + test_params['tof']), rtol=1e-2)

        epoch = Time(j_date, format='jd', scale='tdb')
        initial = Orbit.from_classical(Earth, *test_params['orbit'], epoch=epoch)
        r, v = cowell(initial, np.linspace(0, tof, 400), rtol=1e-10, ad=third_body,
                      k_third=body.k.to(u.km**3 / u.s**2).value, third_body=body_r)

        incs, raans, argps = [], [], []
        for ri, vi in zip(r, v):
            angles = Angle(rv2coe(Earth.k.to(u.km**3 / u.s**2).value, ri, vi)[2:5] * u.rad)  # inc, raan, argp
            angles = angles.wrap_at(180 * u.deg)
            incs.append(angles[0].value)
            raans.append(angles[1].value)
            argps.append(angles[2].value)

        # averaging over 5 last values in the way Curtis does
        inc_f, raan_f, argp_f = np.mean(incs[-5:]), np.mean(raans[-5:]), np.mean(argps[-5:])

        assert_quantity_allclose([(raan_f * u.rad).to(u.deg) - test_params['orbit'][3],
                                  (inc_f * u.rad).to(u.deg) - test_params['orbit'][2],
                                  (argp_f * u.rad).to(u.deg) - test_params['orbit'][4]],
                                 [test_params['raan'], test_params['inc'], test_params['argp']],
                                 rtol=1e-1)

# 需要导入模块: from astropy.coordinates import Angle [as 别名]
# 或者: from astropy.coordinates.Angle import wrap_at [as 别名]
def sky_time(coord, time, rise_set=False, limalt=0*u.deg, site=EarthLocation(0.0, 0.0, 0.0), fuse=TimeDelta(0, format='sec', scale='tai')):
    """
    """
    if type(limalt) != u.quantity.Quantity:
        limalt = limalt*u.deg
    if time.isscalar == True:
        time = Time([time.iso], format='iso', scale='utc')
    coord = coord_pack(coord)
    timeut = time - fuse
    if len(time.shape) == 1:
        timeut = Time([[i] for i in timeut.jd], format='jd', scale='utc')
    timeut.delta_ut1_utc = 0
    timeut.location = site
    ra, ts = mesh_coord(coord, timeut[:,0])
    dif_h_sid = Angle(ra-ts)
    dif_h_sid.wrap_at('180d', inplace=True)
    dif_h_sol = dif_h_sid * (23.0 + 56.0/60.0 + 4.0916/3600.0) / 24.0
    dif = TimeDelta(dif_h_sol.hour*u.h, scale='tai')
    culminacao = timeut + dif
    culminacao.delta_ut1_utc = 0
    culminacao.location = site
    if (site.latitude > 0*u.deg):
        alwaysup = np.where(coord.dec >= 90*u.deg - site.latitude + limalt)
        neverup = np.where(coord.dec <= -(90*u.deg - site.latitude - limalt))
    else:
        alwaysup = np.where(coord.dec <= -(90*u.deg + site.latitude + limalt))
        neverup = np.where(coord.dec >= 90*u.deg + site.latitude - limalt)
    if rise_set == True:
        hangle_lim = np.arccos((np.cos(90.0*u.deg-limalt) - np.sin(coord.dec)*np.sin(site.latitude)) / (np.cos(coord.dec)*np.cos(site.latitude)))
        tsg_lim = Angle(ra + hangle_lim)
        dtsg_lim = tsg_lim - culminacao.sidereal_time('mean')
        dtsg_lim.wrap_at(360 * u.deg, inplace=True)
        dtsg_lim_sol = dtsg_lim * (23.0 + 56.0/60.0 + 4.0916/3600.0) / 24.0
        a = np.where(np.isnan(dtsg_lim_sol))
        dtsg_lim_sol[a] = Angle([48.0]*len(a[0])*u.hour)
        dtsg_np = TimeDelta((dtsg_lim_sol.hour*u.h))
        sunrise = culminacao - dtsg_np
        sunset = culminacao + dtsg_np
        culminacao = culminacao + fuse
        sunrise = sunrise + fuse
        sunset = sunset + fuse
        return culminacao, sunrise, sunset, alwaysup, neverup
    culminacao = culminacao + fuse
    return culminacao, alwaysup, neverup

# 需要导入模块: from astropy.coordinates import Angle [as 别名]
# 或者: from astropy.coordinates.Angle import wrap_at [as 别名]
def test_3rd_body_Curtis(test_params):
    # based on example 12.11 from Howard Curtis
    body = test_params["body"]
    with solar_system_ephemeris.set("builtin"):
        j_date = 2454283.0 * u.day
        tof = (test_params["tof"]).to(u.s).value
        body_r = build_ephem_interpolant(
            body,
            test_params["period"],
            (j_date, j_date + test_params["tof"]),
            rtol=1e-2,
        )

        epoch = Time(j_date, format="jd", scale="tdb")
        initial = Orbit.from_classical(Earth, *test_params["orbit"], epoch=epoch)
        rr, vv = cowell(
            Earth.k,
            initial.r,
            initial.v,
            np.linspace(0, tof, 400) * u.s,
            rtol=1e-10,
            ad=third_body,
            k_third=body.k.to(u.km ** 3 / u.s ** 2).value,
            third_body=body_r,
        )

        incs, raans, argps = [], [], []
        for ri, vi in zip(rr.to(u.km).value, vv.to(u.km / u.s).value):
            angles = Angle(
                rv2coe(Earth.k.to(u.km ** 3 / u.s ** 2).value, ri, vi)[2:5] * u.rad
            )  # inc, raan, argp
            angles = angles.wrap_at(180 * u.deg)
            incs.append(angles[0].value)
            raans.append(angles[1].value)
            argps.append(angles[2].value)

        # averaging over 5 last values in the way Curtis does
        inc_f, raan_f, argp_f = (
            np.mean(incs[-5:]),
            np.mean(raans[-5:]),
            np.mean(argps[-5:]),
        )

        assert_quantity_allclose(
            [
                (raan_f * u.rad).to(u.deg) - test_params["orbit"][3],
                (inc_f * u.rad).to(u.deg) - test_params["orbit"][2],
                (argp_f * u.rad).to(u.deg) - test_params["orbit"][4],
            ],
            [test_params["raan"], test_params["inc"], test_params["argp"]],
            rtol=1e-1,
        )

# 需要导入模块: from astropy.coordinates import Angle [as 别名]
# 或者: from astropy.coordinates.Angle import wrap_at [as 别名]
 def fit_arrays(sci, wht, seg, psf, id=None, platescale=0.06, exptime=0, path='/tmp/', galfit_exec='galfit', gaussian_guess=False, components=[GalfitSersic()], recenter=True, psf_sample=1):
     
     rms = 1/np.sqrt(wht)#*exptime
     if exptime > 0:
         rms = np.sqrt((rms*exptime)**2+sci*exptime*(sci > 0))/exptime
     
     rms[wht == 0] = 1e30
     
     if id is not None:
         mask = ((seg > 0) & (seg != id)) | (wht == 0)
     else:
         mask = wht == 0
     
     sh = sci.shape[0]
     
     fp = open(path+'galfit.feedme','w')
     
     fp.write(GALFIT_IMAGES.format(input=path+'gf_sci.fits', output=path+'gf_out.fits', sigma=path+'gf_rms.fits', psf=path+'gf_psf.fits', mask=path+'gf_mask.fits', xmax=sh, ymax=sh, sh=sh, ps=platescale, psf_sample=psf_sample))
     
     if gaussian_guess:
         fit, q, theta = fit_gauss(sci)
     
     from astropy.coordinates import Angle
     import astropy.units as u
         
     for comp in components:
         if recenter:
             comp.set(pos=[sh/2., sh/2.])
         
         if gaussian_guess:
             comp.set(q=q, pa=Angle.wrap_at(theta*u.rad, 360*u.deg).to(u.deg).value)
             
         fp.write(str(comp))
         
     fp.close()
     
     pyfits.writeto(path+'gf_sci.fits', data=sci, overwrite=True)
     pyfits.writeto(path+'gf_rms.fits', data=rms, overwrite=True)
     pyfits.writeto(path+'gf_mask.fits', data=mask*1, overwrite=True)        
     pyfits.writeto(path+'gf_psf.fits', data=psf, overwrite=True)
     
     for ext in ['out', 'model']:
         if os.path.exists(path+'gf_{0}.fits'.format(ext)):
             os.remove(path+'gf_{0}.fits'.format(ext))
         
     os.system('{0} {1}/galfit.feedme'.format(galfit_exec, path))

# 需要导入模块: from astropy.coordinates import Angle [as 别名]
# 或者: from astropy.coordinates.Angle import wrap_at [as 别名]
 def get_dispersion_PA(self):
     from astropy.coordinates import Angle
     import astropy.units as u
                 
     ### extra tilt of the 1st order grism spectra
     x0 =  self.conf.conf['BEAMA']
     dy_trace, lam_trace = self.conf.get_beam_trace(x=507, y=507, dx=x0, beam='A')
     extra = np.arctan2(dy_trace[1]-dy_trace[0], x0[1]-x0[0])/np.pi*180
     
     # h = self.im_header
     # pa =  Angle((-np.arctan2(h['CD2_2'], h['CD2_1'])/np.pi*180-extra)*u.deg)
     
     ### Distorted WCS
     crpix = self.flt_wcs.wcs.crpix
     xref = [crpix[0], crpix[0]+1]
     yref = [crpix[1], crpix[1]]
     r, d = self.all_pix2world(xref, yref)
     pa =  Angle((extra+np.arctan2(np.diff(r), np.diff(d))[0]/np.pi*180)*u.deg)
     
     self.dispersion_PA = pa.wrap_at(360*u.deg).value

# 需要导入模块: from astropy.coordinates import Angle [as 别名]
# 或者: from astropy.coordinates.Angle import wrap_at [as 别名]
def calcfaixa(vel, data, star, dist, ca, pa, tamanho, step, erro=None, ring=None, atm=None):
    vec = np.arange(0, int(8000/(np.absolute(vel.value))), step)
    g = np.sort(np.concatenate((vec,-vec[1:]), axis=0))
    latlon = {'clat':{'lon':[], 'lat':[], 'lab': [], 'x': [], 'y': [], 'labx': []}, 'lats': {'lon':[], 'lat':[], 'lon2':[], 'lat2':[], 'x': [], 'y': [], 'x2':[], 'y2':[]}}
    if not erro == None:
        latlon['erro'] = {'lon': [], 'lat': [], 'lon2':[], 'lat2':[]}
        err = erro*u.mas
        errd = (dist.to(u.km)*err.to(u.rad)).value*u.km
    if not ring == None:
        latlon['ring'] = {'lon': [], 'lat': [], 'lon2':[], 'lat2':[]}
    if not atm == None:
        latlon['atm'] = {'lon': [], 'lat': [], 'lon2':[], 'lat2':[]}
    pa = Angle(pa)
    pa.wrap_at('180d', inplace=True)
    if pa > 90*u.deg:
        paplus = pa - 180*u.deg
    elif pa < -90*u.deg:
        paplus = pa + 180*u.deg
    else:
        paplus = pa
    for delt in g:
        deltatime = delt*u.s
        datas1 = data + TimeDelta(deltatime)
        datas1.delta_ut1_utc = 0
        lon = star.ra - datas1.sidereal_time('mean', 'greenwich')
        m = Basemap(projection='ortho',lat_0=star.dec.value,lon_0=lon.value,resolution=None)
        a, b = m(lon.value, star.dec.value)
        a = a*u.m
        b = b*u.m
        dista = (dist.to(u.km)*ca.to(u.rad)).value*u.km
        ax = a + dista*np.sin(pa) + (deltatime*vel)*np.cos(paplus)
        by = b + dista*np.cos(pa) - (deltatime*vel)*np.sin(paplus)
        ax2 = ax - tamanho/2*np.sin(paplus)
        by2 = by - tamanho/2*np.cos(paplus)
        ax3 = ax + tamanho/2*np.sin(paplus)
        by3 = by + tamanho/2*np.cos(paplus)
        clon1, clat1 = m(ax.value, by.value, inverse=True)
        if delt == 0:
            latlon['clat']['cxy'] = [ax.value, by.value]
        if clon1 < 1e+30:
            latlon['clat']['lon'].append(clon1)
            latlon['clat']['lat'].append(clat1)
            latlon['clat']['lab'].append(datas1.iso)
        else:
            latlon['clat']['x'].append(ax.value)
            latlon['clat']['y'].append(by.value)
            latlon['clat']['labx'].append(datas1.iso)
        lon1, lat1 = m(ax2.value, by2.value, inverse=True)
        if lon1 < 1e+30:
            latlon['lats']['lon'].append(lon1) 
            latlon['lats']['lat'].append(lat1)
        else:
            latlon['lats']['x'].append(ax2.value) 
            latlon['lats']['y'].append(by2.value)
        lon2, lat2 = m(ax3.value, by3.value, inverse=True)
        if lon2 < 1e+30:
            latlon['lats']['lon2'].append(lon2) 
            latlon['lats']['lat2'].append(lat2)
        else:
            latlon['lats']['x2'].append(ax3.value) 
            latlon['lats']['y2'].append(by3.value)
        if not erro == None:
            ax2 = ax - errd*np.sin(paplus)
            by2 = by - errd*np.cos(paplus)
            ax3 = ax + errd*np.sin(paplus)
            by3 = by + errd*np.cos(paplus)
            lon1, lat1 = m(ax2.value, by2.value, inverse=True)
            if lon1 < 1e+30:
                latlon['erro']['lon'].append(lon1) 
                latlon['erro']['lat'].append(lat1)
            lon2, lat2 = m(ax3.value, by3.value, inverse=True)
            if lon2 < 1e+30:
                latlon['erro']['lon2'].append(lon2) 
                latlon['erro']['lat2'].append(lat2)
        if not ring == None:
            rng = ring*u.km
            ax2 = ax - rng*np.sin(paplus)
            by2 = by - rng*np.cos(paplus)
            ax3 = ax + rng*np.sin(paplus)
            by3 = by + rng*np.cos(paplus)
            lon1, lat1 = m(ax2.value, by2.value, inverse=True)
            if lon1 < 1e+30:
                latlon['ring']['lon'].append(lon1) 
                latlon['ring']['lat'].append(lat1)
            lon2, lat2 = m(ax3.value, by3.value, inverse=True)
            if lon2 < 1e+30:
                latlon['ring']['lon2'].append(lon2) 
                latlon['ring']['lat2'].append(lat2)
        if not atm == None:
            atmo = atm*u.km
            ax2 = ax - atmo*np.sin(paplus)
            by2 = by - atmo*np.cos(paplus)
            ax3 = ax + atmo*np.sin(paplus)
            by3 = by + atmo*np.cos(paplus)
            lon1, lat1 = m(ax2.value, by2.value, inverse=True)
            if lon1 < 1e+30:
                latlon['atm']['lon'].append(lon1) 
                latlon['atm']['lat'].append(lat1)
            lon2, lat2 = m(ax3.value, by3.value, inverse=True)
            if lon2 < 1e+30:
#.........这里部分代码省略.........

# 需要导入模块: from astropy.coordinates import Angle [as 别名]
# 或者: from astropy.coordinates.Angle import wrap_at [as 别名]

#.........这里部分代码省略.........
        plt.savefig(SPLAT_PATH+DOCS_FOLDER+'_images/spt_optical_distribution_{}.png'.format(fname))
        plt.clf()
    # NIR type
        sptarr = numpy.array(nir_spts)
        plt.figure(figsize=(14,6))
        n, bins, patches = plt.hist(sptarr[numpy.where(numpy.logical_and(sptarr >= sptrng[0],sptarr < 20))], bins=len(range(sptrng[0],sptrng[1])), log=True, range=sptrng, facecolor='green', alpha=0.75)
        n, bins, patches = plt.hist(sptarr[numpy.where(numpy.logical_and(sptarr >= 20,sptarr < 30))], bins=len(range(sptrng[0],sptrng[1])), log=True, range=sptrng, facecolor='red', alpha=0.75)
        n, bins, patches = plt.hist(sptarr[numpy.where(numpy.logical_and(sptarr >= 30,sptarr < sptrng[1]))], bins=len(range(sptrng[0],sptrng[1])), log=True, range=sptrng, facecolor='b', alpha=0.75)
        plt.xticks(xticks,labels)
        plt.xlabel('NIR Spectral Type')
        plt.ylabel('log10 Number')
        plt.xlim([sptrng[0]-0.5,sptrng[1]+0.5])
        plt.legend(['M dwarfs ({} sources)'.format(len(sptarr[numpy.where(numpy.logical_and(sptarr >= sptrng[0],sptarr < 20))])),'L dwarfs ({} sources)'.format(len(sptarr[numpy.where(numpy.logical_and(sptarr >= 20,sptarr < 30))])),'T dwarfs ({} sources)'.format(len(sptarr[numpy.where(numpy.logical_and(sptarr >= 30,sptarr < sptrng[1]))]))])
        plt.savefig(SPLAT_PATH+DOCS_FOLDER+'_images/spt_nir_distribution_{}.png'.format(fname))
        plt.clf()
    # Adopted type
        sptarr = numpy.array(spts)
        plt.figure(figsize=(14,6))
        n, bins, patches = plt.hist(sptarr[numpy.where(numpy.logical_and(sptarr >= sptrng[0],sptarr < 20))], bins=len(range(sptrng[0],sptrng[1])), log=True, range=sptrng, facecolor='green', alpha=0.75)
        n, bins, patches = plt.hist(sptarr[numpy.where(numpy.logical_and(sptarr >= 20,sptarr < 30))], bins=len(range(sptrng[0],sptrng[1])), log=True, range=sptrng, facecolor='red', alpha=0.75)
        n, bins, patches = plt.hist(sptarr[numpy.where(numpy.logical_and(sptarr >= 30,sptarr < sptrng[1]))], bins=len(range(sptrng[0],sptrng[1])), log=True, range=sptrng, facecolor='b', alpha=0.75)
        plt.xticks(xticks,labels)
        plt.xlabel('Adopted Spectral Type')
        plt.ylabel('log10 Number')
        plt.xlim([sptrng[0]-0.5,sptrng[1]+0.5])
        plt.legend(['M dwarfs ({} sources)'.format(len(sptarr[numpy.where(numpy.logical_and(sptarr >= sptrng[0],sptarr < 20))])),'L dwarfs ({} sources)'.format(len(sptarr[numpy.where(numpy.logical_and(sptarr >= 20,sptarr < 30))])),'T dwarfs ({} sources)'.format(len(sptarr[numpy.where(numpy.logical_and(sptarr >= 30,sptarr < sptrng[1]))]))])
        plt.savefig(SPLAT_PATH+DOCS_FOLDER+'_images/spt_adopted_distribution_{}.png'.format(fname))
        plt.clf()

# histogram of S/N

# map sources on sky
    raref = Angle(numpy.linspace(0,359.,360)*u.degree)
    raref.wrap_at(180.*u.degree)
    ra = Angle(sall['RA'].filled(numpy.nan)*u.degree)
    ra = ra.wrap_at(180*u.degree)
    dec = Angle(sall['DEC'].filled(numpy.nan)*u.degree)
    rap = Angle(s['RA'].filled(numpy.nan)*u.degree)
    rap = rap.wrap_at(180*u.degree)
    decp = Angle(s['DEC'].filled(numpy.nan)*u.degree)
    fig = plt.figure(figsize=(8,6))
    ax = fig.add_subplot(111, projection="mollweide")
    p1 = ax.scatter(ra.radian, dec.radian,color='r',alpha=0.5,s=10)
    p2 = ax.scatter(rap.radian, decp.radian,color='k',alpha=0.5, s=5)
#    ur = ax.plot(raref.radian,Angle([67.]*len(raref)*u.degree).radian,'k--')
#    ur = ax.plot(raref.radian,Angle([-50.]*len(raref)*u.degree).radian,'k--')
    ax.set_xticklabels(['14h','16h','18h','20h','22h','0h','2h','4h','6h','8h','10h'])
    ax.grid(True)
#    ef = matplotlib.patheffects.withStroke(foreground="w", linewidth=3)
#    axis = ax.axis['lat=0']
#    axis.major_ticklabels.set_path_effects([ef])
 #   axis.label.set_path_effects([ef])
    plt.legend([p1,p2],['All Sources ({})'.format(len(sall)),'Published Sources ({})'.format(len(s))],bbox_to_anchor=(1, 1),bbox_transform=plt.gcf().transFigure)
    fig.savefig(SPLAT_PATH+DOCS_FOLDER+'_images/map_all.png')
    fig.clf()

# map sources on based on spectral type
    fig = plt.figure(figsize=(8,6))
    ax = fig.add_subplot(111, projection="mollweide")
    sm = splat.searchLibrary(spt_range=[10,19.9],spt_type='SPEX')
    sm = sm[sm['OBJECT_TYPE'] == 'VLM']
    sl = splat.searchLibrary(spt_range=[20,29.9],spt_type='SPEX')
    sl = sl[sl['OBJECT_TYPE'] == 'VLM']
    st = splat.searchLibrary(spt_range=[30,39.9],spt_type='SPEX')
    st = st[st['OBJECT_TYPE'] == 'VLM']
    ra = Angle(sm['RA'].filled(numpy.nan)*u.degree)

# 需要导入模块: from astropy.coordinates import Angle [as 别名]
# 或者: from astropy.coordinates.Angle import wrap_at [as 别名]
 def fit_object(self, id=449, radec=(None, None), size=40, components=[GalfitSersic()], recenter=True, get_mosaic=True, gaussian_guess=False, get_extended=True):
     """
     Fit an object
     """
     if id is not None:
         rd = self.cat['X_WORLD', 'Y_WORLD'][self.cat['NUMBER'] == id][0]
         radec = tuple(rd)
     
     xy = self.wcs.all_world2pix([radec[0]], [radec[1]], 0)
     xy = np.array(xy).flatten()
     xp = np.cast[int](np.round(xy))
     
     slx, sly, wcs_slice = self.DPSF.get_driz_cutout(ra=radec[0], dec=radec[1], get_cutout=False, size=size)
     #drz_cutout = self.get_driz_cutout(ra=ra, dec=dec, get_cutout=True)
     h = self.sci[0].header
     psf = self.DPSF.get_psf(ra=radec[0], dec=radec[1], filter=self.filter.upper(), wcs_slice=wcs_slice, pixfrac=h['D001PIXF'], kernel=h['D001KERN'], get_extended=get_extended)
     
     exptime = h['EXPTIME']
     sci = self.sci[0].data[sly, slx]#*exptime
     wht = self.wht[0].data[sly, slx]
     rms = 1/np.sqrt(wht)#*exptime
     rms = np.sqrt((rms*exptime)**2+sci*exptime*(sci > 0))/exptime
     
     rms[wht == 0] = 1e30
     
     seg = self.seg[0].data[sly, slx]
     if id is not None:
         mask = ((seg > 0) & (seg != id)) | (wht == 0)
     else:
         mask = wht == 0
     
     sh = sci.shape[0]
     path = '/tmp/'
     
     psf_file = '{0}-{1}_psf_{2:05d}.fits'.format(self.root, self.filter, id)
     
     fp = open(path+'galfit.feedme','w')
     
     fp.write(GALFIT_IMAGES.format(input=path+'gf_sci.fits', output=path+'gf_out.fits', sigma=path+'gf_rms.fits', psf=psf_file, mask=path+'gf_mask.fits', xmax=sh, ymax=sh, sh=sh, ps=self.DPSF.driz_pscale, psf_sample=1))
     
     if gaussian_guess:
         fit, q, theta = fit_gauss(sci)
     
     from astropy.coordinates import Angle
     import astropy.units as u
         
     for comp in components:
         if recenter:
             comp.set(pos=[sh/2., sh/2.])
         
         if gaussian_guess:
             comp.set(q=q, pa=Angle.wrap_at(theta*u.rad, 360*u.deg).to(u.deg).value)
             
         fp.write(str(comp))
         
     fp.close()
     
     pyfits.writeto(path+'gf_sci.fits', data=sci, overwrite=True)
     pyfits.writeto(path+'gf_rms.fits', data=rms, overwrite=True)
     pyfits.writeto(path+'gf_mask.fits', data=mask*1, overwrite=True)
     #pyfits.writeto(path+'gf_psf.fits', data=psf[1].data, overwrite=True)
     
     pyfits.writeto(psf_file, data=psf[1].data, overwrite=True)
     
     for ext in ['out', 'model']:
         if os.path.exists(path+'gf_{0}.fits'.format(ext)):
             os.remove(path+'gf_{0}.fits'.format(ext))
         
     os.system('{0} {1}/galfit.feedme'.format(self.galfit_exec, path))
     
     #out = pyfits.open('/tmp/gf_out.fits')
     
     # Get in DRZ frame
     gf_file = glob.glob('galfit.[0-9]*')[-1]
     lines = open(gf_file).readlines()
     for il, line in enumerate(lines):
         if line.startswith('A)'):
             lines[il] = 'A) {0}     # Input data image (FITS file)\n'.format(self.sci.filename())
         if line.startswith('B)'):
             lines[il] = 'B) {0}-{1}_galfit_{2:05d}.fits      # Output data image block\n'.format(self.root, self.filter, id)
         if line.startswith('P) 0'):
             lines[il] = "P) 1                   # Choose: 0=optimize, 1=model, 2=imgblock, 3=subcomps\n"
         
         if line.startswith('H)'):
             out_sh = self.sci[0].data.shape
             lines[il] = "H) 1    {0}   1    {1}   # Image region to fit (xmin xmax ymin ymax)\n".format(out_sh[1], out_sh[0])
             
         if line.startswith(' 1)'):
             xy = np.cast[float](line.split()[1:3])
             lines[il] = ' 1) {0}  {1}  1 1  #  Position x, y\n'.format(xy[0]+slx.start, xy[1]+sly.start)
             
     fp = open('galfit.{0}.{1:05d}'.format(self.filter, id),'w')
     fp.writelines(lines)
     fp.close()
     
     os.system('{0} galfit.{1}.{2:05d}'.format(self.galfit_exec, self.filter, id))
     
     model = pyfits.open('{0}-{1}_galfit_{2:05d}.fits'.format(self.root, self.filter, id), mode='update')
     model[0].data /= self.sci[0].header['EXPTIME']
     model.flush()
#.........这里部分代码省略.........

# 需要导入模块: from astropy.coordinates import Angle [as 别名]
# 或者: from astropy.coordinates.Angle import wrap_at [as 别名]
    ydata=despolygon['dec']

    plt.plot(xdata, ydata, label='DES Round13-poly')
    plotid(progname=True)

    figname='vhs_des_check_sadt_despolygon_' + datestamp + '.png'
    print('Saving: '+figname)
    plt.savefig(plotdir+figname)
    plt.suptitle('')


    if viking:
        xdata=viking['ra']
        ydata=viking['dec']
        angle=Angle(xdata * u.deg)
        angle.wrap_at('180d', inplace=True)
        xdata=angle.degree

        plt.plot(xdata, ydata, 'sr',
            ms=7.0, markeredgecolor='r', alpha=0.2, label='VIKING')
        plt.suptitle(dqcfile_viking)

        figname='vhs_des_check_viking_1_' + datestamp + '.png'
        print('Saving: '+figname)
        plt.savefig(plotdir+figname)
        plt.suptitle('')

    if vhs_progress:
        print('Plot VHS obprogress data')
        plot_vhs_obprogress(vhs_obprogress=vhs_obprogress)
        raw_input("Enter any key to continue: ")

# 需要导入模块: from astropy.coordinates import Angle [as 别名]
# 或者: from astropy.coordinates.Angle import wrap_at [as 别名]
def pspec(psd2, nbins=None, return_stddev=False, binsize=1.0,
          logspacing=True, max_bin=None, min_bin=None, return_freqs=True,
          theta_0=None, delta_theta=None, boot_iter=None):
    '''
    Calculate the radial profile using scipy.stats.binned_statistic.

    Parameters
    ----------
    psd2 : np.ndarray
        2D Spectral power density.
    nbins : int, optional
        Number of bins to use. If None, it is calculated based on the size
        of the given arrays.
    return_stddev : bool, optional
        Return the standard deviations in each bin.
    binsize : float, optional
        Size of bins to be used. If logspacing is enabled, this will increase
        the number of bins used by the inverse of the given binsize.
    logspacing : bool, optional
        Use logarithmically spaces bins.
    max_bin : float, optional
        Give the maximum value to bin to.
    min_bin : float, optional
        Give the minimum value to bin to.
    return_freqs : bool, optional
        Return spatial frequencies.
    theta_0 : `~astropy.units.Quantity`, optional
        The center angle of the azimuthal mask. Must have angular units.
    delta_theta : `~astropy.units.Quantity`, optional
        The width of the azimuthal mask. This must be given when
        a `theta_0` is given. Must have angular units.
    boot_iter : int, optional
        Number of bootstrap iterations for estimating the standard deviation
        in each bin. Require `return_stddev=True`.

    Returns
    -------
    bins_cents : np.ndarray
        Centre of the bins.
    ps1D : np.ndarray
        1D binned power spectrum.
    ps1D_stddev : np.ndarray
        Returned when return_stddev is enabled. Standard deviations
        within each of the bins.
    '''

    yy, xx = make_radial_arrays(psd2.shape)

    dists = np.sqrt(yy**2 + xx**2)
    if theta_0 is not None:

        if delta_theta is None:
            raise ValueError("Must give delta_theta.")

        theta_0 = theta_0.to(u.rad)
        delta_theta = delta_theta.to(u.rad)

        theta_limits = Angle([theta_0 - 0.5 * delta_theta,
                              theta_0 + 0.5 * delta_theta])

        # Define theta array
        thetas = Angle(np.arctan2(yy, xx) * u.rad)

        # Wrap around pi
        theta_limits = theta_limits.wrap_at(np.pi * u.rad)

    if nbins is None:
        nbins = int(np.round(dists.max() / binsize) + 1)

    if return_freqs:
        yy_freq, xx_freq = make_radial_freq_arrays(psd2.shape)

        freqs_dist = np.sqrt(yy_freq**2 + xx_freq**2)

        zero_freq_val = freqs_dist[np.nonzero(freqs_dist)].min() / 2.
        freqs_dist[freqs_dist == 0] = zero_freq_val

    if max_bin is None:
        if return_freqs:
            max_bin = 0.5
        else:
            max_bin = dists.max()

    if min_bin is None:
        if return_freqs:
            min_bin = 1.0 / min(psd2.shape)
        else:
            min_bin = 0.5

    if logspacing:
        bins = np.logspace(np.log10(min_bin), np.log10(max_bin), nbins + 1)
    else:
        bins = np.linspace(min_bin, max_bin, nbins + 1)

    if return_freqs:
        dist_arr = freqs_dist
    else:
        dist_arr = dists

    if theta_0 is not None:
#.........这里部分代码省略.........