Matplotlib Plotting mean and standard dev values on skyplot using astropy from hdf5 file

Matplotlib Plotting mean and standard dev values on skyplot using astropy from hdf5 file,matplotlib,hdf5,astropy,polar-coordinates,fits,Matplotlib,Hdf5,Astropy,Polar Coordinates,Fits,I am trying to create a skyplot(using astropy) containing mean and standard dev values from a hdf5 file. The link to the data is https://wwwmpa.mpa-garching.mpg.de/~ensslin/research/data/faraday2020.html (Faraday Sky 2020). I have programmed the following code until now, where data is read from the hdf5 file to ggl and ggb, after which values are converted to galactic coordinates (l and b) in gb and gl. I need these values to be plotted in the skyplot. from astropy import units as u from astrop

I am trying to create a skyplot(using astropy) containing mean and standard dev values from a hdf5 file. The link to the data is https://wwwmpa.mpa-garching.mpg.de/~ensslin/research/data/faraday2020.html (Faraday Sky 2020). I have programmed the following code until now, where data is read from the hdf5 file to ggl and ggb, after which values are converted to galactic coordinates (l and b) in gb and gl. I need these values to be plotted in the skyplot.

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()

However, I am getting the following error:

'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

Please help me on how to fix this.


#1

This is an extension to my previous answer. The original post wanted to plot Mean and Standard Deviation of the Faraday Sky 2020 data on an astropy skyplot. The referenced data source (from Radboud University) only included the mean and standard deviation. The associated longitude and latitude coordinates were obtained from the NASA Goddard LAMBDA-Tools site. The code below shows how to merge the data from both files into a single HDF5 file. For convenience, links to the data sources are repeated here:

The resulting file is named: "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'][:]

#2

I see why you are having problems plotting this data. The data in the linked file (faraday2020.hdf5) is only the mean and standard deviation of the reconstructed Faraday sky. See this note on the linked page: "All maps are presented in Galactic at a HEALPix resolution of Nside=512 and are stored in RING ordering scheme. The units are rad/m2." In other words, you need to get the skyplot coordinates from another source.

A little Googling found the coordinates on the NASA Goddard LAMBDA-Tools site: HEALPix Pixel Coordinates. Specifically, you want this file for NSide=512 / Galactic / Ring Pixel Ordering: pixel_coords_map_ring_galactic_res9.fits

So, first problem solved. Next you need to read the FITS formatted file to get the coordinates. Astropy has the 'fits' module to do that. See code below.

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

The code above gives you l_rad and b_rad for your skyplot coordinates. Next, you need to merge in the code I gave you earlier to read the Farady Sky Mean and StdDev.

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

Finally, plot both sets of data with matplotlib. I changed the plot to use a scatterplot, color coding the markers with c=faraday_sky_mean (the mean values). You can do the same with faraday_sky_stddev to get the standard deviation values.

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')

Put it all together, and you will get the images below. I think this is accurate (but know nothing about astrophysics, so not 100% sure). This should get you pointed in the right direction. Good luck.

Mean Standard Deviation


#3

Check the printed output for the Max/Min values. You should have gotten: Max Mean=2857.1, Min Mean=-1406.1 and Max StdDev=696.6, Min StdDev=0.0206. Note that gl = ggl * u.degree doesn't change the magnitude, just the python class (from float for ggl[0] to astropy.units.quantity.Quantity for gl[0]).

#4

Thank you so much for your help, really appreciate it!

#5

Glad that helped. I am going to vote to close your 2 previous questions (since this one covers everything in 1 question and answer).