Issues with coordinate system using UCoMP images with SunPy

[mgalloy]

Issue reported in discussion at the UCoMP Workshop:

First, If I directly use the header of the header of the primary extension and any data extension to create a Sunpy Map. I got an error because it seems Sunpy or the Astropy WCS module does not recognize the "CRUNIT" key. So I copied the "CRUNIT1" to "CUNIT1" and solved the problem.

Second, when I tried to reproject an AIA map to the UCoMP WCS created from the FITS header, I found it missed a few keywords like CTYPE, NAXIS1 and NAXIS2. The FITS header does not provide the observer information, for example, obsgeo_x, obsgeo_y, and obsgeo_z of MLSO either. Therefore I cannot do this projection. Fortunately, it seems that when I load UCoMP data into a Sunpy Map, Sunpy will automatically correct or refill most of these missing keywords.

Example Jupyter notebook attached.

ucomp_sunpy.pdf

Tasks

• [x] change CRUNIT{1,2} to CUNIT{1,2}
• [x] add CTYPE{1,2}
• [x] add DSUN_OBS, HGLN_OBS, and HGLT_OBS to specify Earth observing location
• [x] add comment for CTYPE{1,2}
• [x] add PCi_j keywords like KCor
• [x] get UCoMP file to successfully display with AIA image in above Jupyter notebook
• [x] we should add INHERIT FITS keyword in extensions (set to T), so that when users read an extension, they get the primary header values in the extension header (needed for display particularly)
• [x] set MJD-OBS and MJD-END corresponding to DATE-OBS and DATE-END (Modified Julian data which is JD - 2,400,000.5)

Questions

• [x] how to specify Earth-based observer with DSUN_OBS, HGLN_OBS, and HGLT_OBS
• [x] how to get aia.lev1.193A_2022-04-07T18_07_04.84Z.image_lev1.fits

[!NOTE] Time estimate: Some of these changes would be quite easy, but if we try to specify everything it could take much longer, 4-16 hours.

[mgalloy]

Looks like we probably need CTYPE1 and CTYPE2 keywords and should change CRUNIT1/CRUNIT2 to CUNIT1/CUNIT2. See page 31 of the FITS standard v4.0.

WCSAXES – [integer; default: NAXIS, or larger of WCS indices i or j]. Number of axes in the WCS description. This keyword, if present, must precede all WCS keywords except NAXIS in the HDU. The value of WCSAXES may exceed the number of pixel axes for the HDU.

CTYPEi – [string; indexed; default: '␣' (i.e. a linear, undefined axis)]. Type for the Intermediate-coordinate Axis i. Any coordinate type that is not covered by this Standard or an officially recognized FITS convention shall be taken to be linear. All non-linear coordinate system names must be expressed in ‘4–3’ form: the first four characters specify the coordinate type, the fifth character is a hyphen (‘-’), and the remaining three characters specify an algorithm code for computing the world coordinate value. Coordinate types with names of fewer than four characters are padded on the right with hyphens, and algorithm codes with fewer than three characters are padded on the right with blanks11. Algorithm codes should be three characters.

CUNITi – [string; indexed; default: '␣' (i.e., undefined)]. Physical units of CRVAL and CDELT for Axis i. Note that units should always be specified (see Sect. 4.3). Units for celestial coordinate systems defined in this Standard must be degrees.

CRPIXj – [floating point; indexed; default: 0.0]. Location of the reference point in the image for Axis j corresponding to r j in Eq. (9). Note that the reference point may lie outside the image and that the first pixel in the image has pixel coordinates (1.0, 1.0, . . .).

CRVALi – [floating point; indexed; default: 0.0]. World-coordinate value at the reference point of Axis i.

CDELTi – [floating point; indexed; default: 1.0]. Increment of the world coordinate at the reference point for Axis i. The value must not be zero.

CROTAi – [floating point; indexed; default: 0.0]. The amount of rotation from the standard coordinate system to a different coordinate system. Further use of this of this keyword is deprecated, in favor of the newer formalisms that use the CDi j or PCi j keywords to define the rotation.

PCi_j – [floating point; defaults: 1.0 when i = j, 0.0 otherwise]. Linear transformation matrix between Pixel Axes j and Intermediate-coordinate Axes i. The PCi j matrix must not be singular.

CDi_j – [floating point; defaults: 0.0, but see below]. Linear transformation matrix (with scale) between Pixel Axes j and Intermediate-coordinate Axes i. This nomenclature is equivalent to PCi j when CDELTi is unity. The CDi j matrix must not be singular. Note that the CDi j formalism is an exclusive alternative to PCi j, and the CDi j and PCi j keywords must not appear together within an HDU.

[mgalloy]

From [Thompson]:

In principal, any standard coordinate system could be used to specify the position of the observer, but it is recommended that at least the keywords DSUN_OBS, HGLN_OBS, and HGLT_OBS be included to specify the Stonyhurst heliographic coordinates of the observer.

