SkyPlot heatmaps aren't rendered correctly when values exist with RA close to 0 or 360 degrees

[derekocallaghan]

I've noticed this issue when visualizing localizations with EquatorialPlot, where the probability heatmap isn't rendered correctly when the centroid RA is close to 0 or 360 degrees. This can be illustrated using a local copy of glg_healpix_all_bn220130968_v01.fit:

from gdt.core.plot.sky import EquatorialPlot
from gdt.missions.fermi.gbm.localization import GbmHealPix

gbm_healpix = GbmHealPix.open("./glg_healpix_all_bn220130968_v01.fit")
gbm_healpix.centroid
(357.87735849056605, -50.480044262442945)

eqplot = EquatorialPlot()
eqplot.add_localization(gbm_healpix, gradient=True)

The following plot is generated:

According to the GBM catalog, the correct skymap should look like:

The issue appears to be caused by the transformation of specified lon/lat arrays to SkyCoord in SkyMap.plot_heatmap():

gdt-core/src/gdt/core/plot/sky.py

Lines 456 to 467 in 91020a3

	coords = SkyCoord(x.flatten(), y.flatten(), frame='gcrs', unit='deg')
	try:
	coords = coords.transform_to(self._astropy_frame)
	except:
	pass
	lon, lat = get_lonlat(coords)
	lon = lon.reshape(x.shape)
	lat = lat.reshape(y.shape)
	heatmap = SkyHeatmap(lon.value, lat.value, heatmap, self.ax,
	frame=self._frame, flipped=self._flipped, **kwargs)

In the case of localizations, the lon_array specified to plot_heatmap() will contain values in $[0, 360]$. However, the call to SkyCoord() will convert the last 360 value to 0. E.g.:

from astropy.coordinates import SkyCoord
SkyCoord(ra=[0, 180, 360], dec=[45, 45, 45], unit="deg")

<SkyCoord (ICRS): (ra, dec) in deg
    [(  0., 45.), (180., 45.), (  0., 45.)]>

This appears to cause the issue when subsequently rendering the heatmap. Note that it's only visible when centroids straddle RA=0/360 degrees, other localisation centroids are rendered correctly. I've also noticed a similar issue with the contour lineswhere although it doesn't impact the visualisation as much as the heatmaps, the lines are broken, and I assume it's also due to the use of SkyCoord.

I have a local workaround that overrides plot_heatmap() to exclude the call to SkyCoord(), where the specified lon_array and lat_array are passed directly to SkyHeatmap, similar to the previous GBM Data Tools implementation:

from gdt.core.plot.plot import SkyHeatmap

class FixEquatorialPlot(EquatorialPlot):
    def plot_heatmap(self, heatmap, lon_array, lat_array, **kwargs):
        """Plot a heatmap on the sky.
           This is a workaround for the issue in the parent SkyPlot.plot_heatmap(),
           where the meshgrid results in a lon_array starting/ending with 0 degrees.
           This causes issues where the heatmap generated from a HEALPix location
           is populated at both RA=0/360 degrees, resulting in an unexpected output
           from the subsequent call to Matplotlib's pcolormesh().

           The old GBM Data Tools didn't generate a meshgrid, and just passes the
           specified coordinate arrays unchanged.

        Args:
            heatmap (np.array): A 2D array of values
            ra_array (np.array): The array of longitudinal gridpoints
            dec_array (np.array): The array of latitudinal gridpoints
            **kwargs: Options to pass to :class:`~gdt.plot.plot.SkyHeatmap`
        """
        heatmap = SkyHeatmap(
            lon_array,
            lat_array,
            heatmap,
            self.ax,
            frame=self._frame,
            flipped=self._flipped,
            **kwargs
        )

        return heatmap

This results in a correct rendering of the heatmap, similar to that contained in the catalo:

fixeqplot = FixEquatorialPlot()
fixeqplot.add_localization(gbm_healpix, gradient=True)

I can create a PR with this fix but I wanted to check first whether the SkyCoord usage was required for other reasons.

[derekocallaghan]

Hi @AdamGoldstein-USRA,

Just wanted to check whether the suggested workaround that excludes the call to SkyCoord() is okay with you? If so, I'll create a PR.

Thanks

[BillCleveland-USRA]

Met with Adam. Agree this needs further research to make sure the workaround won't break any other functionality in GDT.

[AdamGoldstein-USRA]

@derekocallaghan Please checkout PR #38 and see if it works for you.

Your proposed fix would work only for EquatorialPlot, and the problem would still be present for GalacticPlot and SpacecraftPlot. The intent is to allow plotting in other frames (currently Galactic and spacecraft), and after you identified the root cause of the issue, I think the fix in PR #38 addresses that for any frame.

[derekocallaghan]

Thanks @AdamGoldstein-USRA, that was just a quick workaround to get the correct plot here. I've tried the modifed lib.py and sky.py and the plot is rendered as expected with the following, so I'll close this issue:

eqplot = EquatorialPlot()
eqplot.add_localization(gbm_healpix, gradient=True)

Fixed by #38