使用 hdf5 文件中的 astropy 在 skyplot 上绘制平均值和标准 dev 值

Question

我正在尝试创建一个包含来自 hdf5 文件的平均值和标准开发值的天空图（使用 astropy）。数据链接为https://wwwmpa.mpa-garching.mpg.de/~ensslin/research/data/faraday2020.html（Faraday Sky 2020）。 到目前为止，我已经编写了以下代码，其中数据从 hdf5 文件读取到 ggl 和 ggb，然后将值转换为 gb 和 gl 中的银河坐标（l 和 b）。我需要在天空图中绘制这些值。

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

dat = []

ggl=[]

ggb=[]

with h5py.File('faraday2020.hdf5','r') as hdf:
    print(list(hdf.keys()))
    faraday_sky_mean = hdf['faraday_sky_mean'][:]
    faraday_sky_std = hdf['faraday_sky_std'][:]
    
print(faraday_sky_mean.shape, faraday_sky_mean.dtype) 
print(f'Max Mean={max(faraday_sky_mean)}, Min Mean={min(faraday_sky_mean)}') 
print(faraday_sky_std.shape, faraday_sky_std.dtype) 
print(f'Max StdDev={max(faraday_sky_std)}, Min StdDev={min(faraday_sky_std)}') 

ggl = faraday_sky_mean.tolist()
print(len(ggl),type(ggl[0]))
ggb = faraday_sky_std.tolist()
print(len(ggb),type(ggb[0]))

gl = ggl * u.degree
gb = ggb * u.degree


c = SkyCoord(l=gl,b=gb, frame='galactic', unit = (u.deg, u.deg)) #, 

l_rad = c.l.wrap_at(180 * u.deg).radian
b_rad = c.b.radian

###
plt.figure(figsize=(8,4.2))
plt.subplot(111, projection="aitoff")

plt.title("Mean and standard dev", y=1.08, fontsize=20)
plt.grid(True)

P1=plt.plot(l_rad, b_rad,c="blue", s=220, marker="h", alpha=0.7) #, 

plt.subplots_adjust(top=0.95, bottom=0.0)
plt.xlabel('l (deg)', fontsize=20)
plt.ylabel('b (deg)', fontsize=20)

plt.subplots_adjust(top=0.95, bottom=0.0)
plt.show()

但是，我收到以下错误：

'got {}'.format(angles.to(u.degree)))

ValueError: Latitude angle(s) must be within -90 deg <= angle <= 90 deg, got [1.12490771 0.95323024 0.99124631 ... 4.23648627 4.28821608 5.14498169] deg

请帮助我解决这个问题。

Answer 1

这是对我之前答案的扩展。最初的帖子想要在天体星图上绘制法拉第天空 2020 数据的均值和标准偏差。引用的数据源（来自 Radboud 大学）仅包括平均值和标准差。相关的经度和纬度坐标是从 NASA Goddard LAMBDA-Tools 站点获得的。下面的代码显示了如何将两个文件中的数据合并到一个 HDF5 文件中。为方便起见，此处重复了指向数据源的链接：

• 链接到Faraday Sky 2020 data
• 链接到HEALPix Pixel Coordinates

生成的文件名为：“faraday2020_with_coords.h5”。

from astropy.io import fits
import h5py

fits_file = 'pixel_coords_map_ring_galactic_res9.fits'
faraday_file = 'faraday2020.hdf5'

with fits.open(fits_file) as hdul, \    
     h5py.File(faraday_file,'r') as h5r, \
     h5py.File('faraday2020_with_coords.h5','w') as h5w:

    arr = hdul[1].data

    dt = [('LONGITUDE', float), ('LATITUDE', float), \
          ('faraday_sky_mean', float), ('faraday_sky_std', float) ]
                             
    ds = h5w.create_dataset('skyplotdata', shape=(arr.shape[0],), dtype=dt)
    ds['LONGITUDE'] = arr['LONGITUDE'][:]
    ds['LATITUDE']  = arr['LATITUDE'][:]
    ds['faraday_sky_mean'] = h5r['faraday_sky_mean'][:]
    ds['faraday_sky_std']  = h5r['faraday_sky_std'][:]

Answer 2

我明白为什么您在绘制此数据时遇到问题。链接文件（faraday2020.hdf5）中的数据只是重建法拉第天空的平均值和标准差。请参阅链接页面上的此注释：“所有地图均以 HEALPix Nside=512 的分辨率在银河中呈现，并存储在 RING 排序方案中。单位为 rad/m2。”换句话说，您需要从另一个来源获取天空图坐标。

在 NASA Goddard LAMBDA-Tools 网站上搜索了一下坐标：HEALPix Pixel Coordinates。具体来说，您需要此文件用于 NSide=512 / Galactic / Ring Pixel Ordering：pixel_coords_map_ring_galactic_res9.fits

所以，第一个问题解决了。接下来，您需要读取 FITS 格式的文件以获取坐标。 Astropy 有“fits”模块来做到这一点。请参阅下面的代码。

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

filename='pixel_coords_map_ring_galactic_res9.fits'

with fits.open(filename) as hdul:
    print(hdul.info())
    arr = hdul[1].data
    print(arr.shape)
# Returns:
# (3145728,)
    print(arr.dtype)
# Returns:
# dtype((numpy.record, [('LONGITUDE', '>f4'), ('LATITUDE', '>f4')]))

ggl = arr['LONGITUDE'][:].tolist()
ggb = arr['LATITUDE'][:].tolist() 

gl = ggl * u.degree
gb = ggb * u.degree

c = SkyCoord(l=gl,b=gb, frame='galactic', unit = (u.deg, u.deg))  

l_rad = c.l.wrap_at(180 * u.deg).radian
b_rad = c.b.radian

上面的代码为您的天空图坐标提供了l_rad 和b_rad。接下来，您需要合并我之前给您的代码以读取 Farady Sky Mean 和 StdDev。

with h5py.File('faraday2020.hdf5','r') as hdf:
    faraday_sky_mean = hdf['faraday_sky_mean'][:]
    faraday_sky_std = hdf['faraday_sky_std'][:]

最后，用 matplotlib 绘制两组数据。我将绘图更改为使用散点图，用c=faraday_sky_mean（平均值）对标记进行颜色编码。您可以对faraday_sky_stddev 执行相同操作以获取标准差值。

plt.figure(figsize=(8,4.2))
plt.subplot(111, projection="aitoff")

plt.title("Mean", y=1.08, fontsize=20)
plt.grid(True)
# P1=plt.plot(l_rad, b_rad,c="blue", marker="h", alpha=0.7) #, s=220) 
P2 = plt.scatter(l_rad, b_rad, s=20, c=faraday_sky_mean, cmap='hsv')

plt.subplots_adjust(top=0.95, bottom=0.0)
plt.xlabel('l (deg)', fontsize=20)
plt.ylabel('b (deg)', fontsize=20)

plt.subplots_adjust(top=0.95, bottom=0.0)
plt.show()
print('DONE')

把它们放在一起，你会得到下面的图片。我认为这是准确的（但对天体物理学一无所知，所以不能 100% 确定）。这应该让你指向正确的方向。祝你好运。