[Thompson]: http://jsoc.stanford.edu/relevant_papers/ThompsonWT_coordSys.pdf "Coordinate systems for solar image data "

[mgalloy]

Commit 3f764f6 fixed errors for the first Map display, but I do get warnings still:

WARNING: SunpyMetadataWarning: Missing metadata for observer: assuming Earth-based observer.
For frame 'heliographic_stonyhurst' the following metadata is missing: hgln_obs,dsun_obs,hglt_obs
For frame 'heliographic_carrington' the following metadata is missing: crlt_obs,crln_obs,dsun_obs
 [sunpy.map.mapbase]

I still need to "project an AIA Sunpy Map to the UCoMP WCS" cell in the notebook. To do that step, I need the file:

aia.lev1.193A_2022-04-07T18_07_04.84Z.image_lev1.fits

[bberkeyU]

@mgalloy I have poked around the AIA sites and played with Sunpys JSOC FIDO service, and I cant figure out how to find or download a file named: aia.lev1.193A_2022-04-07T18_07_04.84Z.image_lev1.fits

It doesn't look like AIA created a 193 file at 18:07:04:84 UTC. I have found some 193 AIA image_lev1.fits files which I think should satisfy the test of co-alignment between UCoMP and AIA. Files are available in /home/berkey/aia on the unix hosts.

[mgalloy]

OK, I used one your files to get farther in the Jupyter notebook. I think we need the Earth-based observing location keywords to get through it, though.

[jburkepile]

Looking through Bill Thompson's paper and using sun.pro, I think we need these values for the 3 keywords Mike has listed about:

• [x] DSUN_OBS = distance to the Sun from observer in meters. Get this from sun.pro with the 'DIST' keyword (or variable) and convert to meters. 1 AU = 149597870700 meters
• [x] HGLN_OBS = Stonyhurst heliographic longitude - the observer position. Set HGLN_OBS = 0 (zero) to indicate location to be at Earth
• [x] HGLT_OBS = heliosgraphic latitude - as seen by the observer. Set this equal to the B-angle which we already get from the ephemeris.

[mgalloy]

OK, I am getting farther in the notebook, but there are still warnings and issues.

I still get the following warnings:

>>> aia_map_reproject = aia_map.reproject_to(WCS(header_1074))
WARNING: FITSFixedWarning: The WCS transformation has more axes (2) than the image it is associated with (0) [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 59823.767727 from DATE-OBS.
Set MJD-END to 59823.769417 from DATE-END'. [astropy.wcs.wcs]

followed by a long stack trace, ending with:

ValueError: Need to specify shape_out when specifying output_projection as WCS object

[jburkepile]

@mgalloy This warning: "The WCS transformation has more axes (2) than the image it is associated with (0)" makes me think the code is reading in the primary header which has zero axes. Can you run the code reading in the extension header? If you are already reading in the extension then disregard my comment

[mgalloy]

The extension header doesn't have any WCS information in it though, we need to use the primary header for that.

[mgalloy]

OK, using the workaround in the notebook:

fig = plt.figure(figsize=(8,8), constrained_layout=True)
norm_aia = ImageNormalize(aia_map_reproject.data, stretch=AsinhStretch(), vmin=0)
ax = fig.add_subplot(projection=ucomp_map_0)
ucomp_map_0.plot(axes=ax, cmap='sdoaia193', norm=norm_1074)
aia_map_reproject.plot(axes=ax, cmap='sdoaia193', norm=norm_aia)

I get the following image:

So something is still wrong with our WCS info. The notebook got the following image with the workaround:

[jburkepile]

MIke: Let's work on this at our regular weekly meeting tomorrow (Wed.) morning. Would be good to be able to complete this and close the ticket.

[mgalloy]

I wasn't using the same date for the UCoMP and AIA images. Now it works!

[mgalloy]

The keyword inheritance (INHERIT set to 'T') would make the notebook doing the above simpler, but evidently Astropy does not support it.

[bberkeyU]

How do we calculate the keywords in the "COMMENT --- World Coordinate System (WCS) info ---" section?

Are we using a particular program (like sun.pro) that we should document as a COMMENT like we do in the Ephemeris info section?

[jburkepile]

For definitions of world coordinate system see this paper: https://doi.org/10.1051/0004-6361:20054262

Coordinate systems for solar image data by W.T. Thompson, 2006

[bberkeyU]

But is there a software package we rely on to tell us the distance to the sun right now?

For the Ephemeris info (ie SOLAR_P0, CAR_ROT ) we document that we use sun.pro to calulate those values.

[mgalloy]

sun.pro is giving the distance to sun in AU, which is used to compute DSUN_OBS, and the value for HGLT_OBS (same as SOLAR_B0). I think those are the only two WCS values that use ephemeris information.

[mgalloy]

OK, there is a comment for the WCS section now too.