Sky region should wrap at map edges

[jlenain]

Hi all,

First of all, many thanks for the good work!

I am experiencing an issue when trying to plot a CircleSkyRegion at a map edge in Aitoff projection. I followed the issue #76 to properly project a circle on the sky in this projection.

My test displays a polygon which is filled across the whole longitude band, but empty where it is expected to be filled. I attach here a plot illustrating this, as well as a minimal working example:

#!/bin/env python
import matplotlib.pyplot as plt
import astropy.coordinates as c
import astropy.units as u
from gammapy.maps import WcsMapND
from regions import CircleSkyRegion

fig = plt.figure()
m = WcsMapND.create(binsz=1.0, width=(360., 180.), skydir=(0.0, 0.0),
                    coordsys='GAL', proj='AIT')
m.fill_poisson(1.0)
wcs = m.geom.wcs
ax = fig.add_subplot(111, projection=wcs)
ax.imshow(m.data, cmap='plasma', alpha=0.2)

center = c.SkyCoord(l=180.2*u.degree, b=61.*u.degree, frame='galactic')
rad = 20.0*u.degree
region = CircleSkyRegion(center=center, radius=c.Angle(rad))
region.to_pixel(wcs=wcs, mode='full').plot(ax, color='green')

plt.show()

Is there a known workaround ?

Many thanks in advance.

[jlenain]

Hi again, As a bonus, I forgot to mention that the region should be displayed in the Northern hemisphere, not the Southern one, if I did not miss something.

[astrofrog]

@jlenain - could you try with the latest developer version of regions just to make sure this is still relevant for now? I think we decided to temporarily remove the 'spherical circle' region, though when we add it back we should take into account this issue (which is not trivial)

[jlenain]

@astrofrog - Many thanks for your reply. I indeed tested this using v0.2 of regions in the first place. Testing with the latest one (0.3.dev496), I can say that the mode argument of the CircleSkyRegion.to_pixel method was removed, so, no appropriate sky projection any more. With the following modification:

region.to_pixel(wcs=wcs).plot(ax, color='green') the result is: which is even more puzzling to me. I remind you that the tested position for the center of the region is l=180.2°, b=61.0°. Many thanks in any case!

[astrofrog]

About the vertical flip - you are missing origin='lower' in your call to imshow.

I'll try and work soon on adding back the spherical circle class.

[jlenain]

@astrofrog - Oh my, thanks a lot for pointing this out. And thanks again for your efforts.

[cdeil]

@jlenain - So what you want is that the sky circle shows up correct, with two separated patches at different edges of the image, right?

I don't think there is a simple solution at the moment. I might be wrong, but I think no package currently has the code built-in to take your circle parameters and to generate more than one matplotlib patch objects (probably polygon, but maybe MPL also has bezier curves or something else). @astrofrog @adonath - Is that true, or can MPL do this?

So at the moment to implement this, I think you have to first polygonise the circle (for this you can use the code from Tom from the SphericalCircle class you were using earlier), then convert the polygon to pixels, and then write some custom lines to split the polygon in two parts, then plot it with matplotlib. If you really need this plot, and have trouble getting it to work, let me know and I'll give it a shot later today.

We needed something similar for the HESS survey paper, to plot the "not visible to HESS" region on the Figure 1 plot. There we opted for making a mask and using a MPL contour call to plot the region: https://gist.github.com/cdeil/9497c96a369cfc878d24c1405de1c310 This might also be a workable solution for your application.

[jlenain]

Hi @cdeil, many thanks for your input. I think for the time being, I will work around with your proposed solution, as per the HESS survey paper. Many thanks for pointing out your script.