import sunpy.coordinates causes bug in generating synthetic raster images

[yjzhu-solar]

Though sunpy is not a required package for euispice_coreg, I personally often use some handy functions in sunpy.coordinates, including the propagate_with_solar_surface() context manager. However, by importing sunpy.coordinates, it will change the output behaviors of WCS.pixel_to_world() if it is in helioprojective coordinates. For example:

with fits.open('../../src/SPICE/20221024/solo_L2_spice-n-ras_20221024T231535_V07_150995398-000.fits') as hdul:
    wcs_spice = WCS(hdul[0].header).dropaxis(2)
    x, y, t = np.meshgrid(np.arange(0,2),
                        np.arange(0, 2),
                        np.arange(0, 2))
    print(wcs_spice.pixel_to_world(x,y,t))

gives the following without sunpy.coordinates

[<Quantity [[[0.19955875, 0.19955875],
            [0.20066612, 0.20066612]],

           [[0.19953378, 0.19953378],
            [0.20064116, 0.20064116]]] deg>, <Quantity [[[-0.01288007, -0.01288007],
            [-0.01278912, -0.01278912]],

           [[-0.01257609, -0.01257609],
            [-0.01248514, -0.01248514]]] deg>, <Time object: scale='utc' format='mjd' value=[[[59877.10270249 59877.10271406]
  [59877.10200515 59877.10201672]]

 [[59877.10270249 59877.10271406]
  [59877.10200515 59877.10201672]]]>]

Note that the output is a list with a length of 3 (x,y,t). On the other hand, with sunpy.coordinates the output will be organized into SkyCoord objects, so the length of the output list becomes 2 (SkyCoord, t).

[<SkyCoord (Helioprojective: obstime=2022-10-25T00:51:59.753, rsun=695700.0 km, observer=<HeliographicStonyhurst Coordinate (obstime=2022-10-25T00:51:59.753, rsun=695700.0 km): (lon, lat, radius) in (deg, deg, m)
    (-50.59290147, 6.48623548, 5.92481369e+10)>): (Tx, Ty) in arcsec
    [[[(718.41148829, -46.3682358 ), (718.41148829, -46.3682358 )],
      [(722.39804299, -46.04082854), (722.39804299, -46.04082854)]],

     [[(718.32161058, -45.27392775), (718.32161058, -45.27392775)],
      [(722.30816532, -44.94652046), (722.30816532, -44.94652046)]]]>, <Time object: scale='utc' format='mjd' value=[[[59877.10270249 59877.10271406]
  [59877.10200515 59877.10201672]]

 [[59877.10270249 59877.10271406]
  [59877.10200515 59877.10201672]]]>]

This difference will cause bugs in indexing the helioprojective longitudes and latitudes in SPICEComposedMapBuilder._prepare_spectro_data() in line 234 of map_builder.py:

 longitude_spice, latitude_spice, utc_spice = w_xyt.pixel_to_world(x, y, t)

It would be helpful to notify the users either do not import sunpy.coordinates in the python file or jupyter notebook using the SPICEComposedMapBuilder or change the code accordingly.

[adolliou]

Hi Yingjie,

Thank you for pointing this out ! Yes you are correct. This a known issue with regards to sunpy and astropy. You can check the following issue for more details : sunpy/sunpy#6790

For now I added a quick fix in map_builder. I tested it and it seems to work even with sunpy.map imported. Please let me know if there is an issue. It should also be fixed for the alignment code.

Bests,
Antoine