HI4PI：德国埃菲尔斯伯格的100米望远镜和澳大利亚Parkes的64米望远镜，共同绘制了北半球和南半球中性氢的详细结构，它揭示了银河系及周边气体分布大尺度结构的大量细节。

数据下载链接：https://lambda.gsfc.nasa.gov/product/foreground/fg_hi4pi_get.html

业余射电天文交流群：926020514

PlotHI4PI

绘制HI4PI全天区的巡天图

from astropy.io import fits
import numpy as np
import matplotlib.pyplot as plt
from astropy.coordinates import SkyCoord
from astropy import units as u

hdul = fits.open(r'/home/walnut/FUDS/HI4PI/NHI_HPX.fits') #这里改成NHI_HPX.fits文件放置的路径

l = hdul[1].data['GLON']
b = hdul[1].data['GLAT']
hi_density = hdul[1].data['NHI']

coords = SkyCoord(l=l * u.degree, b=b * u.degree,frame = 'galactic')

'''
center_l = (30) * u.degree  
center_b = (0) * u.degree  
delta_l = 5.0* u.degree
delta_b = 5.0* u.degree

cent_coords = SkyCoord(l=center_l, b=center_b ,frame = 'galactic')
indices = np.where(
    (coords.l >= (center_l - delta_l)) & (coords.l <= (center_l + delta_l)) &
    (coords.b >= (center_b - delta_b)) & (coords.b <= (center_b + delta_b))
)
'''
fig = plt.figure(figsize=(18, 10))
ax = fig.add_subplot(111, projection='mollweide')

sc = ax.scatter(coords.l.wrap_at(180 * u.deg).radian, coords.b.radian, c = hi_density, cmap='viridis', s = 5)
cbar = plt.colorbar(sc, ax=ax, label='NHI')
ax.set_xlabel('l (degrees)')
ax.set_ylabel('b (degrees)')

plt.grid(True)
plt.show()

PlotHI4PI-银道坐标指定天区

import numpy as np
import matplotlib.pyplot as plt
from astropy.coordinates import SkyCoord
from astropy import units as u

l = hdul[1].data['GLON']
b = hdul[1].data['GLAT']
hi_density = hdul[1].data['NHI']

coords = SkyCoord(l=l * u.degree, b=b * u.degree,frame = 'galactic')

center_l = 90 * u.degree  #中心点的银经
center_b = 0 * u.degree   #中心点的银纬

delta_l = 10.0* u.degree  #显示范围
delta_b = 10.0* u.degree  #显示范围

indices = np.where(
    (coords.l >= (center_l - delta_l)) & (coords.l <= (center_l + delta_l)) &
    (coords.b >= (center_b - delta_b)) & (coords.b <= (center_b + delta_b))
)

plt.figure(figsize=(18,10))
plt.gca().invert_xaxis()
plt.scatter(coords.l[indices], coords.b[indices], c=hi_density[indices], cmap='viridis', s=50)
plt.colorbar(label='NHI')
plt.xlabel('l[deg]')
plt.ylabel('b[deg]')
plt.grid(True)
plt.show()

PlotHI4PI-赤道坐标指定天区

import numpy as np
import matplotlib.pyplot as plt
from astropy.coordinates import SkyCoord
from astropy import units as u

ra = hdul[1].data['RA2000']
dec = hdul[1].data['DEC2000']
hi_density = hdul[1].data['NHI']
coords = SkyCoord(ra=ra * u.degree, dec=dec * u.degree)

star_data = {
    'M31': {'RA': 10.6848, 'Dec': 41.26916666666666},
    'M33': {'RA': 23.4620690621800, 'Dec': 30.6601751119800},
    'LMC': {'RA': 80.89375, 'Dec': -69.75611111111111},
    'SMC': {'RA': 13.186666666666667, 'Dec': -72.8286111111111},
}

star_name_to_find = 'M31' #从这里改想查看的目标天体的名字

if star_name_to_find in star_data:
    center_ra = star_data[star_name_to_find]['RA']* u.degree
    center_dec = star_data[star_name_to_find]['Dec']* u.degree

else:
    print(f"找不到天体体 {star_name_to_find} 的信息。")


#center_ra = 15.5 * u.degree  
#center_dec = 35 * u.degree  

delta_ra = 5.0* u.degree  #显示范围
delta_dec = 5.0* u.degree

indices = np.where(
    (coords.ra >= (center_ra - delta_ra)) & (coords.ra <= (center_ra + delta_ra)) &
    (coords.dec >= (center_dec - delta_dec)) & (coords.dec <= (center_dec + delta_dec))
)

plt.figure(figsize=(18,10))
plt.gca().invert_xaxis()
plt.scatter(coords.ra[indices], coords.dec[indices], c=hi_density[indices], cmap='viridis', s=50)
plt.colorbar(label='NHI')
plt.xlabel('RA')
plt.ylabel('Dec')
plt.grid(True)
plt.show()

使用天球坐标系显示指定天区：

import matplotlib.ticker as mticker
import cartopy.mpl.ticker as cticker 
import cartopy.crs as ccrs

def convert_to_minus_180_to_180(degrees):
    degrees = np.mod(degrees, 360)
    degrees[degrees > 180] -= 360
    return degrees

fig, ax = plt.subplots(figsize=(18, 10),subplot_kw={'projection': ccrs.Mollweide()})

