Sunpy map swaps WCS Naxis

[Vishal Upendran]

Hi all,

I was trying to use some raster data from iris (at one wavelength alone), and found the map object to be potentially useful. However, I ran into an issue while trying to obtain the `bottom_left` and `top_right` coords of the image. It seems that sunpy map reverses the WCS NAXIS. Code snippet and results below for reference:

head = "XTENSION= 'IMAGE ' / IMAGE extension BITPIX = 16 / Number of bits per data pixel NAXIS = 3 / Number of data axes NAXIS1 = 210 / NAXIS2 = 1095 / NAXIS3 = 400 / PCOUNT = 0 / No Group Parameters GCOUNT = 1 / One Data Group BSCALE = 0.25 / BZERO = 7992 / CDELT1 = 0.0254399999976 / CDELT2 = 0.166350 / CDELT3 = 0.349202072411 / CRPIX1 = 1.00000 / CRPIX2 = 548.000 / CRPIX3 = 200.000 / CRVAL1 = 1390.87174787 / CRVAL2 = -167.374 / CRVAL3 = 469.431 / CTYPE1 = 'WAVE ' / CTYPE2 = 'HPLT-TAN' / CTYPE3 = 'HPLN-TAN' / CUNIT1 = 'Angstrom' / CUNIT2 = 'arcsec ' / CUNIT3 = 'arcsec ' / PC1_1 = 1.00000000000 / PC1_2 = 0.00000000000 / PC2_1 = 0.00000000000 / PC2_2 = 0.999938666821 / PC3_1 = 0.00000000000 / PC3_2 = 0.00537130899584 / PC3_3 = 0.999938666821 / PC2_3 = -0.0236694471641 /"

raster_data = np.ones([1095, 400]) #Raster data at one wavelength

wcs_rast = WCS(head)  #WCS from head

rast_map = mp.Map(raster_data,wcs_rast.dropaxis(0)) #Drop the wavelength axis.

print(rast_map.wcs)
-> WCS Keywords Number of WCS axes: 2 CTYPE : 'HPLT-TAN' 'HPLN-TAN' CRVAL : -0.046492777777778 0.1303975 CRPIX : 548.0 200.0 PC1_1 PC1_2 : 0.999938666821 -0.0236694471641 PC2_1 PC2_2 : 0.00537130899584 0.999938666821 CDELT : 4.6208333333333e-05 9.7000575669722e-05 NAXIS : 400 1095

print(wcs_rast.dropaxis(0))

-> Number of WCS axes: 2 CTYPE : 'HPLT-TAN' 'HPLN-TAN' CRVAL : -0.04649277777777778 0.1303975 CRPIX : 548.0 200.0 PC1_1 PC1_2 : 0.999938666821 -0.0236694471641 PC2_1 PC2_2 : 0.00537130899584 0.999938666821 CDELT : 4.620833333333333e-05 9.700057566972222e-05 NAXIS : 1095 400

Now, when I try to get the coordinates:

print(rast_map.bottom_left_coord)

->(Tx, Ty) in arcsec (-257.57830319, 398.91801151)

print(rast_map.wcs.pixel_to_world(0*u.pix, 0*u.pix))

-> (Tx, Ty) in arcsec (398.91801151, -257.57830319)

print(rast_map.wcs.pixel_to_world(1095*u.pix, 400*u.pix))

->(Tx, Ty) in arcsec (540.64417529, -77.01121223)

print(rast_map.pixel_to_world(*(u.Quantity(rast_map.dimensions)-1*u.pix))) #Sunpy code

-> (Tx, Ty) in arcsec (-195.51605002, 781.66998537)

print(rast_map.top_right_coord)

-> (Tx, Ty) in arcsec (-195.51605002, 781.66998537)

Am I doing something wrong, or is this an issue in map? I have also attached the above code if someone would want to run it.

Cheers,

Vishal

[David Stansby]

Hi Vishal,

Thanks for the detailed report! I've double checked, and I think you have found a bug in sunpy. Essentially at one point it is assumed that CTYPE1 is longitude, and CTYPE2 is latitude. I've opened a fix here: https://github.com/sunpy/sunpy/pull/4700, which we hope to get backported to the 2.x release. This might take a little while, but I'll let you know when it's fixed. Thanks again for the report.

All the best,

David

--
You received this message because you are subscribed to the Google Groups "SunPy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sunpy+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/sunpy/b41aa16c-1fbd-45fc-93c9-81d5cc38a8e4n%40googlegroups.com.

[Vishal Upendran]

Hi David,

Thanks for the quick reply! I would imagine if CTYPE1 is assumed to be longitude and such, wouldn't the corresponding CRVALs, CDELTs, etc change accordingly? Why is NAXIS alone changed?

And thanks for the quick PR. I shall keep following the conversation there.

Thanks,

Vishal

To view this discussion on the web visit https://groups.google.com/d/msgid/sunpy/CAGm9cHB212wKa7s9DEJdaPV5LMpbddR2WhDnOWS-aBU%3DrhsXwg%40mail.gmail.com.

--

Vishal Upendran,

Senior Research Fellow, Inter-University Centre for Astronomy and Astrophysics (IUCAA),

Pune, India - 411007

[David Stansby]

Hi Vishal,

Sorry for the delay; I think running .dropaxis(0) on the WCS object is an operation that takes part purely in astropy, so I think the best thing there would be to either email their mailing list (https://mail.python.org/mailman/listinfo/astropy) or open an issue on their bug tracker. We're hopefully just about to release a bugfix release that will fix the coordinates being swapped in sunpy.

Cheers,

David

To view this discussion on the web visit https://groups.google.com/d/msgid/sunpy/CAEjf1m14n3PVWPqaqNfoHDUanwauQOba992aq_aJW%3DO6adU-iw%40mail.gmail.com.

[Vishal Upendran]

Hi David,

Thanks for your reply. I shall investigate, and see if an issue needs to be submitted to astropy!

Cheers,

Vishal

You received this message because you are subscribed to a topic in the Google Groups "SunPy" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/sunpy/oiwGCLoZENI/unsubscribe.
To unsubscribe from this group and all its topics, send an email to sunpy+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/sunpy/CAGm9cHCxEBQT2dLWkd1PDNvRiKFa-Pf%3DvkXSWq-_MLXzY2iJrg%40mail.gmail.com.