Back (or double?) great circles drawn when used draw_grid from any position.

[dpshelio]

Description

We get the lines behind the sun when using draw_grid.

Expected behavior

That only the front lines ares shown and the backs are not drawn! (or at least @Cadair said: weeeird, bug report time)

Actual behavior

From SDO:

From other places:

Steps to Reproduce

1. Get the observer location

sdo_pos = get_horizons_coord('SDO', '2019-03-14 00:10:00')

2. Generate a map from such observer

def sun_seen_from(position, name):
    
    # Find out the size of the sun
    solar_r = (constants.radius / position.radius).to(u.arcsec, equivalencies=u.dimensionless_angles())
    th = np.linspace(0 * u.deg, 360 * u.deg)

    # Create a FD map from there
    half_image = 100
    datax, datay = np.mgrid[-half_image:half_image, -half_image:half_image]
    data = 1 - np.sqrt(datax**2 + datay**2)
    coord = SkyCoord(0*u.arcsec, 0*u.arcsec, obstime=position.obstime, observer=position, frame=frames.Helioprojective)
    size = solar_r / half_image / u.pixel
    header = sunpy.map.header_helper.make_fitswcs_header(data, coord,
                                                         reference_pixel=u.Quantity([half_image, half_image]*u.pixel),
                                                         scale=u.Quantity([size, size]),
                                                         instrument='sunlike', detector=f'Sun seen from {name}')

    return sunpy.map.Map(data, header)


sdo_map = sun_seen_from(sdo_pos, name='SDO')

3. Plot that map with the draw_grid

fig = plt.figure(figsize=(7,8))
ax = fig.add_subplot(1, 1, 1, projection=sdo_test_map)
sdo_test_map.plot(axes=ax)
sdo_test_map.draw_grid(axes=ax)
fig.savefig('Earth_sun.png')
plt.show()

4. Repeat for other positions around the heliosphere.

System Details

• SunPy Version: v1.0 🎉
• Python Version: 3.6.7
• OS information: Linux using miniconda (conda v.4.5.11)

[ayshih]

(I haven't yet tried what's in the original post.)

I haven't been able to reproduce this issue when using a Map object from an actual file (with manually modifying the observer coordinate to be sure). No back-side grid lines appear.

>>> import astropy.units as u
>>> import sunpy.data.sample
>>> from sunpy.map import Map
>>> m = Map(sunpy.data.sample.EIT_195_IMAGE)
>>> m.meta['hec_z'] *= 1000
>>> m.peek(clip_interval=[5, 99]*u.percent, draw_grid=True)

[Cadair]

So in collaboration with @astrofrog I think this is something to do with how the coordinates round trip HGS -> pixel -> HGS. WCSAxes uses how that transform behaves (i.e. does input match output) to work out of the lines should be drawn. I have so far got this down to the following behaviour:

A custom map:

Point on rear of sun
<SkyCoord (HeliographicStonyhurst: obstime=None): (lon, lat, radius) in (deg, deg, km)
    (-180., 0., 695700.)>
pixel coord
PixelPair(x=<Quantity 27.74920045 pix>, y=<Quantity 104.42559215 pix>)
WCS native coord
<SkyCoord (Helioprojective: obstime=2019-06-03T20:28:40.013, rsun=695700.0 m, observer=<HeliographicStonyhurst Coordinate (obstime=2019-06-03T20:28:40.013): (lon, lat, radius) in (deg, deg, m)
    (-133.44079521, -3.26643583, 2.44214968e+11)>): (Tx, Ty) in arcsec
    (-427.47626507, 23.06636179)>
HGS coord
<SkyCoord (HeliographicStonyhurst: obstime=2019-06-03T20:28:40.013): (lon, lat, radius) in (deg, deg, km)
    (nan, nan, nan)>

AIA_171_IMAGE map:

Point on rear of sun
<SkyCoord (HeliographicStonyhurst: obstime=None): (lon, lat, radius) in (deg, deg, km)
    (-180., 0., 695700.)>
pixel coord
PixelPair(x=<Quantity 510.13144745 pix>, y=<Quantity 511.24706341 pix>)
WCS native coord
<SkyCoord (Helioprojective: obstime=2011-06-07T06:33:02.770, rsun=696000000.0 m, observer=<HeliographicStonyhurst Coordinate (obstime=2011-06-07T06:33:02.770): (lon, lat, radius) in (deg, deg, m)
    (359.99593689, 0.04787238, 1.51846026e+11)>): (Tx, Ty) in arcsec
    (-0.06671055, 0.78599687)>
HGS coord
<SkyCoord (HeliographicStonyhurst: obstime=2011-06-07T06:33:02.770): (lon, lat, radius) in (deg, deg, km)
    (-0.00808742, 0.09528755, 695999.99999955)>

I don't understand why the custom map returns NaNs for the last HPC -> HGS transform. Investigations continue.

[Cadair]

The culprit is rsun=695700.0 m

[ayshih]

Heh, the incredible shrinking Sun!

[ayshih]

Incidentally, the Sun's Schwarzschild radius is 3 km, so shrinking it by three orders of magnitude (to a radius a tenth of Earth's radius) still didn't turn it into a black hole.

[Cadair]

I am glad you provided that information as I was wondering.

It is worth noting that this is a particularly weird WCSAxes failure mode (round trip calculations failing if the result is NaN) so it might also need an upstream bug report.