可以利用astropy.coordinates.SkyCoord直接制定观测计划，首先需要构造源，既可以根据名称让astropy自动查找，也能通过ra和dec来构造。astropy采用CDS name resolver来通过名称获取源的信息，查看CDS (Centre de Données astronomiques de Strasbourg)获取更多信息。源被构造出来之后可以转入其他坐标系。两种方法以及坐标系转换如下所示：
from astropy.time import Time
from astropy.coordinates import SkyCoord, ICRS
import numpy as np
import matplotlib.pyplot as plt
src1 = SkyCoord.from_name('vega')
ra = src1.ra.value # in deg
dec = src1.dec.value # in deg
src2 = SkyCoord(ra=ra, dec=dec, unit='deg', frame=ICRS) # always use ICRS
src2.name = 'mysrc'
print(src1)
print(src2)
src3 = src2.transform_to('galactic')
print(src3)
输出是：
更加详细的文档查看astropy.coordinates.SkyCoord的文档。
构造出源之后可以利用astropy.coordinates.EarthLocation构造观测站位置：
from astropy.coordinates import EarthLocation observer = EarthLocation(lat=39.9042*units.deg, lon=116.4074*units.deg, height=55*units.m) print(observer)
具体参考astropy.coordinates.EarthLocation的文档。
之后可以利用astropy.coordinates.AltAz获取源的az和alt，具体而言有两种做法：
# first
tobs = Time.now()
src1.obstime= tobs
src1.location = observer
print('altaz: ', src1.altaz)
print('az: ', src1.altaz.az)
print('alt: ', src1.altaz.alt)
# second
from astropy.coordinates import AltAz
srcpos = src2.transform_to(AltAz(obstime=tobs, location=observer))
print('az: ', srcpos.az)
print('alt: ', srcpos.alt)
输出为：
altaz:az: 84d21m24.2372s alt: 69d16m23.8372s az: 84d21m24.2372s alt: 69d16m23.8372s
完整代码如下：
# initial source
from astropy.time import Time
from astropy.coordinates import SkyCoord, ICRS
import numpy as np
import matplotlib.pyplot as plt
src1 = SkyCoord.from_name('vega')
ra = src1.ra.value # in deg
dec = src1.dec.value # in deg
src2 = SkyCoord(ra=ra, dec=dec, unit='deg', frame=ICRS) # always use ICRS
src2.name = 'mysrc'
print(src1)
print(src2)
src3 = src2.transform_to('galactic')
print(src3)
# initial observer
from astropy.coordinates import EarthLocation
observer = EarthLocation(lat=39.9042*units.deg, lon=116.4074*units.deg, height=55*units.m)
print(observer)
# calculate az and alt
# first
tobs = Time.now()
src1.obstime= tobs
src1.location = observer
print('altaz: ', src1.altaz)
print('az: ', src1.altaz.az)
print('alt: ', src1.altaz.alt)
# second
from astropy.coordinates import AltAz
srcpos = src2.transform_to(AltAz(obstime=tobs, location=observer))
print('az: ', srcpos.az)
print('alt: ', srcpos.alt)
另一种方法是使用astroplan这个包，这个包提供了更加多样化的功能，并且可以用以上提到的EarthLocation和SkyCoord对象进行初始化，详见astroplan.Observer的文档。针对于天籁射电阵列，我写了一个程序如下，可用于参考（使用 Python2.7编写，Python3.8也可以运行成功）：
import numpy as np
import sys
import matplotlib.pyplot as plt
from astropy.coordinates import SkyCoord, ICRS
from astropy import units
from astroplan import Observer, FixedTarget
from astropy.time import Time, TimeDelta
from astropy.utils.data import download_file
from astropy.utils import iers
# download from mirror, if downloaded, it will not be download again
# in fact this correction is very small, ~1.e-1 seconds
iers.conf.iers_auto_url = 'https://datacenter.iers.org/data/9/finals2000A.all'
# iers.IERS_A_URL = 'https://datacenter.iers.org/data/9/finals2000A.all'
# from astroplan import download_IERS_A
# download_IERS_A()
def load_srctxt(fname):
'''Load src.txt file.
Return source list.
'''
src_dict = {}
with open(fname, 'r') as f:
for line in f:
line = line.strip()
if line.startswith('#') or line == '':
continue
# print line
l = line.split('t')
#s, ra, dec, alt, flux = line.split('t')
s, ra, dec, flux = line.split('t')
s = s.replace(' ', '_')
#src_dict[s] = map(float, (ra, dec, alt, flux))
src_dict[s] = map(float, [ra, dec, flux])
return src_dict
def _source(srcname=None, ra=None, dec=None, unit='deg', frame=ICRS, **kwargs):
assert (srcname is not None) or (ra is not None and dec is not None), 'Must input srcname or ra and dec!'
if srcname is not None:
try:
print('Searching for source %s!'%srcname)
source = SkyCoord.from_name(srcname, frame=frame)
print('Loaded source %s!'%srcname)
return source
except Exception as e:
print(e)
print('Cannot load source %s! Try to use ra and dec to build it!'%srcname)
if srcname is not None:
assert ra is not None and dec is not None, 'Cannot build because ra or dec is None!'
src = SkyCoord(ra=ra, dec=dec, frame=frame, unit=unit, **kwargs)
source = FixedTarget(src, name=srcname)
return source
def get_src(srcname, srcfile='tlskysrc.txt', astropy_prior=False, unit='deg', frame=ICRS, **kwargs):
if srcfile is not None:
srcdict = load_srctxt(srcfile)
else:
srcdict = {}
if srcname in srcdict:
ra, dec, _ = srcdict[srcname]
if not astropy_prior:
print('Build source %s with information in %s!'%(srcname, srcfile))
src = SkyCoord(ra=ra, dec=dec, frame=frame, unit=unit, **kwargs)
return FixedTarget(src, name=srcname)
else:
print('Source %s is not in %s!'%(srcname, srcfile))
ra = dec = None
return _source(srcname, ra=ra, dec=dec, frame=frame, unit=unit, **kwargs)
def tianlai():
#timezone = TimezoneInfo(utc_offset=8*units.hour)
#timezone = 'UTC'
timezone = 'Asia/Shanghai'
obs = Observer(latitude='44d9m9.66s', longitude='91d48m24.72s', elevation=1493.7*units.m, pressure=0, name='TianLai Array', timezone=timezone)
return obs
#from check_beam import get_azalt
#
#srcname = 'Cyg_A'
srcname = 'Cas_A'
src = get_src(srcname)
obs = tianlai()
tt = obs.target_meridian_transit_time(Time.now(), src, which='next')
##tt = Time.now()
print(obs.altaz(tt, src))
#lon = obs.location.lon.to('deg').value
#lat = obs.location.lat.to('deg').value
#dt = tt.to_datetime(timezone=obs.timezone)
dt = obs.astropy_time_to_datetime(tt)
print(dt)
注：iers.conf.iers_auto_url = 'https://datacenter.iers.org/data/9/finals2000A.all'是因为原来的链接不能下载，所以我换了个镜像，如果这个也不行，可以自行更换。
其中tlskysrc.txt为：
#Defined by Jixia Li on 2019/11/06 #Maybe updated in the future. #Object_Name RA DEC Flux750 NorthPole 0.00 90.0 0 3C_010 6.3083 64.14431 62 3C_058 31.39683 64.82633 34 3C_123 69.26823 29.6705 76 NRAO_206 80.65927 33.46381 39.4 3C_147.1 85.4255 -1.89848 54 3C_147 85.65057 49.85201 34 IC_443 94.15586 22.53166 190 Hydra_A 139.52362 -12.09554 81 3C_273 187.27792 2.05239 47 M_087 187.70593 12.39112 353 3C_295 212.8355 52.20277 37 Hercules_A 252.78395 4.99259 88 3C_353 260.11733 -0.97962 88 NRAO_569 278.84711 -7.32547 80 NRAO_576 280.31136 -5.13509 61 3C_392 284.04304 1.31606 242 3C_397 286.8871 7.14251 41 NRAO_598 287.57728 9.06852 45 3C_398 287.7722 9.09415 45 3C_400 290.74167 14.19729 673 3C_403.2 298.56181 32.89873 75 Cyg_A 299.86815 40.73392 2980 NRAO_621 300.43736 33.29001 53 NRAO_650 318.08928 52.39816 48 Cas_A 350.866417 58.811778 4858 M_001 83.633212 22.01446 1000 3C_061.1 35.646024 86.318380 8.3 PSR_B0329+54(*) 53.248 54.58 1.500 PSR_B0531+21(*) 83.63 22.01 0.55 #PSR_J0534+22(*) 83.63 22.01 0.55 #PSR_J0332+5434 53.0 54.56 0.71452 PSR_J0341+5711 55.25 57.18 0.3647 PSR_B0950+08 148.288 7.93 0.25307 PSR_B1133+16 174.013 15.85 0.257 PSR_B1642-03 251.258 -3.30 0.38769 PSR_B1929+10 293.058 10.99 0.303 PSR_B1933+16 293.949 16.28 0.242 PSR_B1937+21 294.908 21.58 0.240 PSR_B2016+28 304.517 28.67 0.314 PSR_B2111+46 318.350 46.74 0.230 #PSR_B1749-28 268.246 -28.11 0.56256 #NGC_5907 228.97 56.33 0 SGR_1935+2154 293.732 21.897 0 FRB201124 76.98475 26.07726 0 #src0 293.73 21.90 1000 #this_one 281.63 -2.99 1.0
初始化观测者的时候可以输入观测者的时区，建议使用字符串初始化，可用的时区对应的字符串可以用以下代码查看：
from pytz import all_timezones
for tz in all_timezones:
print(tz)
