我们从Python开源项目中，提取了以下17个代码示例，用于说明如何使用astropy.units.degree()。
def __init__(self,x,y, frame = FK5): ''' init function Parameters --------- x : float, first coordinate of the source f : float, second coordinate of the source frame : Astropy coordinate frame ICRS, Galactic, FK4, FK5 , see astropy for more information ''' self.frame = frame try : self.skycoord = SkyCoord(x, y, frame=frame) except: self.skycoord = SkyCoord(x*u.degree, y*u.degree, frame=frame) try : self.X = self.skycoord.ra self.Y = self.skycoord.dec except: self.X = self.skycoord.l self.Y = self.skycoord.b self.to_string = self.skycoord.to_string self.to_string = self.skycoord.to_string self.frame = self.skycoord.frame.name
def overlap_check(ra, dec, radius): """ A very rough implementation, not using any distance or coverage calculation for efficiency """ for i, re in enumerate(overlap_regions): if (ra >= re[0] and ra <= re[1] and dec >= re[2] and dec <= re[3] and radius <= re[4]): #raise Exception("Sorry the region '{0}' does not contain any GLEAM sources".format(overlap_names[i])) raise Exception(overlap_err) icrs_coord = SkyCoord(ra * u.degree, dec * u.degree, frame='icrs') gg = icrs_coord.galactic b = float(str(gg.to_string().split(',')[0]).split()[1]) if (b < 9 and b > -9): #raise Exception("Sorry Galactic Plane does not contain any GLEAM sources.") raise Exception(overlap_err)
def matching(objra,objdec,ras,decs,boxsize=1): """ objra and objdec are in degrees ras and decs are numpy arrays of all the catalog values (also in degrees) """ obj = coord.ICRSCoordinates(ra=objra, dec=objdec, unit=(u.degree, u.degree)) distances = [] for i in range(len(ras)): if abs(ras[i] - objra) > 0.1 or abs(decs[i] - objdec) > 0.1: distances.append(9999) else: catobj = coord.ICRSCoordinates(ra=ras[i], dec=decs[i], unit=(u.degree, u.degree)) distances.append(obj.separation(catobj).arcsecs) distances = np.array(distances) minindex = np.argmin(distances) mindist = np.amin(distances) b = distances < boxsize boxnum = len(distances[b]) return [minindex,mindist,boxnum]
def match_to_database(source_list, image_list, cat, Ast, imscale, thresh=8.): ''' RA, Decdetection catalogs and vizier (astroquery) match? vo conesearch? ''' matches = [] for sources, image in zip(source_list, image_list): ifuslot = sources[1] Ast.get_ifuslot_projection(ifuslot, imscale, image.shape[1]/2., image.shape[0]/2.) for source in sources[0]: ra, dec = Ast.tp_ifuslot.wcs.wcs_pix2world(source['xcentroid'], source['ycentroid'],1) c = SkyCoord(ra, dec, unit=(unit.degree, unit.degree)) idx, d2d, d3d = match_coordinates_sky(c, cat, nthneighbor=1) if d2d.arcsec[0] < thresh: dra = cat.ra[idx].deg - float(ra) ddec = cat.dec[idx].deg - float(dec) matches.append([int(ifuslot),int(idx), float(ra), float(dec), cat.ra[idx].deg, cat.dec[idx].deg, dra, ddec, d2d.arcsec[0]]) return matches
def xloop(i): c = SkyCoord(ra=df2['alpha_j2000'][i]*u.degree, dec=df2['delta_j2000'][i]*u.degree) #print i val = c.cartesian.x return float(val)
def zloop(i): c = SkyCoord(ra=df2['alpha_j2000'][i]*u.degree, dec=df2['delta_j2000'][i]*u.degree) #print i val = c.cartesian.z return float(val) #print val #print df2['mmm'][i]
def yloop(i): c = SkyCoord(ra=df2['alpha_j2000'][i]*u.degree, dec=df2['delta_j2000'][i]*u.degree) #print i val = c.cartesian.y return float(val)
def SetRADEC(self,RA,DEC): c_icrs = CH.CoordinatesHandler(RA * u.degree, DEC * u.degree,ICRS) self.config["target"]['ra'] = c_icrs.X.degree self.config["target"]['dec'] =c_icrs.Y.degree self.config["target"]['l'] =c_icrs.skycoord.galactic.l.value self.config["target"]['b'] =c_icrs.skycoord.galactic.b.value self._set_center()
def SetLB(self,L,B): c_gal = CH.CoordinatesHandler(L * u.degree, B * u.degree,Galactic) self.config["target"]['ra'] = c_gal.skycoord.icrs.ra.value self.config["target"]['dec'] = c_gal.skycoord.icrs.dec.value self.config["target"]['l'] =c_gal.X.degree self.config["target"]['b'] =c_gal.Y.degree self._set_center()
def parallatic_angle(ha, dec, geolat): ''' Parallactic angle of a source in degrees Parameters ---------- ha : array_like Hour angle, in hours dec : float Declination, in degrees geolat : float Observatory declination, in degrees Returns ------- pa : array_like Parallactic angle values ''' pa = -np.arctan2(-np.sin(ha), np.cos(dec) * np.tan(geolat) - np.sin(dec) * np.cos(ha)) if (dec >= geolat): pa[ha < 0] += 360*units.degree return np.degrees(pa)
def query_radecl(ra, decl, filtersystem='sloan_2mass', field_deg2=1.0, usebinaries=True, extinction_sigma=0.1, magnitude_limit=26.0, maglim_filtercol=4, trilegal_version=1.6, extraparams=None, forcefetch=False, cachedir='~/.astrobase/trilegal-cache', verbose=True, timeout=60.0, refresh=150.0, maxtimeout=700.0): ''' This runs the TRILEGAL query for decimal equatorial coordinates. ''' # convert the ra/decl to gl, gb radecl = SkyCoord(ra=ra*u.degree, dec=decl*u.degree) gl = radecl.galactic.l.degree gb = radecl.galactic.b.degree return query_galcoords(gl, gb, filtersystem=filtersystem, field_deg2=field_deg2, usebinaries=usebinaries, extinction_sigma=extinction_sigma, magnitude_limit=magnitude_limit, maglim_filtercol=maglim_filtercol, trilegal_version=trilegal_version, extraparams=extraparams, forcefetch=forcefetch, cachedir=cachedir, verbose=verbose, timeout=timeout, refresh=refresh, maxtimeout=maxtimeout)
def xmatch(cat1,cat2,maxdist=2, colRA1='RA',colDec1='DEC',epoch1=2000., colRA2='RA',colDec2='DEC',epoch2=2000., colpmRA2='pmra',colpmDec2='pmdec', swap=False): """ NAME: xmatch PURPOSE: cross-match two catalogs (incl. proper motion in cat2 if epochs are different) INPUT: cat1 - First catalog cat2 - Second catalog maxdist= (2) maximum distance in arcsec colRA1= ('RA') name of the tag in cat1 with the right ascension in degree in cat1 (assumed to be ICRS) colDec1= ('DEC') name of the tag in cat1 with the declination in degree in cat1 (assumed to be ICRS) epoch1= (2000.) epoch of the coordinates in cat1 colRA2= ('RA') name of the tag in cat2 with the right ascension in degree in cat2 (assumed to be ICRS) colDec2= ('DEC') name of the tag in cat2 with the declination in degree in cat2 (assumed to be ICRS) epoch2= (2000.) epoch of the coordinates in cat2 colpmRA2= ('pmra') name of the tag in cat2 with the proper motion in right ascension in degree in cat2 (assumed to be ICRS; includes cos(Dec)) [only used when epochs are different] colpmDec2= ('pmdec') name of the tag in cat2 with the proper motion in declination in degree in cat2 (assumed to be ICRS) [only used when epochs are different] swap= (False) if False, find closest matches in cat2 for each cat1 source, if False do the opposite (important when one of the catalogs has duplicates) OUTPUT: (index into cat1 of matching objects, index into cat2 of matching objects, angular separation between matching objects) HISTORY: 2016-09-12 - Written - Bovy (UofT) 2016-09-21 - Account for Gaia epoch 2015 - Bovy (UofT) """ if ('ref_epoch' in cat1.dtype.fields and numpy.fabs(epoch1-2015.) > 0.01)\ or ('ref_epoch' in cat2.dtype.fields and \ numpy.fabs(epoch2-2015.) > 0.01): warnings.warn("You appear to be using a Gaia catalog, but are not setting the epoch to 2015., which may lead to incorrect matches") depoch= epoch2-epoch1 if depoch != 0.: # Use proper motion to get both catalogs at the same time dra=cat2[colpmRA2]/numpy.cos(cat2[colDec2]/180.*numpy.pi)\ /3600000.*depoch ddec= cat2[colpmDec2]/3600000.*depoch else: dra= 0. ddec= 0. mc1= acoords.SkyCoord(cat1[colRA1],cat1[colDec1], unit=(u.degree, u.degree),frame='icrs') mc2= acoords.SkyCoord(cat2[colRA2]-dra,cat2[colDec2]-ddec, unit=(u.degree, u.degree),frame='icrs') if swap: idx,d2d,d3d = mc2.match_to_catalog_sky(mc1) m1= numpy.arange(len(cat2)) else: idx,d2d,d3d = mc1.match_to_catalog_sky(mc2) m1= numpy.arange(len(cat1)) mindx= d2d < maxdist*u.arcsec m1= m1[mindx] m2= idx[mindx] if swap: return (m2,m1,d2d[mindx]) else: return (m1,m2,d2d[mindx])
def compute_times(frames_info): ''' Compute the various timestamps associated to frames Parameters ---------- frames_info : dataframe The data frame with all the information on science frames ''' # get necessary values time_start = frames_info['DATE-OBS'].values time_end = frames_info['DET FRAM UTC'].values time_delta = (time_end - time_start) / frames_info['DET NDIT'].values.astype(np.int) DIT = np.array(frames_info['DET SEQ1 DIT'].values.astype(np.float)*1000, dtype='timedelta64[ms]') # calculate UTC time stamps idx = frames_info.index.get_level_values(1).values ts_start = time_start + time_delta * idx ts = time_start + time_delta * idx + DIT/2 ts_end = time_start + time_delta * idx + DIT # calculate mjd geolon = coord.Angle(frames_info['TEL GEOLON'].values[0], units.degree) geolat = coord.Angle(frames_info['TEL GEOLAT'].values[0], units.degree) geoelev = frames_info['TEL GEOELEV'].values[0] utc = Time(ts_start.astype(str), scale='utc', location=(geolon, geolat, geoelev)) mjd_start = utc.mjd utc = Time(ts.astype(str), scale='utc', location=(geolon, geolat, geoelev)) mjd = utc.mjd utc = Time(ts_end.astype(str), scale='utc', location=(geolon, geolat, geoelev)) mjd_end = utc.mjd # update frames_info frames_info['TIME START'] = ts_start frames_info['TIME'] = ts frames_info['TIME END'] = ts_end frames_info['MJD START'] = mjd_start frames_info['MJD'] = mjd frames_info['MJD END'] = mjd_end
def load_fits_file(self, msg): filename = msg if isinstance(filename, str) or isinstance(filename, unicode): if filename.lower().endswith('.fits') or filename.lower().endswith('.fits.gz'): if os.path.isfile(filename): # TODO: be more flexible about hdulist where image data is NOT just [0].data # TODO also, in case of extended fits files need to deal with additional header info # following try/except handles situation when autoloading files tries to autoload a file # before it's been fully written to disk. max_n_tries = 5 pause_time_between_tries_sec = 1. cur_try = 0 not_yet_successful = True while (cur_try < max_n_tries) and not_yet_successful: try: self.cur_fits_hdulist = self.load_hdulist_from_fitsfile(filename) self.load_numpy_array(self.cur_fits_hdulist[0].data, is_fits_file=True) not_yet_successful = False except: # I've only seen ValueError, but might as well catch for all errors and re-try time.sleep(pause_time_between_tries_sec) cur_try += 1 self.cur_fitsfile_basename = os.path.basename(filename) self.cur_fitsfile_path = os.path.abspath(os.path.dirname(filename)) self.set_window_title() if (hasattr(self.primary_image_panel, 'cur_fits_header_dialog') and self.primary_image_panel.cur_fits_header_dialog.is_dialog_still_open): raw_header_str = self.cur_fits_hdulist[0].header.tostring() header_str = (('\n'.join([raw_header_str[i:i+80] for i in np.arange(0, len(raw_header_str), 80) if raw_header_str[i:i+80] != " "*80])) + '\n') self.primary_image_panel.cur_fits_header_dialog.SetTitle(self.cur_fitsfile_basename) self.primary_image_panel.cur_fits_header_dialog.text.SetValue(header_str) self.primary_image_panel.cur_fits_header_dialog.last_find_index = 0 self.primary_image_panel.cur_fits_header_dialog.on_search(None) # TODO: better error handling for if WCS not available or partially available try: w = wcs.WCS(self.cur_fits_hdulist[0].header) # TODO: (urgent) need to check ones/arange in following, do I have this reversed? a = w.all_pix2world( np.outer(np.ones(self.raw_image.shape[-2]), np.arange(self.raw_image.shape[-1])), np.outer(np.arange(self.raw_image.shape[-2]), np.ones(self.raw_image.shape[-1])), 0) self.image_radec = ICRS(a[0]*units.degree, a[1]*units.degree) except: # just ignore radec if anything at all goes wrong. self.image_radec = None wx.CallAfter(pub.sendMessage, 'fitsfile-loaded', msg=filename) else: raise Error("Cannot find file: {}".format(filename)) else: raise Error("Requested filename ({}) does not end with .fits, .fits.gz, " + "or other capitalization of those".format(filename)) else: raise Error("load_fits_file requires string input, not type: {}".format(type(filename)))
def plot_dust_overlay(nh3_image_fits,h2_image_fits,region,plot_param,v_min,v_max,maskLim,obsMaskFits): text_size = 14 b18_text_size = 20 # Contour parameters (currently NH3 moment 0) cont_color='black' cont_lw = 0.6 cont_levs=2**np.arange( 0,20)*plot_param['w11_step'] # Masking of small (noisy) regions selem = np.array([[0,1,0],[1,1,1],[0,1,0]]) LowestContour = cont_levs[0]*0.5 w11_hdu = fits.open(nh3_image_fits) map = w11_hdu[0].data mask = binary_opening(map > LowestContour, selem) MaskedMap = mask*map w11_hdu[0].data = MaskedMap # Labels if region == 'B18': label_colour = 'white' text_size = b18_text_size else: label_colour = 'white' fig=aplpy.FITSFigure(h2_image_fits,figsize=(plot_param['size_x'], plot_param['size_y'])) fig.show_colorscale(cmap='hot',vmin=v_min,vmax=v_max,stretch='log',vmid=v_min-v_min*plot_param['vmid_scale']) fig.set_nan_color('0.95') # Observations mask contour fig.show_contour(obsMaskFits,colors='white',levels=1,linewidths=1.5) # NH3 moment contours fig.show_contour(w11_hdu,colors=cont_color,levels=cont_levs,linewidths=cont_lw) fig.axis_labels.set_font(family='sans_serif',size=text_size) # Ticks fig.tick_labels.set_font(family='sans_serif',size=text_size) fig.ticks.set_color('white') fig.tick_labels.set_xformat('hh:mm:ss') fig.tick_labels.set_style('colons') fig.tick_labels.set_yformat('dd:mm') # Scale bar # magic line of code to obtain scale in arcsec obtained from # http://www.astropy.org/astropy-tutorials/Quantities.html ang_sep = (plot_param['scalebar_size'].to(u.au)/plot_param['distance']).to(u.arcsec, equivalencies=u.dimensionless_angles()) fig.add_scalebar(ang_sep.to(u.degree)) fig.scalebar.set_font(family='sans_serif',size=text_size) fig.scalebar.set_corner(plot_param['scalebar_pos']) fig.scalebar.set(color='white') fig.scalebar.set_label('{0:4.2f}'.format(plot_param['scalebar_size'])) fig.add_label(plot_param['label_xpos'], plot_param['label_ypos'], '{0}\n{1}'.format(region,'500 $\mu$m'), relative=True, color=label_colour, horizontalalignment=plot_param['label_align'], family='sans_serif',size=text_size) fig.save( 'figures/{0}_continuum_image.pdf'.format(region),adjust_bbox=True,dpi=200)#, bbox_inches='tight') fig.close()
def plot_abundance(nh3_cont_fits,nh3_col_hdu,h2_col_hdu,region,plot_pars,maskLim,obsMaskFits): text_size = 14 b18_text_size = 20 if region == 'B18': text_size = b18_text_size # Get protostellar locations ra_prot, de_prot = get_prot_loc(region) # Contour parameters (currently NH3 moment 0) cont_color='0.6' cont_lw = 0.6 cont_levs=2**np.arange( 0,20)*plot_param['w11_step'] # Calculate abundance log_xnh3 = nh3_col_hdu[0].data - np.log10(h2_col_hdu.data) log_xnh3_hdu = fits.PrimaryHDU(log_xnh3,nh3_col_hdu[0].header) log_xnh3_hdu.writeto('../testing/{0}/parameterMaps/{0}_XNH3_{1}.fits'.format(region,file_extension),clobber=True) fig=aplpy.FITSFigure(log_xnh3_hdu,figsize=(plot_param['size_x'], plot_param['size_y'])) fig.show_colorscale(cmap='YlOrRd_r',vmin=plot_param['xnh3_lim'][0],vmax=plot_param['xnh3_lim'][1]) #fig.set_nan_color('0.95') # Observations mask contour fig.show_contour(obsMaskFits,colors='white',levels=1,linewidths=1.5) # NH3 moment contours # Masking of small (noisy) regions selem = np.array([[0,1,0],[1,1,1],[0,1,0]]) LowestContour = cont_levs[0]*0.5 w11_hdu = fits.open(nh3_cont_fits) map = w11_hdu[0].data mask = binary_opening(map > LowestContour, selem) MaskedMap = mask*map w11_hdu[0].data = MaskedMap fig.show_contour(w11_hdu,colors=cont_color,levels=cont_levs,linewidths=cont_lw) # Ticks fig.ticks.set_color('black') fig.tick_labels.set_font(family='sans_serif',size=text_size) fig.tick_labels.set_xformat('hh:mm:ss') fig.tick_labels.set_style('colons') fig.tick_labels.set_yformat('dd:mm') # Scale bar ang_sep = (plot_param['scalebar_size'].to(u.au)/plot_param['distance']).to(u.arcsec, equivalencies=u.dimensionless_angles()) fig.add_colorbar() fig.colorbar.show(box_orientation='horizontal', width=0.1, pad=0.0, location='top', ticks=[-10,-9.5,-9,-8.5,-8,-7.5,-7,-6.5]) fig.colorbar.set_font(family='sans_serif',size=text_size) fig.add_scalebar(ang_sep.to(u.degree)) fig.scalebar.set_font(family='sans_serif',size=text_size) fig.scalebar.set_corner(plot_param['scalebar_pos']) fig.scalebar.set(color='black') fig.scalebar.set_label('{0:4.2f}'.format(plot_param['scalebar_size'])) label_colour = 'black' fig.add_label(plot_param['label_xpos'], plot_param['label_ypos'], '{0}\n{1}'.format(region,r'$\mathrm{log} \ X(\mathrm{NH}_3)$'), relative=True, color=label_colour, horizontalalignment=plot_param['label_align'], family='sans_serif',size=text_size) fig.save( 'figures/{0}_xnh3_image.pdf'.format(region),adjust_bbox=True,dpi=200)#, bbox_inches='tight') # Add protostars fig.show_markers(ra_prot,de_prot,marker='*',s=50, c='white',edgecolors='black',linewidth=0.5,zorder=4) fig.save( 'figures/{0}_xnh3_image_prot.pdf'.format(region),adjust_bbox=True,dpi=200) fig.close()
