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
python类degree()的实例源码
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 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()
