In [37]:
import numpy as np
import matplotlib.pyplot as plt
x = np.arange(1,5,0.1)
y1 = np.sin(x)
y2 = np.cos(x)
fig = plt.figure(figsize=(9,6))
ax = fig.gca()
ax.plot(x,y1,'r-')
ax.plot(x,y2,'bo')
plt.show()
创建 Figure：
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(4,5))
<Figure size 288x360 with 0 Axes>
以上代码新建了一个 figure 实例。参数 figszie=(4,5) 代表图像的宽、高分别为 400、500 像素。
创建 Axes：在 figure 实例中增加一个 axes 实例：
ax = fig.gca()
gca 代表 get current axes，是最简单的办法。
也可以用：
ax = fig.add_subplot(111)
其中 111 代表本 axes 实例是一行、一列的第一个元素。或者
ax = fig.add_axes([0.15,0.15,0.8,0.8])
用这种方法可以指定坐标轴的大小和位置。 [0.15,0.15,0.8,0.8] 代表到边界的距离和坐标轴的宽、高占整个图片宽高的百分比。
axes 类的 plot() 方法是最常用的绘制平面图的方法。例如：
import numpy as np
import matplotlib.pyplot as plt
x = np.arange(1,5,0.1)
y1 = np.sin(x)
y2 = np.cos(x)
fig = plt.figure(figsize=(9,6))
ax = fig.gca()
ax.plot(x,y1,'r-')
ax.plot(x,y2,'bo')
plt.show()
axes.plot() 方法支持的参数有：
fmt='ro-'。这里支持线和点的标记结合，比如 -（直线）、--（虚线）、:（点线）markerfacecolor='b'markeredgecolor='b'markeredgewidth=1markersize=2.5linewidth=0.6alpha=0.5label='sin(x)'Matplotlib 支持的颜色代码有：
Matplotlib 支持的标记的形状有：
除此之外还有 x（叉）、d（小菱形）、D（大菱形）、h（竖六边形） H（横六边形）
Matplotlib 支持的线的形状、线的端点的形状、虚线的端点的形状有：
画直方图可以用 axes.hist() 方法。例如：
import numpy as np
x = np.random.rand(100)
fig = plt.figure(figsize=(9,6))
ax = fig.gca()
ax.hist(x, bins=10)
plt.show()
画阶梯图可以用 axes.step() 方法。例如：
import numpy as np
x = np.arange(1,5,0.1)
y1 = np.sin(x)
y2 = np.cos(x)
fig = plt.figure(figsize=(9,6))
ax = fig.gca()
ax.step(x, y1, where='mid')
plt.show()
其中 where 可取值有 ‘pre’, ‘post’, ‘mid’ 分别表示阶梯的“平台”在这个 x 值之前、之后还是中间。
Axes 的 errorbar() 方法绘制误差棒。例如：
import numpy as np
x = np.arange(1,5,0.1)
y = np.sin(x)
errors = np.random.rand(x.size)
fig = plt.figure(figsize=(9,6))
ax = fig.gca()
ax.errorbar(x, y, errors, fmt='ro-')
plt.show()
errorbar() 方法的完整调用：
matplotlib.pyplot.errorbar(x, y, yerr=None, xerr=None, fmt='', ecolor=None, elinewidth=None, capsize=None, barsabove=False, lolims=False, uplims=False, xlolims=False, xuplims=False, errorevery=1, capthick=None, *, data=None, **kwargs)
支持的参数：
rcParams["errorbar.capsize"] 取值可以用 axes 对象中的 plot_date() 方法画时间序列图。
字符串转换成时间序列有两种方法：
dateutil 包的 parser.parse() 函数将字符串处理成由 datetime 对象组成的列表。接着可以用 pylab 包的 date2num() 函数或者 matplotlib 包 dates 模块里的 date2num() 把时间序列转换为浮点数，用作 plot_date() 的输入。例如：import dateutil.parser
import matplotlib.dates as mdates
import numpy as np
x_lst = np.arange(7)
y_lst = np.sin(x_lst)
date_lst = ['2018-01-01T22:45:56', '2018-01-01T23:56:56', '2018-01-02T01:12:21', '2018-01-02T02:30:14',
'2018-01-02T03:33:21', '2018-01-02T04:35:45', '2018-01-02T05:40:11']
dates = [dateutil.parser.parse(s) for s in date_lst]
datenums = mdates.date2num(dates)
fig = plt.figure(figsize=(9,6))
ax = fig.gca()
ax.plot_date(datenums, y_lst ,'r-')
plt.show()
matplotlib 包 dates 模块里的 datestr2num() 函数把字符串序列转换成浮点数另一种方法是把 x 或 y 转换成浮点数后，先用普通的 plot() 方法画图，而后再用 xaxis_date() 或 yaxis_date() 方法把对应的轴转换为时间轴。例如：import matplotlib.dates as mdates
date = ['3 Jan 2013', '4 Jan 2013', '5 Jan 2013', '6 Jan 2013', '7 Jan 2013']
time = ['16:44:00', '16:45:00', '16:46:00', '16:47:00', '16:48:00']
x = mdates.datestr2num(date)
y = mdates.datestr2num(time)
fig = plt.figure(figsize=(9,6))
ax = fig.gca()
ax.plot(x, y, 'ro-')
ax.xaxis_date()
ax.yaxis_date()
plt.show()
格式化时间轴：在 plot_date() 之后可以用 fig.autofmt_xdate() 自动格式化时间轴。
指定时间轴刻度的大小需要使用 matplotlib 包的 dates 模块中的一些列时间间隔对象。例如：
import matplotlib.dates as mdates
x_lst = np.arange(7)
y_lst = np.sin(x_lst)
date_lst = ['2018-01-01T22:45:56', '2018-01-01T23:56:56', '2018-01-02T01:12:21', '2018-01-02T02:30:14',
'2018-01-02T03:33:21', '2018-01-02T04:35:45', '2018-01-02T05:40:11']
dates = [dateutil.parser.parse(s) for s in date_lst]
datenums = mdates.date2num(dates)
fig = plt.figure(figsize=(9,6))
ax = fig.gca()
ax.plot(datenums, y_lst, 'r-')
ax.xaxis_date()
ax.xaxis.set_major_locator(mdates.HourLocator())
ax.xaxis.set_minor_locator(mdates.MinuteLocator(byminute=np.arange(0,60,10)))
plt.show()
上面的例子表示大刻度为每个小时一个，小刻度为每10分钟一个。类似的 locator 还有：
例如，YearLocator() 表示刻度位于每年的 1 月 1 日处；
YearLocator(5, month=7, day=4) 表示刻度位于每年的 7 月 4 日。
其余的用法略有区别，例如 MinuteLocator(byminute=np.arange(0, 10, 60))
表示刻度为每 10 分钟一个，位于 0, 10, 20, ... 50 处。
可以用 axes.scatter() 方法画不同颜色的散点图。
ax.scatter(x_lst, logg_lst, c=c_lst, s=radius_lst, alpha=0.75, cmap='jet')
其中各输入参数：
lw=0edgecolor='r' 把边框设为红色用 axes.imshow() 方法画平面二维图像
axes.imshow()
aspect='auto'Axes() 对象的 quiver() 方法可以绘制箭头。它的调用有 4 种方法：
quiver(U, V, **kw)
quiver(U, V, C, **kw)
quiver(X, Y, U, V, **kw)
quiver(X, Y, U, V, C, **kw)其中 X,Y 是箭头尾部的数据值（可以是 1-D 或者 2-D 数组）；U,V 是箭头矢量；C 是颜色。可选参数有：
可以用 figure.colorbar() 方法在 axes.imshow() 或者 axes.scatter() 画的图上添加颜色棒 (colorbar)
fig = plt.figure(figsize=(12,8))
ax = fig.gca()
data = np.random.rand(50,50)*50000
cax = ax.imshow(data)
cbar = fig.colorbar(cax)
plt.show()
在这一节中，我们结合天文绘图中的一些实例来演示Matplotlib的用法
首先导入需要用到的python软件包
import pandas as pd
from astropy.table import Table
import matplotlib.pyplot as plt
import numpy as np
接下来我们用已知的系外行星表格为例，演示散点图和天空投影图的绘制方法。
首先到 http://exoplanet.eu/ 网站上下载系外行星表格，文件名为 exoplanet.eu_catalog.csv。这是一个csv格式的文件，包含已知的系外行星的各种参数。
用 pandas 的 read_csv() 函数读入这个文件，并用 astropy.table 模块转换成 Astropy 格式的表格：
df = pd.read_csv('exoplanet.eu_catalog.csv')
t = Table.from_pandas(df)
print(t)
print(t.colnames)
# name planet_status ... star_alternate_names
---------- ------------- ... ----------------------------------------------
11 Com b Confirmed ... --
11 Oph b Confirmed ... Oph 1622-2405, Oph 11A
11 UMi b Confirmed ... --
14 And b Confirmed ... --
14 Her b Confirmed ... --
14 Her c Confirmed ... --
16 Cyg B b Confirmed ... --
... ... ... ...
tau Gem b Confirmed ... --
ups And b Confirmed ... --
ups And c Confirmed ... --
ups And d Confirmed ... --
ups And e Confirmed ... --
ups Leo b Confirmed ... --
zet Del B Confirmed ... HD 196180, HIP 101589, 2MASS J20351852+1440272
Length = 5156 rows
['# name', 'planet_status', 'mass', 'mass_error_min', 'mass_error_max', 'mass_sini', 'mass_sini_error_min', 'mass_sini_error_max', 'radius', 'radius_error_min', 'radius_error_max', 'orbital_period', 'orbital_period_error_min', 'orbital_period_error_max', 'semi_major_axis', 'semi_major_axis_error_min', 'semi_major_axis_error_max', 'eccentricity', 'eccentricity_error_min', 'eccentricity_error_max', 'inclination', 'inclination_error_min', 'inclination_error_max', 'angular_distance', 'discovered', 'updated', 'omega', 'omega_error_min', 'omega_error_max', 'tperi', 'tperi_error_min', 'tperi_error_max', 'tconj', 'tconj_error_min', 'tconj_error_max', 'tzero_tr', 'tzero_tr_error_min', 'tzero_tr_error_max', 'tzero_tr_sec', 'tzero_tr_sec_error_min', 'tzero_tr_sec_error_max', 'lambda_angle', 'lambda_angle_error_min', 'lambda_angle_error_max', 'impact_parameter', 'impact_parameter_error_min', 'impact_parameter_error_max', 'tzero_vr', 'tzero_vr_error_min', 'tzero_vr_error_max', 'k', 'k_error_min', 'k_error_max', 'temp_calculated', 'temp_calculated_error_min', 'temp_calculated_error_max', 'temp_measured', 'hot_point_lon', 'geometric_albedo', 'geometric_albedo_error_min', 'geometric_albedo_error_max', 'log_g', 'publication', 'detection_type', 'mass_detection_type', 'radius_detection_type', 'alternate_names', 'molecules', 'star_name', 'ra', 'dec', 'mag_v', 'mag_i', 'mag_j', 'mag_h', 'mag_k', 'star_distance', 'star_distance_error_min', 'star_distance_error_max', 'star_metallicity', 'star_metallicity_error_min', 'star_metallicity_error_max', 'star_mass', 'star_mass_error_min', 'star_mass_error_max', 'star_radius', 'star_radius_error_min', 'star_radius_error_max', 'star_sp_type', 'star_age', 'star_age_error_min', 'star_age_error_max', 'star_teff', 'star_teff_error_min', 'star_teff_error_max', 'star_detected_disc', 'star_magnetic_field', 'star_alternate_names']
接下来，我们选取其中用视向速度方法（Radial Velocity）和凌星方法（Transit）发现的系外行星，并提取他们的RA、Dec
m = t['detection_type']=='Radial Velocity'
t1 = t[m]
ra = t1['ra']
dec = t1['dec']
Vmag = t1['mag_v']
m2 = t['detection_type']=='Primary Transit'
t2 = t[m2]
ra2 = t2['ra']
dec2 = t2['dec']
Vmag2 = t2['mag_v']
接下来把这些系外行星的坐标用不同颜色画在一张图上
fig = plt.figure(figsize=(14,7))
ax = fig.gca()
ax.scatter(ra, dec, alpha=0.5, s=5)
ax.scatter(ra2, dec2, alpha=0.5, s=5)
ax.set_xlabel('RA (deg)')
ax.set_ylabel('Dec (deg)')
plt.show()
下面改用 Hammer 投影：
fig = plt.figure(figsize=(12,6))
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection='hammer')
ax.plot(np.pi-np.deg2rad(ra), np.deg2rad(dec), 'o', ms=3, alpha=0.6)
ax.plot(np.pi-np.deg2rad(ra2), np.deg2rad(dec2), 'o', ms=3, alpha=0.6)
plt.show()
银道面是一些银经（l）从 0 到 360 deg 递增的点，但是银纬均为 0。可以用 Astropy 中的 SkyCoord 对象初始化一系列银道面的坐标，然后将其转换到赤道坐标上。
from astropy.coordinates import SkyCoord
gc = SkyCoord(np.arange(360), 0, frame='galactic', unit='deg')
eq = gc.transform_to('icrs')
ra_lst = eq.ra.deg
dec_lst = eq.dec.deg
接下来画在上面的图上
fig = plt.figure(figsize=(12,6))
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection='hammer')
ax.plot(np.pi-np.deg2rad(ra), np.deg2rad(dec), 'o', ms=3, alpha=0.6)
ax.plot(np.pi-np.deg2rad(ra2), np.deg2rad(dec2), 'o', ms=3, alpha=0.6)
ax.plot(np.pi-np.deg2rad(ra_lst), np.deg2rad(dec_lst))
plt.show()
可以看到上方有一条额外的横线，这是因为在从银道坐标向赤道坐标转换的过程中出现了从0到360 deg的跳变
for _ra, _dec in zip(ra_lst, dec_lst):
print(_ra, _dec)
266.4049882865447 -28.936177761791473 266.9955176679128 -28.08135168736493 267.5767060161091 -27.224037841089295 268.1490022696926 -26.364362477475595 268.7128381626685 -25.5024462620737 269.26862915603925 -24.63840462743747 269.81677533218704 -23.77234810755297 270.3576622528005 -22.904382652207385 270.8916617812854 -22.03460992269079 271.41913287076875 -21.163127570140684 271.9404223189431 -20.290029497760873 272.4558654910891 -19.41540610807084 272.9657870126846 -18.53934453627135 273.47050143305205 -17.66192887074451 273.97031386151957 -16.783240361644236 274.4655205775834 -15.903357618473493 274.956409616554 -15.022356797490446 275.443261332159 -14.140311779733736 275.92634893755695 -13.257294340410404 276.4059390261891 -12.373374310345426 276.88229207387474 -11.48861973015191 277.3556629235181 -10.603096997743211 277.82630125376977 -9.716871009774811 278.29445203295046 -8.830005297571905 278.7603559595139 -7.942562158071049 279.22424989029605 -7.054602780278212 279.6863672577668 -6.166187367722703 280.1469384774753 -5.277375257365946 280.6061913468523 -4.388225035405868 281.06435143651026 -3.4987946504018 281.5216424751636 -2.6091415241310076 281.9782867292705 -1.7193226605764367 282.43450537848526 -0.829394753435806 282.89051888799634 0.06058570746566121 283.3465473788177 0.9505623304818671 283.80281099709407 1.840478716798579 284.259530283478 2.7302783550610554 284.71692654363744 3.6199045156046323 285.17522222095494 4.509300143915023 285.6346412724869 5.398407752940556 286.09540954925956 6.287169313873157 286.5577551819927 7.175526145006639 287.02190897335515 8.063418798271503 287.48810479787585 8.950786943033199 287.95658001065556 9.837569246726982 288.4275758660451 10.723703251886231 288.9013379474855 11.609125249102664 289.3781166097293 12.493770145435926 289.85816743469553 13.37757132776689 290.3417517022372 14.260460520562733 290.8291368771373 15.142367637493535 291.320597113675 16.02322062630844 291.8164137791409 16.902945306344822 292.3168759977044 17.781465198006536 292.82228121606835 18.658701343506113 293.3329357923663 19.534572118121435 293.84915560977703 20.408993031169693 294.3712667163406 21.281876515849557 294.89960599246234 22.15313170704686 295.4345218475784 23.022664206140213 295.976374947429 23.890375831778837 296.5255389733405 24.75616435553788 297.0824014148435 25.619923221284736 297.6473643968569 26.48154124701468 298.22084554252706 27.34090230783456 298.8032788726349 28.197884998691183 299.3951157422478 29.052362275354906 299.9968258149968 29.904201072081324 300.60889807499103 30.753261894284506 301.23184187591585 31.599398384465232 301.8661880262955 32.4424568595495 302.512489909207 33.282275817706584 303.1713246338876 34.118685412636765 303.84329421566446 34.95150689324701 304.5290267794107 35.78055200657473 305.2291777802729 36.60562236177787 305.9444312336805 37.42650875299223 306.6755009445863 38.24299043886812 307.42313172345973 39.0548343766507 308.1881005737029 39.86179440876841 308.97121783181626 40.663610400057614 309.7733282377545 41.46000732399134 310.59531190840033 42.25069429661651 311.43808518187967 43.03536355735449 312.3026012944694 43.81368939641419 313.1898508450357 44.58532702932762 314.10086199422443 45.349911420083856 315.036700336949 46.10705605554057 315.99846837705775 46.856351675279996 316.98730452241205 47.59736496288937 318.0043815070175 48.32963720684089 319.0509041344361 49.05268294177055 320.1281062236728 49.7659885840692 321.23724662538973 50.46901107935159 322.379604163124 51.16117658361253 323.55647134182135 51.84187920474915 324.7691466553278 52.51047983664887 326.01892531664 53.16630512420765 327.30708823116873 53.80864660441273 328.6348890357846 54.436760075904246 330.00353903714 55.049865257064695 331.4141899041588 55.64714580042597 332.86791400438676 56.22774973868711 334.3656823250085 56.79079044443061 335.9083399895893 57.335348191106654 337.49657947351216 57.860472406274155 339.13091173638287 58.365184708559596 340.8116356288618 58.84848281631341 342.5388060930252 59.309345407437824 344.3122018555012 59.74673799526114 346.13129350407814 60.15961986369133 347.9952130303879 60.54695207550242 349.9027260988443 60.90770653023386 351.8522084468951 61.24087600317357 353.84162791279533 61.54548504540372 355.86853360255543 61.820601569005376 357.9300536270824 62.06534888431203 0.022902648209691012 62.27891790156743 2.1434001609123046 62.46057916216022 4.287500012557363 62.60969432977033 6.4508311362392785 62.72572675402351 8.628748885938778 62.80825072250175 10.816395750803329 62.85695904358509 13.008769647053144 62.871668652954114 15.200797493950706 62.85232400866973 17.387411424977106 62.79899812924392 19.56362480457165 62.71189122970544 21.724605234214096 62.59132701473461 23.865741937774327 62.43774678744911 25.982705292539205 62.251701619880926 28.07149677932926 62.03384290051911 30.12848821148306 61.78491162147891 32.15044971302585 61.50572679123264 34.134566499087825 61.197173359051526 36.07844502375452 60.86019001698652 37.980109472029135 60.49575720840557 39.83798986773294 60.10488562363062 41.65090324563002 59.68860540805453 43.41802940255796 59.24795625087803 45.13888271532724 58.78397846716109 46.81328141339644 58.297705135174716 48.44131554380656 57.79015530701964 50.023314685320756 57.26232827418853 51.559816276095404 56.71519884144903 53.051535228736284 56.14971354179037 54.49933532863094 55.56678771146578 55.90420275268767 54.96730333640625 57.267221909647965 54.35210757839957 58.589553691078194 53.72201189035915 59.87241613331321 53.07779163375647 61.11706742313291 52.420186116992014 62.32479113123326 51.7498989804038 63.49688352474446 51.06759886117136 64.63464279025995 50.3739202791202 65.73935998937385 49.669464692027 66.81231156714722 48.95480167624669 67.85475323815822 48.23047019517276 68.86791508309094 47.496979924121426 69.85299769978133 46.75481260566137 70.81116926515999 46.00442341419401 71.7435633777629 45.246242312749445 72.6512775638052 44.48067538853924 73.53537234277759 43.70810615685175 74.39687076083838 42.92889682543728 75.23675831173807 42.14338951367192 76.05598317553378 41.3519074225532 76.85545671487732 40.55475595302824 77.63605417720434 39.752223771323294 78.39861555874239 38.944583820878634 79.1439465929495 38.13209428123125 79.8728198318594 37.314999474757904 80.58597579391446 36.493530722627646 81.28412415629028 35.66790715163271 81.96794497352619 34.83833645379704 82.6380899075459 34.00501560081409 83.29518345694454 33.168131515460225 83.93982417579302 32.32786170217454 84.57258587422368 31.484374839003653 85.1940187947557 30.63783133308915 85.80465075974658 29.78838384182976 86.40498828654475 28.93617776179148 86.99551766791282 28.081351687364922 87.57670601610907 27.2240378410893 88.14900226969256 26.364362477475606 88.71283816266855 25.5024462620737 89.26862915603932 24.638404627437485 89.81677533218705 23.772348107552965 90.35766225280055 22.904382652207392 90.89166178128544 22.034609922690787 91.41913287076878 21.16312757014069 91.94042231894313 20.290029497760862 92.4558654910891 19.415406108070844 92.96578701268464 18.53934453627135 93.47050143305206 17.66192887074451 93.9703138615196 16.783240361644243 94.46552057758343 15.903357618473493 94.956409616554 15.022356797490456 95.44326133215907 14.140311779733736 95.92634893755692 13.25729434041041 96.40593902618912 12.373374310345424 96.88229207387477 11.48861973015191 97.35566292351811 10.603096997743204 97.82630125376978 9.716871009774811 98.29445203295047 8.830005297571915 98.76035595951394 7.942562158071049 99.22424989029604 7.054602780278221 99.68636725776679 6.1661873677226975 100.14693847747536 5.277375257365946 100.60619134685236 4.388225035405862 101.06435143651031 3.498794650401806 101.52164247516362 2.6091415241309983 101.97828672927052 1.719322660576437 102.43450537848526 0.8293947534358155 102.89051888799634 -0.06058570746566121 103.34654737881772 -0.9505623304818517 103.8028109970941 -1.8404787167985823 104.25953028347801 -2.7302783550610497 104.71692654363747 -3.6199045156046354 105.17522222095499 -4.509300143915012 105.63464127248692 -5.398407752940567 106.09540954925957 -6.287169313873158 106.55775518199272 -7.175526145006628 107.02190897335512 -8.063418798271503 107.48810479787588 -8.95078694303319 107.95658001065557 -9.837569246726991 108.42757586604515 -10.723703251886228 108.90133794748553 -11.609125249102666 109.3781166097294 -12.49377014543592 109.85816743469559 -13.377571327766898 110.34175170223726 -14.260460520562734 110.82913687713727 -15.142367637493528 111.32059711367502 -16.023220626308422 111.81641377914093 -16.902945306344833 112.3168759977044 -17.781465198006543 112.82228121606833 -18.658701343506106 113.33293579236633 -19.53457211812141 113.84915560977704 -20.408993031169707 114.3712667163406 -21.281876515849557 114.89960599246238 -22.153131707046857 115.4345218475784 -23.0226642061402 115.976374947429 -23.890375831778826 116.52553897334049 -24.7561643555379 117.08240141484352 -25.619923221284747 117.64736439685689 -26.481541247014672 118.22084554252704 -27.34090230783454 118.80327887263498 -28.197884998691197 119.39511574224782 -29.052362275354902 119.99682581499681 -29.90420107208131 120.60889807499105 -30.7532618942845 121.23184187591588 -31.599398384465214 121.8661880262956 -32.44245685954953 122.512489909207 -33.28227581770659 123.17132463388754 -34.11868541263676 123.84329421566449 -34.951506893247 124.52902677941076 -35.780552006574744 125.22917778027293 -36.60562236177788 125.94443123368053 -37.42650875299223 126.67550094458629 -38.24299043886811 127.42313172345973 -39.05483437665069 128.18810057370297 -39.86179440876844 128.97121783181632 -40.663610400057614 129.7733282377545 -41.460007323991334 130.59531190840028 -42.25069429661648 131.4380851818797 -43.03536355735451 132.30260129446938 -43.81368939641419 133.1898508450357 -44.585327029327615 134.10086199422446 -45.349911420083835 135.03670033694902 -46.10705605554055 135.99846837705783 -46.85635167528 136.9873045224121 -47.59736496288937 138.0043815070175 -48.329637206840886 139.05090413443608 -49.05268294177054 140.12810622367283 -49.7659885840692 141.2372466253898 -50.46901107935159 142.379604163124 -51.16117658361252 143.55647134182138 -51.84187920474914 144.76914665532775 -52.510479836648834 146.01892531664 -53.16630512420767 147.3070882311688 -53.80864660441272 148.6348890357846 -54.43676007590423 150.00353903713997 -55.04986525706469 151.4141899041588 -55.647145800425996 152.86791400438682 -56.22774973868711 154.36568232500852 -56.79079044443061 155.9083399895893 -57.335348191106654 157.4965794735121 -57.86047240627414 159.13091173638296 -58.365184708559596 160.8116356288618 -58.84848281631341 162.53880609302524 -59.30934540743781 164.3122018555012 -59.746737995261114 166.13129350407814 -60.15961986369133 167.9952130303879 -60.54695207550242 169.90272609884434 -60.907706530233845 171.85220844689502 -61.24087600317357 173.84162791279527 -61.545485045403716 175.86853360255552 -61.820601569005376 177.9300536270824 -62.06534888431203 180.02290264820965 -62.27891790156743 182.1434001609122 -62.46057916216021 184.28750001255736 -62.60969432977034 186.4508311362393 -62.72572675402351 188.62874888593873 -62.80825072250175 190.8163957508033 -62.85695904358511 193.00876964705304 -62.871668652954114 195.20079749395072 -62.85232400866973 197.38741142497707 -62.79899812924392 199.5636248045716 -62.71189122970542 201.724605234214 -62.59132701473461 203.86574193777435 -62.43774678744911 205.9827052925392 -62.25170161988091 208.0714967793292 -62.03384290051911 210.12848821148302 -61.784911621478926 212.15044971302572 -61.50572679123264 214.13456649908787 -61.197173359051526 216.0784450237545 -60.86019001698651 217.98010947202908 -60.49575720840557 219.8379898677329 -60.104885623630636 221.65090324563005 -59.688605408054535 223.41802940255795 -59.24795625087803 225.1388827153272 -58.78397846716109 226.81328141339642 -58.2977051351747 228.44131554380647 -57.79015530701965 230.02331468532074 -57.262328274188526 231.55981627609538 -56.71519884144904 233.0515352287362 -56.149713541790376 234.4993353286309 -55.5667877114658 235.90420275268764 -54.96730333640624 237.26722190964804 -54.35210757839956 238.58955369107815 -53.72201189035915 239.87241613331315 -53.077791633756476 241.11706742313285 -52.420186116992056 242.32479113123327 -51.749898980403785 243.49688352474445 -51.06759886117136 244.63464279025993 -50.3739202791202 245.73935998937378 -49.66946469202703 246.81231156714722 -48.954801676246674 247.85475323815825 -48.23047019517274 248.86791508309085 -47.49697992412143 249.85299769978127 -46.75481260566138 250.81116926515998 -46.00442341419403 251.7435633777629 -45.24624231274942 252.65127756380517 -44.48067538853925 253.53537234277758 -43.708106156851755 254.39687076083828 -42.9288968254373 255.23675831173807 -42.14338951367189 256.05598317553375 -41.351907422553175 256.8554567148773 -40.55475595302825 257.63605417720436 -39.7522237713233 258.39861555874234 -38.94458382087866 259.14394659294953 -38.13209428123123 259.8728198318594 -37.31499947475791 260.5859757939144 -36.49353072262765 261.28412415629026 -35.66790715163274 261.9679449735262 -34.83833645379703 262.63808990754586 -34.00501560081407 263.2951834569445 -33.16813151546023 263.93982417579304 -32.32786170217454 264.57258587422365 -31.484374839003685 265.19401879475566 -30.637831333089142 265.80465075974655 -29.788383841829763
可以用 numpy 中的 roll() 函数将这个数组滚动到跳变所在的位置
imaxdiff = np.abs(np.diff(ra_lst)).argmax()
print(imaxdiff)
ra_lst = np.roll(ra_lst, -imaxdiff-1)
dec_lst = np.roll(dec_lst, -imaxdiff-1)
for _ra, _dec in zip(ra_lst, dec_lst):
print(_ra, _dec)
116 0.022902648209691012 62.27891790156743 2.1434001609123046 62.46057916216022 4.287500012557363 62.60969432977033 6.4508311362392785 62.72572675402351 8.628748885938778 62.80825072250175 10.816395750803329 62.85695904358509 13.008769647053144 62.871668652954114 15.200797493950706 62.85232400866973 17.387411424977106 62.79899812924392 19.56362480457165 62.71189122970544 21.724605234214096 62.59132701473461 23.865741937774327 62.43774678744911 25.982705292539205 62.251701619880926 28.07149677932926 62.03384290051911 30.12848821148306 61.78491162147891 32.15044971302585 61.50572679123264 34.134566499087825 61.197173359051526 36.07844502375452 60.86019001698652 37.980109472029135 60.49575720840557 39.83798986773294 60.10488562363062 41.65090324563002 59.68860540805453 43.41802940255796 59.24795625087803 45.13888271532724 58.78397846716109 46.81328141339644 58.297705135174716 48.44131554380656 57.79015530701964 50.023314685320756 57.26232827418853 51.559816276095404 56.71519884144903 53.051535228736284 56.14971354179037 54.49933532863094 55.56678771146578 55.90420275268767 54.96730333640625 57.267221909647965 54.35210757839957 58.589553691078194 53.72201189035915 59.87241613331321 53.07779163375647 61.11706742313291 52.420186116992014 62.32479113123326 51.7498989804038 63.49688352474446 51.06759886117136 64.63464279025995 50.3739202791202 65.73935998937385 49.669464692027 66.81231156714722 48.95480167624669 67.85475323815822 48.23047019517276 68.86791508309094 47.496979924121426 69.85299769978133 46.75481260566137 70.81116926515999 46.00442341419401 71.7435633777629 45.246242312749445 72.6512775638052 44.48067538853924 73.53537234277759 43.70810615685175 74.39687076083838 42.92889682543728 75.23675831173807 42.14338951367192 76.05598317553378 41.3519074225532 76.85545671487732 40.55475595302824 77.63605417720434 39.752223771323294 78.39861555874239 38.944583820878634 79.1439465929495 38.13209428123125 79.8728198318594 37.314999474757904 80.58597579391446 36.493530722627646 81.28412415629028 35.66790715163271 81.96794497352619 34.83833645379704 82.6380899075459 34.00501560081409 83.29518345694454 33.168131515460225 83.93982417579302 32.32786170217454 84.57258587422368 31.484374839003653 85.1940187947557 30.63783133308915 85.80465075974658 29.78838384182976 86.40498828654475 28.93617776179148 86.99551766791282 28.081351687364922 87.57670601610907 27.2240378410893 88.14900226969256 26.364362477475606 88.71283816266855 25.5024462620737 89.26862915603932 24.638404627437485 89.81677533218705 23.772348107552965 90.35766225280055 22.904382652207392 90.89166178128544 22.034609922690787 91.41913287076878 21.16312757014069 91.94042231894313 20.290029497760862 92.4558654910891 19.415406108070844 92.96578701268464 18.53934453627135 93.47050143305206 17.66192887074451 93.9703138615196 16.783240361644243 94.46552057758343 15.903357618473493 94.956409616554 15.022356797490456 95.44326133215907 14.140311779733736 95.92634893755692 13.25729434041041 96.40593902618912 12.373374310345424 96.88229207387477 11.48861973015191 97.35566292351811 10.603096997743204 97.82630125376978 9.716871009774811 98.29445203295047 8.830005297571915 98.76035595951394 7.942562158071049 99.22424989029604 7.054602780278221 99.68636725776679 6.1661873677226975 100.14693847747536 5.277375257365946 100.60619134685236 4.388225035405862 101.06435143651031 3.498794650401806 101.52164247516362 2.6091415241309983 101.97828672927052 1.719322660576437 102.43450537848526 0.8293947534358155 102.89051888799634 -0.06058570746566121 103.34654737881772 -0.9505623304818517 103.8028109970941 -1.8404787167985823 104.25953028347801 -2.7302783550610497 104.71692654363747 -3.6199045156046354 105.17522222095499 -4.509300143915012 105.63464127248692 -5.398407752940567 106.09540954925957 -6.287169313873158 106.55775518199272 -7.175526145006628 107.02190897335512 -8.063418798271503 107.48810479787588 -8.95078694303319 107.95658001065557 -9.837569246726991 108.42757586604515 -10.723703251886228 108.90133794748553 -11.609125249102666 109.3781166097294 -12.49377014543592 109.85816743469559 -13.377571327766898 110.34175170223726 -14.260460520562734 110.82913687713727 -15.142367637493528 111.32059711367502 -16.023220626308422 111.81641377914093 -16.902945306344833 112.3168759977044 -17.781465198006543 112.82228121606833 -18.658701343506106 113.33293579236633 -19.53457211812141 113.84915560977704 -20.408993031169707 114.3712667163406 -21.281876515849557 114.89960599246238 -22.153131707046857 115.4345218475784 -23.0226642061402 115.976374947429 -23.890375831778826 116.52553897334049 -24.7561643555379 117.08240141484352 -25.619923221284747 117.64736439685689 -26.481541247014672 118.22084554252704 -27.34090230783454 118.80327887263498 -28.197884998691197 119.39511574224782 -29.052362275354902 119.99682581499681 -29.90420107208131 120.60889807499105 -30.7532618942845 121.23184187591588 -31.599398384465214 121.8661880262956 -32.44245685954953 122.512489909207 -33.28227581770659 123.17132463388754 -34.11868541263676 123.84329421566449 -34.951506893247 124.52902677941076 -35.780552006574744 125.22917778027293 -36.60562236177788 125.94443123368053 -37.42650875299223 126.67550094458629 -38.24299043886811 127.42313172345973 -39.05483437665069 128.18810057370297 -39.86179440876844 128.97121783181632 -40.663610400057614 129.7733282377545 -41.460007323991334 130.59531190840028 -42.25069429661648 131.4380851818797 -43.03536355735451 132.30260129446938 -43.81368939641419 133.1898508450357 -44.585327029327615 134.10086199422446 -45.349911420083835 135.03670033694902 -46.10705605554055 135.99846837705783 -46.85635167528 136.9873045224121 -47.59736496288937 138.0043815070175 -48.329637206840886 139.05090413443608 -49.05268294177054 140.12810622367283 -49.7659885840692 141.2372466253898 -50.46901107935159 142.379604163124 -51.16117658361252 143.55647134182138 -51.84187920474914 144.76914665532775 -52.510479836648834 146.01892531664 -53.16630512420767 147.3070882311688 -53.80864660441272 148.6348890357846 -54.43676007590423 150.00353903713997 -55.04986525706469 151.4141899041588 -55.647145800425996 152.86791400438682 -56.22774973868711 154.36568232500852 -56.79079044443061 155.9083399895893 -57.335348191106654 157.4965794735121 -57.86047240627414 159.13091173638296 -58.365184708559596 160.8116356288618 -58.84848281631341 162.53880609302524 -59.30934540743781 164.3122018555012 -59.746737995261114 166.13129350407814 -60.15961986369133 167.9952130303879 -60.54695207550242 169.90272609884434 -60.907706530233845 171.85220844689502 -61.24087600317357 173.84162791279527 -61.545485045403716 175.86853360255552 -61.820601569005376 177.9300536270824 -62.06534888431203 180.02290264820965 -62.27891790156743 182.1434001609122 -62.46057916216021 184.28750001255736 -62.60969432977034 186.4508311362393 -62.72572675402351 188.62874888593873 -62.80825072250175 190.8163957508033 -62.85695904358511 193.00876964705304 -62.871668652954114 195.20079749395072 -62.85232400866973 197.38741142497707 -62.79899812924392 199.5636248045716 -62.71189122970542 201.724605234214 -62.59132701473461 203.86574193777435 -62.43774678744911 205.9827052925392 -62.25170161988091 208.0714967793292 -62.03384290051911 210.12848821148302 -61.784911621478926 212.15044971302572 -61.50572679123264 214.13456649908787 -61.197173359051526 216.0784450237545 -60.86019001698651 217.98010947202908 -60.49575720840557 219.8379898677329 -60.104885623630636 221.65090324563005 -59.688605408054535 223.41802940255795 -59.24795625087803 225.1388827153272 -58.78397846716109 226.81328141339642 -58.2977051351747 228.44131554380647 -57.79015530701965 230.02331468532074 -57.262328274188526 231.55981627609538 -56.71519884144904 233.0515352287362 -56.149713541790376 234.4993353286309 -55.5667877114658 235.90420275268764 -54.96730333640624 237.26722190964804 -54.35210757839956 238.58955369107815 -53.72201189035915 239.87241613331315 -53.077791633756476 241.11706742313285 -52.420186116992056 242.32479113123327 -51.749898980403785 243.49688352474445 -51.06759886117136 244.63464279025993 -50.3739202791202 245.73935998937378 -49.66946469202703 246.81231156714722 -48.954801676246674 247.85475323815825 -48.23047019517274 248.86791508309085 -47.49697992412143 249.85299769978127 -46.75481260566138 250.81116926515998 -46.00442341419403 251.7435633777629 -45.24624231274942 252.65127756380517 -44.48067538853925 253.53537234277758 -43.708106156851755 254.39687076083828 -42.9288968254373 255.23675831173807 -42.14338951367189 256.05598317553375 -41.351907422553175 256.8554567148773 -40.55475595302825 257.63605417720436 -39.7522237713233 258.39861555874234 -38.94458382087866 259.14394659294953 -38.13209428123123 259.8728198318594 -37.31499947475791 260.5859757939144 -36.49353072262765 261.28412415629026 -35.66790715163274 261.9679449735262 -34.83833645379703 262.63808990754586 -34.00501560081407 263.2951834569445 -33.16813151546023 263.93982417579304 -32.32786170217454 264.57258587422365 -31.484374839003685 265.19401879475566 -30.637831333089142 265.80465075974655 -29.788383841829763 266.4049882865447 -28.936177761791473 266.9955176679128 -28.08135168736493 267.5767060161091 -27.224037841089295 268.1490022696926 -26.364362477475595 268.7128381626685 -25.5024462620737 269.26862915603925 -24.63840462743747 269.81677533218704 -23.77234810755297 270.3576622528005 -22.904382652207385 270.8916617812854 -22.03460992269079 271.41913287076875 -21.163127570140684 271.9404223189431 -20.290029497760873 272.4558654910891 -19.41540610807084 272.9657870126846 -18.53934453627135 273.47050143305205 -17.66192887074451 273.97031386151957 -16.783240361644236 274.4655205775834 -15.903357618473493 274.956409616554 -15.022356797490446 275.443261332159 -14.140311779733736 275.92634893755695 -13.257294340410404 276.4059390261891 -12.373374310345426 276.88229207387474 -11.48861973015191 277.3556629235181 -10.603096997743211 277.82630125376977 -9.716871009774811 278.29445203295046 -8.830005297571905 278.7603559595139 -7.942562158071049 279.22424989029605 -7.054602780278212 279.6863672577668 -6.166187367722703 280.1469384774753 -5.277375257365946 280.6061913468523 -4.388225035405868 281.06435143651026 -3.4987946504018 281.5216424751636 -2.6091415241310076 281.9782867292705 -1.7193226605764367 282.43450537848526 -0.829394753435806 282.89051888799634 0.06058570746566121 283.3465473788177 0.9505623304818671 283.80281099709407 1.840478716798579 284.259530283478 2.7302783550610554 284.71692654363744 3.6199045156046323 285.17522222095494 4.509300143915023 285.6346412724869 5.398407752940556 286.09540954925956 6.287169313873157 286.5577551819927 7.175526145006639 287.02190897335515 8.063418798271503 287.48810479787585 8.950786943033199 287.95658001065556 9.837569246726982 288.4275758660451 10.723703251886231 288.9013379474855 11.609125249102664 289.3781166097293 12.493770145435926 289.85816743469553 13.37757132776689 290.3417517022372 14.260460520562733 290.8291368771373 15.142367637493535 291.320597113675 16.02322062630844 291.8164137791409 16.902945306344822 292.3168759977044 17.781465198006536 292.82228121606835 18.658701343506113 293.3329357923663 19.534572118121435 293.84915560977703 20.408993031169693 294.3712667163406 21.281876515849557 294.89960599246234 22.15313170704686 295.4345218475784 23.022664206140213 295.976374947429 23.890375831778837 296.5255389733405 24.75616435553788 297.0824014148435 25.619923221284736 297.6473643968569 26.48154124701468 298.22084554252706 27.34090230783456 298.8032788726349 28.197884998691183 299.3951157422478 29.052362275354906 299.9968258149968 29.904201072081324 300.60889807499103 30.753261894284506 301.23184187591585 31.599398384465232 301.8661880262955 32.4424568595495 302.512489909207 33.282275817706584 303.1713246338876 34.118685412636765 303.84329421566446 34.95150689324701 304.5290267794107 35.78055200657473 305.2291777802729 36.60562236177787 305.9444312336805 37.42650875299223 306.6755009445863 38.24299043886812 307.42313172345973 39.0548343766507 308.1881005737029 39.86179440876841 308.97121783181626 40.663610400057614 309.7733282377545 41.46000732399134 310.59531190840033 42.25069429661651 311.43808518187967 43.03536355735449 312.3026012944694 43.81368939641419 313.1898508450357 44.58532702932762 314.10086199422443 45.349911420083856 315.036700336949 46.10705605554057 315.99846837705775 46.856351675279996 316.98730452241205 47.59736496288937 318.0043815070175 48.32963720684089 319.0509041344361 49.05268294177055 320.1281062236728 49.7659885840692 321.23724662538973 50.46901107935159 322.379604163124 51.16117658361253 323.55647134182135 51.84187920474915 324.7691466553278 52.51047983664887 326.01892531664 53.16630512420765 327.30708823116873 53.80864660441273 328.6348890357846 54.436760075904246 330.00353903714 55.049865257064695 331.4141899041588 55.64714580042597 332.86791400438676 56.22774973868711 334.3656823250085 56.79079044443061 335.9083399895893 57.335348191106654 337.49657947351216 57.860472406274155 339.13091173638287 58.365184708559596 340.8116356288618 58.84848281631341 342.5388060930252 59.309345407437824 344.3122018555012 59.74673799526114 346.13129350407814 60.15961986369133 347.9952130303879 60.54695207550242 349.9027260988443 60.90770653023386 351.8522084468951 61.24087600317357 353.84162791279533 61.54548504540372 355.86853360255543 61.820601569005376 357.9300536270824 62.06534888431203
接下来再次画这个图，就可以去掉多余的横线。
fig = plt.figure(figsize=(12,6))
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection='hammer')
ax.plot(np.pi-np.deg2rad(ra_lst), np.deg2rad(dec_lst))
ax.plot(np.pi-np.deg2rad(ra), np.deg2rad(dec), 'o', ms=3, alpha=0.6)
ax.plot(np.pi-np.deg2rad(ra2), np.deg2rad(dec2), 'o', ms=3, alpha=0.6)
ax.grid(True)
plt.show()
同理，我们也可以将黄道面画在这张图上。
ec = SkyCoord(np.arange(360), 0, frame='barycentrictrueecliptic', unit='deg')
eq = ec.transform_to('icrs')
ra_lst = eq.ra.deg
dec_lst = eq.dec.deg
fig = plt.figure(figsize=(12,6))
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection='hammer'
ax.plot(np.pi-np.deg2rad(ra), np.deg2rad(dec), 'o', ms=3, alpha=0.6)
ax.plot(np.pi-np.deg2rad(ra2), np.deg2rad(dec2), 'o', ms=3, alpha=0.6)
ax.plot(np.pi-np.deg2rad(ra_lst), np.deg2rad(dec_lst))
ax.grid(True)
plt.show()
在这一节中，我们以TESS光变曲线为例，演示 WCS 投影图的做法
TESS 卫星全称是 Transiting Exoplanet Survey Satellite，目的是用凌星方法在全天搜寻太阳系外行星，2018年发射升空。截止到 2023年6月份，TESS 已经扫描了天空 64 个 24° × 96° 的区域（称为 sector），发布了 40 多万颗目标星的 2min 间隔的光变曲线。数据发布在 https://archive.stsci.edu/tess/bulk_downloads/bulk_downloads_ffi-tp-lc-dv.html
下面我们用一颗已知的带有系外行星的恒星 XO-2N （TIC 356473034）为例，演示光变曲线和 TESS 发布的 target pixel 文件的绘制方法。
首先导入需要的软件包
import numpy as np
import astropy.io.fits as fits
import matplotlib.pyplot as plt
读入光变曲线文件。
fname = 'tess2019357164649-s0020-0000000356473034-0165-s_lc.fits'
data = fits.getdata(fname)
print(data.dtype)
(numpy.record, [('TIME', '>f8'), ('TIMECORR', '>f4'), ('CADENCENO', '>i4'), ('SAP_FLUX', '>f4'), ('SAP_FLUX_ERR', '>f4'), ('SAP_BKG', '>f4'), ('SAP_BKG_ERR', '>f4'), ('PDCSAP_FLUX', '>f4'), ('PDCSAP_FLUX_ERR', '>f4'), ('QUALITY', '>i4'), ('PSF_CENTR1', '>f8'), ('PSF_CENTR1_ERR', '>f4'), ('PSF_CENTR2', '>f8'), ('PSF_CENTR2_ERR', '>f4'), ('MOM_CENTR1', '>f8'), ('MOM_CENTR1_ERR', '>f4'), ('MOM_CENTR2', '>f8'), ('MOM_CENTR2_ERR', '>f4'), ('POS_CORR1', '>f4'), ('POS_CORR2', '>f4')])
光变曲线数据存储在这个FITS文件的二进制表格中，'TIME' 列代表时间，'PDCSAP_FLUX' 是流量，'QUALITY' 代表数据点的质量。QUALITY=0 代表数据点正常。
t_lst = data['TIME']
f_lst = data['PDCSAP_FLUX']
q_lst = data['QUALITY']
绘制光变曲线：
m = q_lst==0
t_lst = t_lst[m]
f_lst = f_lst[m]
fig = plt.figure(figsize=(14, 6))
ax = fig.gca()
ax.plot(t_lst, f_lst, lw=0.5)
plt.show()
TP (Target Pixel) 文件包含了目标星周围若干个像素的时间序列文件，有助于判断光变是否来自于目标星，还是周围的背景星。
filename = 'tess2019357164649-s0020-0000000356473034-0165-s_tp.fits'
hdulst = fits.open(filename)
print(hdulst[1].data.dtype)
(numpy.record, [('TIME', '>f8'), ('TIMECORR', '>f4'), ('CADENCENO', '>i4'), ('RAW_CNTS', '>i4', (11, 11)), ('FLUX', '>f4', (11, 11)), ('FLUX_ERR', '>f4', (11, 11)), ('FLUX_BKG', '>f4', (11, 11)), ('FLUX_BKG_ERR', '>f4', (11, 11)), ('QUALITY', '>i4'), ('POS_CORR1', '>f4'), ('POS_CORR2', '>f4')])
这个FITS文件的第二个 HDU 的'FLUX'列包含了目标星周围的图像数据。
t_lst = hdulst[1].data['TIME']
img_lst = hdulst[1].data['FLUX']
选取 img_lst 中的第一幅进行绘图
fig = plt.figure(figsize=(6,6))
ax = fig.gca()
ax.imshow(img_lst[0])
<matplotlib.image.AxesImage at 0x7f2768b44250>
import astropy.io.fits as fits
from astropy.wcs import WCS
head1 = hdulst[1].header
data = hdulst[1].data
print(head1)
img = data['FLUX'][0]
ny, nx = img.shape
print(ny, nx)
XTENSION= 'BINTABLE' / marks the beginning of a new HDU BITPIX = 8 / array data type NAXIS = 2 / number of array dimensions NAXIS1 = 2448 / length of first array dimension NAXIS2 = 18954 / length of second array dimension PCOUNT = 0 / group parameter count (not used) GCOUNT = 1 / group count (not used) TFIELDS = 11 / number of table fields TTYPE1 = 'TIME ' / column title: data time stamps TFORM1 = 'D ' / column format: 64-bit floating point TUNIT1 = 'BJD - 2457000, days' / column units: Barycenter corrected TESS JulianTDISP1 = 'D14.7 ' / column display format TTYPE2 = 'TIMECORR' / column title: barycentric correction TFORM2 = 'E ' / column format: 32-bit floating point TUNIT2 = 'd ' / column units: Days TDISP2 = 'E14.7 ' / column display format TTYPE3 = 'CADENCENO' / column title: unique cadence number TFORM3 = 'J ' / column format: signed 32-bit integer TDISP3 = 'I10 ' / column display format TTYPE4 = 'RAW_CNTS' / column title: raw pixel counts TFORM4 = '121J ' / column format: image of signed 32-bit integers TUNIT4 = 'count ' / column units: count TDISP4 = 'I8 ' / column display format TDIM4 = '(11,11) ' / column dimensions: pixel aperture array TNULL4 = -1 / column null value indicator WCSN4P = 'PHYSICAL' / table column WCS name WCAX4P = 2 / table column physical WCS dimensions 1CTY4P = 'RAWX ' / table column physical WCS axis 1 type, CCD col 2CTY4P = 'RAWY ' / table column physical WCS axis 2 type, CCD row 1CUN4P = 'PIXEL ' / table column physical WCS axis 1 unit 2CUN4P = 'PIXEL ' / table column physical WCS axis 2 unit 1CRV4P = 1631 / table column physical WCS ax 1 ref value 2CRV4P = 254 / table column physical WCS ax 2 ref value 1CDL4P = 1.0 / table column physical WCS a1 step 2CDL4P = 1.0 / table column physical WCS a2 step 1CRP4P = 1 / table column physical WCS a1 reference 2CRP4P = 1 / table column physical WCS a2 reference WCAX4 = 2 / number of WCS axes 1CTYP4 = 'RA---TAN' / right ascension coordinate type 2CTYP4 = 'DEC--TAN' / declination coordinate type 1CRPX4 = 5.443758692689471 / [pixel] reference pixel along image axis 1 2CRPX4 = 5.272636267197981 / [pixel] reference pixel along image axis 2 1CRVL4 = 117.0267107898172700 / [deg] right ascension at reference pixel 2CRVL4 = 50.2249536906459700 / [deg] declination at reference pixel 1CUNI4 = 'deg ' / physical unit in column dimension 2CUNI4 = 'deg ' / physical unit in row dimension 1CDLT4 = -0.005613958177205 / [deg] pixel scale in RA dimension 2CDLT4 = 0.005613958177205 / [deg] pixel scale in DEC dimension 11PC4 = 0.9881493277952451 / linear transformation matrix element cos(th) 12PC4 = 0.21324264501531376 / linear transformation matrix element -sin(th) 21PC4 = 0.24022028831990944 / linear transformation matrix element sin(th) 22PC4 = -0.9601532517855615 / linear transformation matrix element cos(th) TTYPE5 = 'FLUX ' / column title: calibrated pixel flux TFORM5 = '121E ' / column format: image of 32-bit floating point TUNIT5 = 'e-/s ' / column units: electrons per second TDISP5 = 'E14.7 ' / column display format TDIM5 = '(11,11) ' / column dimensions: pixel aperture array WCSN5P = 'PHYSICAL' / table column WCS name WCAX5P = 2 / table column physical WCS dimensions 1CTY5P = 'RAWX ' / table column physical WCS axis 1 type, CCD col 2CTY5P = 'RAWY ' / table column physical WCS axis 2 type, CCD row 1CUN5P = 'PIXEL ' / table column physical WCS axis 1 unit 2CUN5P = 'PIXEL ' / table column physical WCS axis 2 unit 1CRV5P = 1631 / table column physical WCS ax 1 ref value 2CRV5P = 254 / table column physical WCS ax 2 ref value 1CDL5P = 1.0 / table column physical WCS a1 step 2CDL5P = 1.0 / table column physical WCS a2 step 1CRP5P = 1 / table column physical WCS a1 reference 2CRP5P = 1 / table column physical WCS a2 reference WCAX5 = 2 / number of WCS axes 1CTYP5 = 'RA---TAN' / right ascension coordinate type 2CTYP5 = 'DEC--TAN' / declination coordinate type 1CRPX5 = 5.443758692689471 / [pixel] reference pixel along image axis 1 2CRPX5 = 5.272636267197981 / [pixel] reference pixel along image axis 2 1CRVL5 = 117.0267107898172700 / [deg] right ascension at reference pixel 2CRVL5 = 50.2249536906459700 / [deg] declination at reference pixel 1CUNI5 = 'deg ' / physical unit in column dimension 2CUNI5 = 'deg ' / physical unit in row dimension 1CDLT5 = -0.005613958177205 / [deg] pixel scale in RA dimension 2CDLT5 = 0.005613958177205 / [deg] pixel scale in DEC dimension 11PC5 = 0.9881493277952451 / linear transformation matrix element cos(th) 12PC5 = 0.21324264501531376 / linear transformation matrix element -sin(th) 21PC5 = 0.24022028831990944 / linear transformation matrix element sin(th) 22PC5 = -0.9601532517855615 / linear transformation matrix element cos(th) TTYPE6 = 'FLUX_ERR' / column title: 1-sigma calibrated uncertainty TFORM6 = '121E ' / column format: image of 32-bit floating point TUNIT6 = 'e-/s ' / column units: electrons per second (1-sigma) TDISP6 = 'E14.7 ' / column display format TDIM6 = '(11,11) ' / column dimensions: pixel aperture array WCSN6P = 'PHYSICAL' / table column WCS name WCAX6P = 2 / table column physical WCS dimensions 1CTY6P = 'RAWX ' / table column physical WCS axis 1 type, CCD col 2CTY6P = 'RAWY ' / table column physical WCS axis 2 type, CCD row 1CUN6P = 'PIXEL ' / table column physical WCS axis 1 unit 2CUN6P = 'PIXEL ' / table column physical WCS axis 2 unit 1CRV6P = 1631 / table column physical WCS ax 1 ref value 2CRV6P = 254 / table column physical WCS ax 2 ref value 1CDL6P = 1.0 / table column physical WCS a1 step 2CDL6P = 1.0 / table column physical WCS a2 step 1CRP6P = 1 / table column physical WCS a1 reference 2CRP6P = 1 / table column physical WCS a2 reference WCAX6 = 2 / number of WCS axes 1CTYP6 = 'RA---TAN' / right ascension coordinate type 2CTYP6 = 'DEC--TAN' / declination coordinate type 1CRPX6 = 5.443758692689471 / [pixel] reference pixel along image axis 1 2CRPX6 = 5.272636267197981 / [pixel] reference pixel along image axis 2 1CRVL6 = 117.0267107898172700 / [deg] right ascension at reference pixel 2CRVL6 = 50.2249536906459700 / [deg] declination at reference pixel 1CUNI6 = 'deg ' / physical unit in column dimension 2CUNI6 = 'deg ' / physical unit in row dimension 1CDLT6 = -0.005613958177205 / [deg] pixel scale in RA dimension 2CDLT6 = 0.005613958177205 / [deg] pixel scale in DEC dimension 11PC6 = 0.9881493277952451 / linear transformation matrix element cos(th) 12PC6 = 0.21324264501531376 / linear transformation matrix element -sin(th) 21PC6 = 0.24022028831990944 / linear transformation matrix element sin(th) 22PC6 = -0.9601532517855615 / linear transformation matrix element cos(th) TTYPE7 = 'FLUX_BKG' / column title: calibrated background flux TFORM7 = '121E ' / column format: image of 32-bit floating point TUNIT7 = 'e-/s ' / column units: electrons per second TDISP7 = 'E14.7 ' / column display format TDIM7 = '(11,11) ' / column dimensions: pixel aperture array WCSN7P = 'PHYSICAL' / table column WCS name WCAX7P = 2 / table column physical WCS dimensions 1CTY7P = 'RAWX ' / table column physical WCS axis 1 type, CCD col 2CTY7P = 'RAWY ' / table column physical WCS axis 2 type, CCD row 1CUN7P = 'PIXEL ' / table column physical WCS axis 1 unit 2CUN7P = 'PIXEL ' / table column physical WCS axis 2 unit 1CRV7P = 1631 / table column physical WCS ax 1 ref value 2CRV7P = 254 / table column physical WCS ax 2 ref value 1CDL7P = 1.0 / table column physical WCS a1 step 2CDL7P = 1.0 / table column physical WCS a2 step 1CRP7P = 1 / table column physical WCS a1 reference 2CRP7P = 1 / table column physical WCS a2 reference WCAX7 = 2 / number of WCS axes 1CTYP7 = 'RA---TAN' / right ascension coordinate type 2CTYP7 = 'DEC--TAN' / declination coordinate type 1CRPX7 = 5.443758692689471 / [pixel] reference pixel along image axis 1 2CRPX7 = 5.272636267197981 / [pixel] reference pixel along image axis 2 1CRVL7 = 117.0267107898172700 / [deg] right ascension at reference pixel 2CRVL7 = 50.2249536906459700 / [deg] declination at reference pixel 1CUNI7 = 'deg ' / physical unit in column dimension 2CUNI7 = 'deg ' / physical unit in row dimension 1CDLT7 = -0.005613958177205 / [deg] pixel scale in RA dimension 2CDLT7 = 0.005613958177205 / [deg] pixel scale in DEC dimension 11PC7 = 0.9881493277952451 / linear transformation matrix element cos(th) 12PC7 = 0.21324264501531376 / linear transformation matrix element -sin(th) 21PC7 = 0.24022028831990944 / linear transformation matrix element sin(th) 22PC7 = -0.9601532517855615 / linear transformation matrix element cos(th) TTYPE8 = 'FLUX_BKG_ERR' / column title: 1-sigma cal. background uncertainTFORM8 = '121E ' / column format: image of 32-bit floating point TUNIT8 = 'e-/s ' / column units: electrons per second (1-sigma) TDISP8 = 'E14.7 ' / column display format TDIM8 = '(11,11) ' / column dimensions: pixel aperture array WCSN8P = 'PHYSICAL' / table column WCS name WCAX8P = 2 / table column physical WCS dimensions 1CTY8P = 'RAWX ' / table column physical WCS axis 1 type, CCD col 2CTY8P = 'RAWY ' / table column physical WCS axis 2 type, CCD row 1CUN8P = 'PIXEL ' / table column physical WCS axis 1 unit 2CUN8P = 'PIXEL ' / table column physical WCS axis 2 unit 1CRV8P = 1631 / table column physical WCS ax 1 ref value 2CRV8P = 254 / table column physical WCS ax 2 ref value 1CDL8P = 1.0 / table column physical WCS a1 step 2CDL8P = 1.0 / table column physical WCS a2 step 1CRP8P = 1 / table column physical WCS a1 reference 2CRP8P = 1 / table column physical WCS a2 reference WCAX8 = 2 / number of WCS axes 1CTYP8 = 'RA---TAN' / right ascension coordinate type 2CTYP8 = 'DEC--TAN' / declination coordinate type 1CRPX8 = 5.443758692689471 / [pixel] reference pixel along image axis 1 2CRPX8 = 5.272636267197981 / [pixel] reference pixel along image axis 2 1CRVL8 = 117.0267107898172700 / [deg] right ascension at reference pixel 2CRVL8 = 50.2249536906459700 / [deg] declination at reference pixel 1CUNI8 = 'deg ' / physical unit in column dimension 2CUNI8 = 'deg ' / physical unit in row dimension 1CDLT8 = -0.005613958177205 / [deg] pixel scale in RA dimension 2CDLT8 = 0.005613958177205 / [deg] pixel scale in DEC dimension 11PC8 = 0.9881493277952451 / linear transformation matrix element cos(th) 12PC8 = 0.21324264501531376 / linear transformation matrix element -sin(th) 21PC8 = 0.24022028831990944 / linear transformation matrix element sin(th) 22PC8 = -0.9601532517855615 / linear transformation matrix element cos(th) TTYPE9 = 'QUALITY ' / column title: pixel quality flags TFORM9 = 'J ' / column format: signed 32-bit integer TDISP9 = 'B16.16 ' / column display format TTYPE10 = 'POS_CORR1' / column title: column position correction TFORM10 = 'E ' / column format: 32-bit floating point TUNIT10 = 'pixel ' / column units: pixel TDISP10 = 'E14.7 ' / column display format TTYPE11 = 'POS_CORR2' / column title: row position correction TFORM11 = 'E ' / column format: 32-bit floating point TUNIT11 = 'pixel ' / column units: pixel TDISP11 = 'E14.7 ' / column display format INHERIT = T / inherit the primary header EXTNAME = 'PIXELS ' / name of extension EXTVER = 1 / extension version number (not format version) SIMDATA = F / file is based on simulated data TELESCOP= 'TESS ' / telescope INSTRUME= 'TESS Photometer' / detector type OBJECT = 'TIC 356473034' / string version of target id TICID = 356473034 / unique tess target identifier RADESYS = 'ICRS ' / reference frame of celestial coordinates RA_OBJ = 117.0267107898172700 / [deg] right ascension DEC_OBJ = 50.2249536906459700 / [deg] declination EQUINOX = 2000.0 / equinox of celestial coordinate system EXPOSURE= 19.630949790673 / [d] time on source TIMEREF = 'SOLARSYSTEM' / barycentric correction applied to times TASSIGN = 'SPACECRAFT' / where time is assigned TIMESYS = 'TDB ' / time system is Barycentric Dynamical Time (TDB)BJDREFI = 2457000 / integer part of BTJD reference date BJDREFF = 0.00000000 / fraction of the day in BTJD reference date TIMEUNIT= 'd ' / time unit for TIME, TSTART and TSTOP TELAPSE = 26.319569811024 / [d] TSTOP - TSTART LIVETIME= 20.8450992903311400 / [d] TELAPSE multiplied by DEADC TSTART = 1842.507957785229 / observation start time in BTJD TSTOP = 1868.827527364771 / observation stop time in BTJD DATE-OBS= '2019-12-25T00:10:18.369' / TSTART as UTC calendar date DATE-END= '2020-01-20T07:50:29.180' / TSTOP as UTC calendar date DEADC = 0.7920000000000000 / deadtime correction TIMEPIXR= 0.5 / bin time beginning=0 middle=0.5 end=1 TIERRELA= 1.16E-05 / [d] relative time error INT_TIME= 1.980000000000 / [s] photon accumulation time per frame READTIME= 0.020000000000 / [s] readout time per frame FRAMETIM= 2.000000000000 / [s] frame time (INT_TIME + READTIME) NUM_FRM = 60 / number of frames per time stamp TIMEDEL = 0.001388888888888889 / [d] time resolution of data BACKAPP = T / background is subtracted DEADAPP = T / deadtime applied VIGNAPP = T / vignetting or collimator correction applied GAINA = 5.21999979019165 / [electrons/count] CCD output A gain GAINB = 5.210000038146973 / [electrons/count] CCD output B gain GAINC = 5.210000038146973 / [electrons/count] CCD output C gain GAIND = 5.260000228881836 / [electrons/count] CCD output D gain READNOIA= 10.022398948669434 / [electrons] read noise CCD output A READNOIB= 7.7108001708984375 / [electrons] read noise CCD output B READNOIC= 7.4502997398376465 / [electrons] read noise CCD output C READNOID= 10.099200248718262 / [electrons] read noise CCD output D NREADOUT= 48 / number of read per cadence FXDOFF = 209700 / compression fixed offset CDPP0_5 = 750.28845215 / RMS CDPP on 0.5-hr time scales CDPP1_0 = 572.43994141 / RMS CDPP on 1.0-hr time scales CDPP2_0 = 456.91378784 / RMS CDPP on 2.0-hr time scales CROWDSAP= 0.83436364 / Ratio of target flux to total flux in op. ap. FLFRCSAP= 0.71524101 / Frac. of target flux w/in the op. aperture CHECKSUM= '3a9nAT9l5Z9lAZ9l' / HDU checksum updated 2020-03-24T02:36:06Z TMOFST11= 0.0 / (s) readout delay for camera 1 and ccd 1 MEANBLCA= 6653 / [count] FSW mean black level CCD output A MEANBLCB= 6720 / [count] FSW mean black level CCD output B MEANBLCC= 6612 / [count] FSW mean black level CCD output C MEANBLCD= 6458 / [count] FSW mean black level CCD output D END 11 11
可以看到这个 HDU 的头部（Header）中包含了 WCS 的信息。接下来将Header中与WCS相关的关键字整理成一个 dict 并初始化一个WCS对象
wcs_input_dict = {
'CTYPE1': head1['1CTYP5'],
'CUNIT1': head1['1CUNI5'],
'CDELT1': head1['1CDLT5'],
'CRPIX1': head1['1CRPX5'],
'CRVAL1': head1['1CRVL5'],
'NAXIS1': nx,
'CTYPE2': head1['2CTYP5'],
'CUNIT2': head1['2CUNI5'],
'CDELT2': head1['2CDLT5'],
'CRPIX2': head1['2CRPX5'],
'CRVAL2': head1['2CRVL5'],
'NAXIS2': ny,
'PC1_1': head1['11PC5'],
'PC1_2': head1['12PC5'],
'PC2_1': head1['21PC5'],
'PC2_2': head1['22PC5'],
}
wcoord = WCS(wcs_input_dict)
利用这个 WCS 对象重新绘制上面的图
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(8,8))
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection=wcoord)
ax.imshow(img)
<matplotlib.image.AxesImage at 0x7f2765bd3b10>
接下来对上图进行适当的美化，并使用一个蓝色为主基调的颜色表
fig = plt.figure(figsize=(8,8))
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection=wcoord)
ax.imshow(img,cmap='YlGnBu_r')
xcoords = ax.coords[0]
ycoords = ax.coords[1]
xcoords.set_major_formatter('d.ddd')
ycoords.set_major_formatter('d.ddd')
xcoords.ticklabels.set_fontsize(12)
ycoords.ticklabels.set_fontsize(12)
xcoords.set_axislabel('RA (deg)',fontsize=12)
ycoords.set_axislabel('Dec (deg)',fontsize=12)
ax.grid(True, color='w', ls='--', lw=0.5)
plt.show()
TESS 的 TP 文件中也包含了抽取光变曲线所采用的孔径信息，在第三个HDU 中，用一个二维掩板（MASK）的形式存储
data2 = hdulst[2].data
print(data2)
[[257 261 261 261 257 257 257 257 257 257 257] [257 261 261 261 257 257 257 261 261 261 257] [257 257 261 257 257 257 257 257 261 261 257] [257 257 257 257 267 267 267 257 261 261 257] [261 261 257 267 267 267 267 257 261 261 257] [261 261 257 257 257 267 267 257 261 261 257] [261 261 257 257 257 257 257 257 261 257 257] [261 261 257 257 257 257 257 261 257 257 257] [257 257 257 257 257 257 261 261 257 257 257] [257 257 257 257 257 257 257 257 257 257 257] [257 257 257 257 257 257 257 257 257 257 257]]
将这个数组与 2 做“与”运算即可得到测光孔径信息
aperture = data2&2>0
print(aperture)
[[False False False False False False False False False False False] [False False False False False False False False False False False] [False False False False False False False False False False False] [False False False False True True True False False False False] [False False False True True True True False False False False] [False False False False False True True False False False False] [False False False False False False False False False False False] [False False False False False False False False False False False] [False False False False False False False False False False False] [False False False False False False False False False False False] [False False False False False False False False False False False]]
接下来我们定义一个函数将上面这个Mask中标记为 True 的像素转换为像素的边界
import numpy as np
def get_aperture_bound(aperture):
bound_lst = []
ny, nx = aperture.shape
for iy in np.arange(ny):
for ix in np.arange(nx):
if aperture[iy, ix]>0:
line_lst = [(ix, iy, ix, iy+1),
(ix, iy+1, ix+1, iy+1),
(ix+1, iy, ix+1, iy+1),
(ix, iy, ix+1, iy),
]
for line in line_lst:
if line in bound_lst:
index = bound_lst.index(line)
bound_lst.pop(index)
else:
bound_lst.append(line)
return bound_lst
重新绘图，并用红色画出孔径信息
fig = plt.figure(figsize=(8,8))
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection=wcoord)
ax.imshow(img, cmap='YlGnBu_r')
bound_lst = get_aperture_bound(aperture)
for (x1, y1, x2, y2) in bound_lst:
ax.plot([x1-0.5, x2-0.5], [y1-0.5, y2-0.5], 'r-', lw=1)
xcoords = ax.coords[0]
ycoords = ax.coords[1]
xcoords.set_major_formatter('d.ddd')
ycoords.set_major_formatter('d.ddd')
xcoords.ticklabels.set_fontsize(7)
ycoords.ticklabels.set_fontsize(7)
xcoords.set_axislabel('RA (deg)',fontsize=7)
ycoords.set_axislabel('Dec (deg)',fontsize=7)
ax.grid(True, color='w', ls='--', lw=0.5)
因为TESS望远镜的每像素对应的天空张角很大（21"/pixel），我们希望得到目标星周围的其他星的信息，可以用 astroquery包的 Vizier 和 Simbad 类对周围星的情况做查询。
导入软件包：
from astroquery.vizier import Vizier
from astroquery.simbad import Simbad
首先查询目标星的RA、Dec 坐标
Simbad.query_object('XO-2N')
| MAIN_ID | RA | DEC | RA_PREC | DEC_PREC | COO_ERR_MAJA | COO_ERR_MINA | COO_ERR_ANGLE | COO_QUAL | COO_WAVELENGTH | COO_BIBCODE | SCRIPT_NUMBER_ID |
|---|---|---|---|---|---|---|---|---|---|---|---|
| "h:m:s" | "d:m:s" | mas | mas | deg | |||||||
| object | str13 | str13 | int16 | int16 | float32 | float32 | int16 | str1 | str1 | object | int32 |
| BD+50 1471 | 07 48 06.4723 | +50 13 32.920 | 14 | 14 | 0.012 | 0.011 | 90 | A | O | 2020yCat.1350....0G | 1 |
可见，返回一个目标星的表格，其中的 RA 和 DEC 列就是我们想要知道的赤道坐标。
result = Simbad.query_object('XO-2N')
ra = result[0]['RA']
dec = result[0]['DEC']
print(ra, dec)
07 48 06.4723 +50 13 32.920
接下来将Simbad得到的RA、Dec 坐标转换成一个 Astropy 的SkyCoord对象，然后在 TIC (TESS 输入星表)中查询目标星周围 120" 的天体
from astropy.coordinates import SkyCoord
import astropy.units as u
coord = SkyCoord(ra, dec, unit=(u.hourangle, u.deg))
tablelist = Vizier(catalog='IV/39/tic82', columns=['**','_+r']).query_region(coord, radius=120*u.arcsec)
table = tablelist['IV/39/tic82']
print(table)
_r TIC RAJ2000 DEJ2000 ... WDFl ID Sloan
deg deg ...
------- --------- --------------- --------------- ... ---- --------- -----
111.148 356473016 117.04982487654 50.19861339785 ... -1 129540237 Sloan
82.149 741714554 117.04379109132 50.20568263642 ... 0 129625485 Sloan
114.334 356473015 117.03638312960 50.19462167800 ... 0 129540236 Sloan
70.694 356473021 117.02109316838 50.20653178336 ... 0 129540238 Sloan
0.025 356473034 117.02696916943 50.22581142130 ... 0 129540244 Sloan
31.178 356473029 117.03117259168 50.21757161613 ... 0 129540240 Sloan
... ... ... ... ... ... ... ...
99.045 356473042 117.03328117521 50.25302012504 ... 0 129540245 Sloan
105.823 356473030 116.98275276630 50.21783201486 ... 0 129248114 Sloan
101.780 356473036 116.98411762257 50.23270998683 ... 0 129248116 Sloan
103.402 356473035 116.98245961599 50.22954991993 ... 0 129248115 Sloan
55.771 742731700 117.00275893505 50.22581463622 ... 0 129625874 Sloan
118.886 742731701 116.97814698495 50.23652849537 ... 0 129333414 Sloan
76.999 742731699 116.99586975530 50.21796739032 ... 1 129333413 Sloan
Length = 21 rows
WARNING: UnitsWarning: Unit 'Sun' not supported by the VOUnit standard. Did you mean uN? [astropy.units.format.vounit]
这个表格包含的所有列如下：
print(table.colnames)
['_r', 'TIC', 'RAJ2000', 'DEJ2000', 'HIP', 'TYC', 'UCAC4', '_2MASS', 'objID', 'WISEA', 'GAIA', 'APASS', 'KIC', 'S_G', 'Ref', 'r_Pos', 'pmRA', 'e_pmRA', 'pmDE', 'e_pmDE', 'r_pm', 'Plx', 'e_Plx', 'r_Plx', 'GLON', 'GLAT', 'ELON', 'ELAT', 'Bmag', 'e_Bmag', 'Vmag', 'e_Vmag', 'umag', 'e_umag', 'gmag', 'e_gmag', 'rmag', 'e_rmag', 'imag', 'e_imag', 'zmag', 'e_zmag', 'Jmag', 'e_Jmag', 'Hmag', 'e_Hmag', 'Kmag', 'e_Kmag', 'q_2MASS', 'W1mag', 'e_W1mag', 'W2mag', 'e_W2mag', 'W3mag', 'e_W3mag', 'W4mag', 'e_W4mag', 'Gmag', 'e_Gmag', 'Tmag', 'e_Tmag', 'f_Tmag', 'Flag', 'Teff', 's_Teff', 'logg', 's_logg', '__M_H_', 'e__M_H_', 'Rad', 's_Rad', 'Mass', 's_Mass', 'rho', 's_rho', 'LClass', 'Lum', 's_Lum', 'Dist', 's_Dist', 'E_B-V_', 's_E_B-V_', 'Ncont', 'Rcont', 'Disp', 'm_TIC', 'Prior', 'e_E_B-V_', 'E_E_B-V_', 'f_E_B-V_', 'e_Mass', 'E_Mass', 'e_Rad', 'E_Rad', 'e_rho', 'E_rho', 'e_logg', 'E_logg', 'e_Lum', 'E_Lum', 'e_Dist', 'E_Dist', 'r_Dist', 'e_Teff', 'E_Teff', 'r_Teff', 'BPmag', 'e_BPmag', 'RPmag', 'e_RPmag', 'q_Gaia', 'r_Vmag', 'r_Bmag', 'Clist', 'e_RAJ2000', 'e_DEJ2000', 'RAOdeg', 'DEOdeg', 'e_RAOdeg', 'e_DEOdeg', 'RadFl', 'WDFl', 'ID', 'Sloan']
我们从中选取RA、Dec、Tmag 这三列画图，把每一颗星用圆圈标记在刚才的二维图像上，并且用大小表示星等。目标星的星等越亮，圆圈越大。
ra_lst = table['RAJ2000']
dec_lst = table['DEJ2000']
tmag_lst = table['Tmag']
coord_lst = SkyCoord(ra_lst, dec_lst, unit='deg')
用 wcs对象的 all_world2pix 函数把ra, dec坐标转换成图像上的x、y值
fig = plt.figure(figsize=(8,8))
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection=wcoord)
ax.imshow(img, cmap='YlGnBu_r')
bound_lst = get_aperture_bound(aperture)
for (x1, y1, x2, y2) in bound_lst:
ax.plot([x1-0.5, x2-0.5], [y1-0.5, y2-0.5], 'r-', lw=1)
size = np.maximum((22-tmag_lst)*10, 0.1)
x_lst, y_lst = wcoord.all_world2pix(ra_lst, dec_lst, 0)
ax.scatter(x_lst, y_lst, s=size, c='none', ec='r', lw=1)
xcoords = ax.coords[0]
ycoords = ax.coords[1]
xcoords.set_major_formatter('d.ddd')
ycoords.set_major_formatter('d.ddd')
xcoords.ticklabels.set_fontsize(12)
ycoords.ticklabels.set_fontsize(12)
xcoords.set_axislabel('RA (deg)',fontsize=12)
ycoords.set_axislabel('Dec (deg)',fontsize=12)
ax.grid(True, color='w', ls='--', lw=0.5)
plt.show()
由于我们取的星表的查询半径比图像的尺寸略大，为了兼顾所有的星，matplotlib自动把x、y轴的显示范围进行了扩大。可以在画星点之前取的x、y轴的显示范围，待画完后重置x、y轴。
fig = plt.figure(figsize=(8,8))
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection=wcoord)
ax.imshow(img, cmap='YlGnBu_r')
_x1, _x2 = ax.get_xlim()
_y1, _y2 = ax.get_ylim()
bound_lst = get_aperture_bound(aperture)
for (x1, y1, x2, y2) in bound_lst:
ax.plot([x1-0.5, x2-0.5], [y1-0.5, y2-0.5], 'r-', lw=1)
size = np.maximum((22-tmag_lst)*10, 0.1)
x_lst, y_lst = wcoord.all_world2pix(ra_lst, dec_lst, 0)
ax.scatter(x_lst, y_lst, s=size, c='none', ec='r', lw=1)
ax.set_xlim(_x1, _x2)
ax.set_ylim(_y1, _y2)
xcoords = ax.coords[0]
ycoords = ax.coords[1]
xcoords.set_major_formatter('d.ddd')
ycoords.set_major_formatter('d.ddd')
xcoords.ticklabels.set_fontsize(12)
ycoords.ticklabels.set_fontsize(12)
xcoords.set_axislabel('RA (deg)',fontsize=12)
ycoords.set_axislabel('Dec (deg)',fontsize=12)
ax.grid(True, color='w', ls='--', lw=0.5)
plt.show()
由于TESS图像的每像素张角太大，我们也可以使用Astroquery提供的SkyView接口来查询其他望远镜在同一天区拍摄的图像
导入SkyView类：
from astroquery.skyview import SkyView
查询目标星周围360角秒的DSS 巡天图像
paths = SkyView.get_images(position=coord,
survey='DSS',
radius = 360*u.arcsec,
sampler ='Clip',
scaling = 'Log',
pixels = (500, 500),
)
获取巡天图像：
hdulst = paths[0]
hdu = hdulst[0]
dssdata = hdu.data
dsshead = hdu.header
DSS数据的Header中也包含了一个wcs坐标信息，可以帮助我们为DSS图像数据建立赤道坐标系下的投影
wcoord2 = WCS(dsshead)
显示该区域图像：
fig2 = plt.figure(figsize=(8,8))
ax = fig2.add_axes([0.1, 0.1, 0.8, 0.8], projection=wcoord2)
ax.imshow(dssdata)
plt.show()
接下来我们查询这个天区中目标星的Gaia 3数据，并且把自行的大小和方向用箭头画在上面的图上
catid = 'I/355/gaiadr3'
viz = Vizier(catalog=catid, columns=['**','+_r'])
viz.ROW_LIMIT=-1
tablelist = viz.query_region(coord, radius=360*u.arcsec)
table = tablelist[catid]
print(table)
_r DR3Name RA_ICRS ... e_DEJ2000 RADEcorJ2000
deg ... mas
------- --------------------------- --------------- ... --------- ------------
2.493 Gaia DR3 934346809278715776 117.02676263913 ... 0.183818 0.2070
33.403 Gaia DR3 934346740559239296 117.03096897417 ... 0.201551 0.1990
43.682 Gaia DR3 934346774918979328 117.01127879980 ... 1.188749 0.1890
48.436 Gaia DR3 934347114221102976 117.04712252495 ... 9.507419 -0.0052
55.754 Gaia DR3 982385228210176896 117.00276630629 ... 2.711706 0.0046
70.770 Gaia DR3 934346706199502464 117.02109606239 ... 2.873963 0.0576
... ... ... ... ... ...
354.111 Gaia DR3 934345499313725824 117.09400929799 ... 3.857020 0.0742
354.181 Gaia DR3 982384541015414912 116.87340531103 ... 0.624162 0.1779
354.455 Gaia DR3 934352306836841472 117.16180131591 ... 1.522227 0.1555
355.106 Gaia DR3 982384747173843840 116.87858336632 ... 0.473481 0.1226
355.299 Gaia DR3 934344060499650432 117.04913547891 ... 1.240011 0.1160
356.264 Gaia DR3 982391790920203776 117.05783498617 ... 3.915965 0.1323
358.248 Gaia DR3 982383097906403584 116.91122197673 ... 0.876622 0.1121
Length = 201 rows
注意在上面这个例子中将 ROW_LIMIT （返回的表格最大行数）设置为 -1，表示不设置行数上限。
查看Gaia DR3星表的列信息：
print(table.colnames)
['_r', 'DR3Name', 'RA_ICRS', 'DE_ICRS', 'SolID', 'Source', 'RandomI', 'e_RA_ICRS', 'e_DE_ICRS', 'Plx', 'e_Plx', 'RPlx', 'PM', 'pmRA', 'e_pmRA', 'pmDE', 'e_pmDE', 'RADEcor', 'RAPlxcor', 'RApmRAcor', 'RApmDEcor', 'DEPlxcor', 'DEpmRAcor', 'DEpmDEcor', 'PlxpmRAcor', 'PlxpmDEcor', 'pmRApmDEcor', 'NAL', 'NAC', 'NgAL', 'NbAL', 'gofAL', 'chi2AL', 'epsi', 'sepsi', 'Solved', 'APF', 'nueff', 'pscol', 'e_pscol', 'RApscolCorr', 'DEpscolCorr', 'PlxpscolCorr', 'pmRApscolCorr', 'pmDEpscolCorr', 'MatchObsA', 'Nper', 'amax', 'MatchObs', 'NewMatchObs', 'MatchObsrm', 'IPDgofha', 'IPDgofhp', 'IPDfmp', 'IPDfow', 'RUWE', 'SDSk1', 'SDSk2', 'SDSk3', 'SDSk4', 'SDMk1', 'SDMk2', 'SDMk3', 'SDMk4', 'Dup', 'o_Gmag', 'FG', 'e_FG', 'RFG', 'Gmag', 'e_Gmag', 'o_BPmag', 'FBP', 'e_FBP', 'RFBP', 'BPmag', 'e_BPmag', 'o_RPmag', 'FRP', 'e_FRP', 'RFRP', 'RPmag', 'e_RPmag', 'E_BP_RP_', 'NBPcont', 'NBPblend', 'NRPcont', 'NRPblend', 'Mode', 'BP-RP', 'BP-G', 'G-RP', 'RV', 'e_RV', 'n_RV', 'o_RV', 'o_RVd', 'RVNper', 'RVS_N', 'RVgof', 'RVchi2', 'RVTdur', 'RVamp', 'RVtempTeff', 'RVtemplogg', 'RVtemp_Fe_H_', 'Vatmparam', 'Vbroad', 'e_Vbroad', 'o_Vbroad', 'GRVSmag', 'e_GRVSmag', 'o_GRVSmag', 'RVSS_N', 'VarFlag', 'GLON', 'GLAT', 'ELON', 'ELAT', 'QSO', 'Gal', 'NSS', 'XPcont', 'XPsamp', 'RVS', 'EpochPh', 'EpochRV', 'MCMCGSP', 'MCMCMSC', 'And', 'PQSO', 'PGal', 'PSS', 'Teff', 'b_Teff', 'B_Teff', 'logg', 'b_logg', 'B_logg', '__Fe_H_', 'b__Fe_H_', 'B__Fe_H_', 'Dist', 'b_Dist', 'B_Dist', 'A0', 'b_A0', 'B_A0', 'AG', 'b_AG', 'B_AG', 'E_BP-RP_', 'b_E_BP-RP_', 'B_E_BP-RP_', 'Lib', 'HIP', 'dHIP', 'nHIP', 'f_HIP', 'PS1coid', 'PS1', 'dPS1', 'nPS1', 'mPS1', 'f_PS1', 'SDSS13coid', 'SDSS13', 'dSDSS13', 'nSDSS13', 'mSDSS13', 'f_SDSS13', 'SKYM2', 'dSKYM2', 'nSKYM2', 'mSKYM2', 'f_SKYM2', 'TYC2', 'dTYC2', 'f_TYC2', 'TYC2moid', 'nTYC2', 'URAT1', 'dURAT1', 'f_URAT1', 'URAT1oid', 'nURAT1', 'mURAT1', 'AllWISE', 'dAllWISE', 'f_AllWISE', 'AllWISEoid', 'nAllWISE', 'mAllWISE', 'APASS9coid', 'APASS9', 'dAPASS9', 'nAPASS9', 'mAPASS9', 'f_APASS9', 'GSC23', 'dGSC23', 'f_GSC23', 'GSC23coid', 'nGSC23', 'mGSC23', 'RAVE5', 'dRAVE5', 'f_RAVE5', 'RAVE5coid', 'nRAVE5', '_2MASS', 'd2MASS', 'f_2MASS', '_2MASScoid', 'n2MASS', 'm2MASS', 'RAVE6', 'dRAVE6', 'f_RAVE6', 'RAVE6oid', 'nRAVE6', 'RAJ2000', 'DEJ2000', 'e_RAJ2000', 'e_DEJ2000', 'RADEcorJ2000']
我们需要其中的RA、DEC坐标，以及沿着两个方向的自行信息。我们把这几个数据从表格中拿出来，转换成 Numpy的数组，并把其中的空值(masked)变成0.0
ra_lst = list(table['RA_ICRS'])
dec_lst = list(table['DE_ICRS'])
gmag_lst = list(table['Gmag'])
pm_ra = list(table['pmRA'])
pm_dec = list(table['pmDE'])
pm_ra = np.array([0.0 if v is np.ma.masked else v for v in pm_ra])
pm_dec = np.array([0.0 if v is np.ma.masked else v for v in pm_dec])
接下来设置一个适当的时间间隔，计算这个时间后，目标星移动到的新位置（近似）
pmlen = 1000
d_ra = pm_ra*1e-3/3600./np.cos(np.deg2rad(dec_lst))*pmlen
d_dec = pm_dec*1e-3/3600.*pmlen
ra2_lst = ra_lst + d_ra
dec2_lst = dec_lst + d_dec
类似的，需要把现在的坐标、新的位置坐标转换成图像上的x、y位置
x_lst, y_lst = wcoord2.all_world2pix(ra_lst, dec_lst, 0)
x2_lst, y2_lst = wcoord2.all_world2pix(ra2_lst, dec2_lst, 0)
重新画图，这次使用一个新的颜色表，用灰度表示。
fig2 = plt.figure(figsize=(8,8))
ax = fig2.add_axes([0.1, 0.1, 0.8, 0.8], projection=wcoord2)
ax.imshow(dssdata, cmap='gray_r')
_x1, _x2 = ax.get_xlim()
_y1, _y2 = ax.get_ylim()
size = np.maximum((20-np.array(gmag_lst))*10, 0.1)
ax.scatter(x_lst, y_lst, s=size, c='none', ec='r', lw=0.5)
dx_lst = x2_lst - x_lst
dy_lst = y2_lst - y_lst
ax.quiver(x_lst, y_lst, dx_lst, dy_lst, width=1, units='dots',
angles='xy', scale_units='xy', scale=1, color='r')
ax.set_xlim(_x1, _x2)
ax.set_ylim(_y1, _y2)
plt.show()
也可以用matplotlib的anmiation模块绘制视频
import matplotlib.animation as animation
fig = plt.figure(figsize=(14, 14))
def update(i):
fig.clf()
ax1 = fig.add_axes([0.08, 0.10, 0.40, 0.80],
projection=wcoord)
cax = ax1.imshow(data['FLUX'][i],cmap='YlGnBu_r')
_x1, _x2 = ax1.get_xlim()
_y1, _y2 = ax1.get_ylim()
# plot aperture
for (x1, y1, x2, y2) in bound_lst:
ax1.plot([x1-0.5, x2-0.5], [y1-0.5, y2-0.5], 'r-')
# adjust ax1
ax1.text(0.9*_x1+0.1*_x2, 0.1*_y1+0.9*_y2,
't={:9.4f}'.format(data['TIME'][i]),
color='w')
ax1.set_xlim(_x1, _x2)
ax1.set_ylim(_y1, _y2)
xcoords = ax1.coords[0]
ycoords = ax1.coords[1]
xcoords.set_major_formatter('d.ddd')
ycoords.set_major_formatter('d.ddd')
xcoords.set_axislabel('RA (deg)')
ycoords.set_axislabel('Dec (deg)')
return cax,
anim = animation.FuncAnimation(fig, update,
frames=np.arange(len(data)), interval=1, blit=False)
anim.save('XO-2.mp4', fps=25, extra_args=['-vcodec', 'libx264'])