ra = hdul[1].data['RA2000']
dec = hdul[1].data['DEC2000']
hi_density = hdul[1].data['NHI']
coords = SkyCoord(ra=ra * u.degree, dec=dec * u.degree)

star_data = {
    'M31': {'RA': 10.6848, 'Dec': 41.26916666666666},
    'M33': {'RA': 23.4620690621800, 'Dec': 30.6601751119800},
    'LMC': {'RA': 80.89375, 'Dec': -69.75611111111111},
    'SMC': {'RA': 13.186666666666667, 'Dec': -72.8286111111111},
}

star_name_to_find = 'SMC'  #改目标天体的名字

if star_name_to_find in star_data:
    center_ra = star_data[star_name_to_find]['RA']* u.degree
    center_dec = star_data[star_name_to_find]['Dec']* u.degree

else:
    print(f"找不到天体体 {star_name_to_find} 的信息。")

#center_ra = 0 * u.degree  
#center_dec = 0 * u.degree  
delta_ra = 10.0* u.degree  #显示范围
delta_dec = 5.0* u.degree

indices = np.where(
    (coords.ra >= (center_ra - delta_ra)) & (coords.ra <= (center_ra + delta_ra)) &
    (coords.dec >= (center_dec - delta_dec)) & (coords.dec <= (center_dec + delta_dec))
)

projected_data = ax.projection.transform_points(ccrs.PlateCarree(), coords.ra[indices], coords.dec[indices])
x = projected_data[:, 0]
y = projected_data[:, 1]

sc = ax.scatter(x, y, c=hi_density[indices], cmap='viridis', s=50)
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False, linestyle='--', linewidth=0.5)

#ax.set_extent([-10, 10, -10, 10])
cbar = plt.colorbar(sc, ax=ax, label='NHI')
plt.gca().invert_xaxis()
plt.grid(True)
plt.show()

天体升起下落图

import matplotlib.pyplot as plt
import numpy as np
from astropy.visualization import astropy_mpl_style, quantity_support
from astropy.coordinates import AltAz, EarthLocation, SkyCoord
from astropy.time import Time
from matplotlib.ticker import MultipleLocator
from astropy import units as u

star_data = {
    'M31': {'RA': 10.6848, 'Dec': 41.26916666666666, 'frame': 'icrs'}, #赤经，赤纬
    'M33': {'RA': 23.4620690621800, 'Dec': 30.6601751119800, 'frame': 'icrs'},
    'LMC': {'RA': 80.89375, 'Dec': -69.75611111111111, 'frame': 'icrs'},
    'SMC': {'RA': 13.186666666666667, 'Dec': -72.8286111111111, 'frame': 'icrs'},
    'G000': {'l': 0, 'b': 0, 'frame': 'galactic'},                    #银经，银维
    'G030': {'l': 30, 'b': 0, 'frame': 'galactic'},                    
    'G090': {'l': 90, 'b': 0, 'frame': 'galactic'},  
}

star_name_to_find = 'M31'  #从这里改目标天体的名字

if star_name_to_find in star_data:
    if star_data[star_name_to_find]['frame'] == 'icrs':
        center_ra = star_data[star_name_to_find]['RA']* u.degree
        center_dec = star_data[star_name_to_find]['Dec']* u.degree
        star = SkyCoord(ra=center_ra, dec=center_dec, frame="icrs") 
    else:
        center_l = star_data[star_name_to_find]['l']* u.degree
        center_b = star_data[star_name_to_find]['b']* u.degree
        star = SkyCoord(l=center_l, b=center_b, frame="galactic")
else:
    print(f"找不到天体体 {star_name_to_find} 的信息。")


time = Time('2023-12-30 00:00:00') 
utcoffset = 8 * u.hour
time = time - utcoffset # 转换为BJT

#观测时望远镜所在的地理位置，以北京为例
lat = (39.+54./60.0+27./3600.0)*u.degree      # 当地纬度
lon = (116.0+23.0/60.0+17./3600.0)*u.degree   # 当地经度
height =  50.0 * u.m                          # 当地海拔
location = EarthLocation(lat=lat, lon=lon, height=height)

delta_time = np.linspace(0, 24, 288) * u.hour
times = time + delta_time

frame = AltAz(obstime=times, location=location)
staraltazs = star.transform_to(frame)

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(18, 25), sharex=True)
# Plot for Altitude
ax1.plot(delta_time, staraltazs.alt, label=star_name_to_find, c='k')
ax1.set_ylabel('Altitude [deg]')
ax1.set_title((time + 8*u.hour).value[0:10])
ax1.legend(fontsize=15)
ax1.axhline(30, c='r')
ax1.set_xticks(np.arange(0, 25, 1))
ax1.set_yticks(np.arange(-95, 95, 5))
ax1.set_ylim(staraltazs.alt.min()/u.deg*1.1, staraltazs.alt.max()/u.deg*1.1)

# Plot for Azimuth
def convert_to_minus_180_to_180(degrees):
    degrees = np.mod(degrees, 360)
    degrees[degrees > 180] -= 360
    return degrees

az = convert_to_minus_180_to_180(staraltazs.az/u.deg)*u.deg
ax2.plot(delta_time, az, label=star_name_to_find, c='k')
ax2.set_xlabel('Time (CST) [hour]')
ax2.set_ylabel('Azimuth [deg]')
ax2.legend(fontsize=15)
ax2.set_yticks(np.arange(-180, 180, 10))
ax2.set_xticks(np.arange(0, 25, 1))