使用 hdf5 文件中的 astropy 在 skyplot 上绘制平均和标准 dev 值
Plotting mean and standard dev values on skyplot using astropy from hdf5 file
我正在尝试创建一个 skyplot(使用 astropy),其中包含来自 hdf5 文件的平均和标准 dev 值。数据的 link 是 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
请帮我解决这个问题。
我明白你为什么在绘制这些数据时遇到问题。链接文件 (faraday2020.hdf5
) 中的数据只是重建的法拉第天空的平均值和标准差。请参阅链接页面上的注释:“所有地图均以 Nside=512 的 HEALPix 分辨率呈现在 Galactic 中,并存储在 RING 排序方案中。单位为 rad/m2。" 换句话说,您需要从其他来源获取天空图坐标。
稍微谷歌一下就在 NASA Goddard LAMBDA-Tools 网站上找到了坐标:HEALPix Pixel Coordinates. Specifically, you want this file for 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% 确定)。这应该让你指向正确的方向。祝你好运。
这是对我之前回答的扩展。最初的 post 想要在天体星空图上绘制法拉第天空 2020 数据的均值和标准差。参考数据源(来自 Radboud 大学)仅包括均值和标准差。相关的经度和纬度坐标是从 NASA Goddard LAMBDA-Tools 站点获得的。下面的代码显示了如何将两个文件中的数据合并到一个 HDF5 文件中。为方便起见,此处重复数据源链接:
- Link 到 Faraday Sky 2020 data
- Link 到 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'][:]
我正在尝试创建一个 skyplot(使用 astropy),其中包含来自 hdf5 文件的平均和标准 dev 值。数据的 link 是 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
请帮我解决这个问题。
我明白你为什么在绘制这些数据时遇到问题。链接文件 (faraday2020.hdf5
) 中的数据只是重建的法拉第天空的平均值和标准差。请参阅链接页面上的注释:“所有地图均以 Nside=512 的 HEALPix 分辨率呈现在 Galactic 中,并存储在 RING 排序方案中。单位为 rad/m2。" 换句话说,您需要从其他来源获取天空图坐标。
稍微谷歌一下就在 NASA Goddard LAMBDA-Tools 网站上找到了坐标:HEALPix Pixel Coordinates. Specifically, you want this file for 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% 确定)。这应该让你指向正确的方向。祝你好运。
这是对我之前回答的扩展。最初的 post 想要在天体星空图上绘制法拉第天空 2020 数据的均值和标准差。参考数据源(来自 Radboud 大学)仅包括均值和标准差。相关的经度和纬度坐标是从 NASA Goddard LAMBDA-Tools 站点获得的。下面的代码显示了如何将两个文件中的数据合并到一个 HDF5 文件中。为方便起见,此处重复数据源链接:
- Link 到 Faraday Sky 2020 data
- Link 到 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'][:]