利用astropy和astroplan制定观测计划

利用astropy和astroplan制定观测计划

可以利用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)

输出是：

<SkyCoord (ICRS): (ra, dec) in deg
    (279.23473479, 38.78368896)>
<SkyCoord (ICRS): (ra, dec) in deg
    (279.23473479, 38.78368896)>
<SkyCoord (Galactic): (l, b) in deg
    (67.44820814, 19.23725227)>

更加详细的文档查看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:  <SkyCoord (AltAz: obstime=2021-11-08 05:53:57.996401, location=(-2179092.56684014, 4388330.31647065, 4069866.70528184) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron): (az, alt) in deg
    (84.35673256, 69.2732881)>
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)