和地理学差不多,天文的研究对象通常也可看成分布在一个球面上,此球面即为天球。天球的半径并不重要,定为多少其实都没有关系,教科书上说是无穷大,也并非不可,但是此说法通常误导人。在描述天体的时候,如果不需要考虑距离,就把它们投影到天球上(就是取其球坐标中的俯仰角和方位角,而不考虑半径,这才是本质)。
通常我们的书和电脑屏幕都是平面,我们经常需要在上面表示整个世界或者整个天球。不过,即使从拓扑学的角来讲,平面(特别是有限大小的平面区域)和球面也有本质的不同。所以如果要在一个平面区域上表示出球面的一部分或者整个球面就需要一些变换,就是投影。
关于天文里常用的FITS文件中对投影的规范,可以参考http://www.cv.nrao.edu/fits/documents/wcs/wcs.html,里面有一篇文章专门叙述了FITS文件所接受的投影变换的细节。
这里想记的其实只是怎么把一个包含了全体中性氢数据的文件里的数据画出来,并标上坐标。首先要做的是读出数据,目前比较好的一个全体的中性氢分布数据是LAB巡天的数据。假设文件名为lab.fit,那么首先可以用mrdfits读出
IDL>a=mrdfits(‘lab.fit’)
这时a是一个三维数组,为了画全天的分布图,可以画一个频率通道,也可以画总强度。假设画总强度
IDL>b=total(a,3)
下面就是画图了
IDL>LoadCT, 39, NColors=colors, Bottom=1, /Silent
IDL>device,decompose=0
这是定义颜色,然后对数据作投影,画坐标线
IDL>map_set,0,0,/MOLLWEIDE,/ISOTROPIC,/HORIZON,/GRID
IDL>result=map_image(b,Startx,Starty,Xsize,Ysize,compress=1,LATMIN=latmin,LONMIN=lonmin,LATMAX=latmax,LONMAX=lonmax,scale=0.1)
IDL>result=bytscl(result)
IDL>tvscl,result,Startx,Starty,XSIZE=Xsize,YSIZE=Ysize
IDL>lons=indgen(360/45+1)*45-180
IDL>lonnames=strtrim(-lons)
IDL>map_grid,latdel=10,lonnames=lonnames,lons=lons,color=0.30*!d.n_colors,/LABEL,/HORIZON
其中map_image做的事就是对数据作投影,这里用的是Mollweide投影(用map_set设置/MOLLWEIDE),
两个重要的返回值是Xsize和Ysize,如果不将这两个值传递给tvscl,那么在
PS图里就无法将坐标线和投影图对齐,但是对于Xwindow的输出没有影响,原因是PS图像素的大小是可以变的。在标注坐标线的经度的时候要注意,星图和地图正好是反的,所以有lonnames=strtrim(-lons)一句。
先草记这些,要画一幅漂亮的全天分布图其实有很大学问,等我会了再写。
https://blog.sciencenet.cn/blog-117333-373635.html
上一篇:
一些射电望远镜观测时间申请的信息下一篇:
小记IDL读FITS文件(数据标定)