Plotting astronomical scatter plot - mean vs longitude (galactic)

Question

I want to create a scatter plot with the X-Axis as the longitude coordinates in the healpix file https://wwwmpa.mpa-garching.mpg.de/~ensslin/research/data/faraday2020.html (Healpix)

and the Y-Axis as the mean values in hdf5 file https://wwwmpa.mpa-garching.mpg.de/~ensslin/research/data/faraday2020.html (Faraday sky 2020)

Code until now:

from astropy.io import fits                        #libraries
from astropy import units as u
from astropy.coordinates import Galactic
import matplotlib.pyplot as plt
import h5py
from astropy_healpix import HEALPix
import numpy as np

filename='pixel_coords_map_ring_galactic_res9.fits'                       #healpix

hdulist=fits.open(filename) 
nside = hdulist[1].header['NSIDE']
order = hdulist[1].header['ORDERING']
hp = HEALPix(nside=nside, order=order, frame=Galactic())    

print(hdulist[1].header)
print(nside)
print(order)

ggl = hdulist[1].data['LONGITUDE']           #storing coordinate values in ggl and ggb
ggb = hdulist[1].data['LATITUDE'] 
print(ggl)

gl = ggl * u.degree                            #convering to galactic coordinates
gb = ggb * u.degree
print(gl)

c = Galactic(l=gl,b=gb) 
l_rad = c.l.wrap_at(180 * u.deg).radian            #X Axis
b_rad = c.b.radian

with h5py.File('faraday2020.hdf5','r') as hdf:            #importing raw data from hdf5 file
    print(hdf.keys())
    faraday_sky_mean = hdf['faraday_sky_mean'][:]           #Y Axis
    faraday_sky_std = hdf['faraday_sky_std'][:]

I have absolutely no idea how to plot a 2D square Scatter plot, given Longitude and mean are in different formats. Also, I needlongitude to be in galactic coordinates. Please help.

Answer 1

I think you're close. IMHO, this scatter plot is easier than plotting with both skyplot coordinates (projection="aitoff"). The process is similar to the answers I posted on your earlier question: Plot mean and standard dev values on skyplot using astropy from hdf5 file. You just need some minor teaks to the function parameters.

I modified your code to create a 2D scatter plot. Here's a quick summary of the differences:

• Changed from astropy.coordinates import SkyCoord (instead of HEALPix)
• Changed matplot type (removeprojection=)
• Changes y-variable from b_rad to faraday_sky_mean on scatter plot.
• Deleted c=faraday_sky_mean from plt.scatter() so data points are not color coded.

See code below.

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

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

with fits.open(fits_file) as hdulist:
    nside = hdulist[1].header['NSIDE']
    order = hdulist[1].header['ORDERING']
    #hp = HEALPix(nside=nside, order=order, frame=Galactic())    
    
    #print(hdulist[1].header)
    #print(nside)
    #print(order)
    
    ggl = hdulist[1].data['LONGITUDE']           #storing coordinate values in ggl and ggb
    ggb = hdulist[1].data['LATITUDE'] 
    #print(ggl)
    
    gl = ggl * u.degree                            #convering to galactic coordinates
    gb = ggb * u.degree
    #print(gl)
    
    #c = Galactic(l=gl,b=gb) 
    c = SkyCoord(l=gl,b=gb, frame='galactic', unit = (u.deg, u.deg))  
    l_rad = c.l.wrap_at(180 * u.deg).radian            #X Axis
    b_rad = c.b.radian
    print(len(l_rad))

with h5py.File(faraday_file,'r') as hdf:            #importing raw data from hdf5 file
    #print(hdf.keys())
    faraday_sky_mean = hdf['faraday_sky_mean'][:]           #Y Axis
    print(len(faraday_sky_mean))
    faraday_sky_std = hdf['faraday_sky_std'][:]
    
plt.figure(figsize=(8,4.2))
plt.subplot(111)

plt.title("Mean", y=1.08, fontsize=20)
plt.grid(True)
P2 = plt.scatter(l_rad, faraday_sky_mean, s=20, cmap='hsv')

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

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