import healpy as hp
import matplotlib.pyplot as plt
from ligo.skymap.io import read_sky_map
import ligo.skymap.plot  # registers the "astro ..." projections

url = ("https://gracedb.ligo.org/api/superevents/S190728q/files/"
       "GW190728_064510_PublicationSamples_flattened.fits.gz,0")

prob, meta = read_sky_map(url)          # prob[i]: probability of pixel i
nside = hp.npix2nside(len(prob))

ax = plt.axes(projection="astro hours mollweide")
ax.imshow_hpx((prob, "ICRS"), cmap="cylon")
ax.grid()