Improve IAU/COSMO documentation in enmap.read_map

[zatkins2]

• pixell version: 0.12.1
• Python version: 3.9.5
• Operating System: Ubuntu 20.04

Description

Currently enmap.read_map flips the U component of an ndmap in an opaque way, depending on the contents of the fits header. The fits header needs to have specific cards/values for the introspection to flip the component. These conditions (all evaluated in enmap.get_stokes_flips) should be documented. Also, the default behavior of the function (to convert to COSMO) should be documented, or amended and documented if we want to run with IAU going forward.

Even better, it would be nice to have some user arguments supply optional information about the data, e.g. which axes might need to be flipped, and what the desired output convention is, which could help handle the case of "simple" fits headers (as are produced with the enmap.write_map defaults) when the automatic introspection of the polarization convention, according to the conditions, might fail.

What I Did

Playing with fits header of a file with a known polarization convention. I believe all these conditions (see enmap.get_stokes_flips) must be true for enmap.read_fits to flip the polarization convention:

1. There must be a "CTYPEn" card with value "STOKES", where "n" is the "NAXIS" index of each stokes-like component, for each stokes-like component.
2. If included in the header, "CRPIXn", "CRVALn", and "CRDELTn" (with "n" as above) must be set so that (3-crval)/cdelt+crpix - 1 gives the 0-indexed index of the U component in "NAXISn".
3. If "CRPIXn", "CRVALn", and "CRDELTn" are not included in the header, the U component in "NAXISn" must be at index 2 of the 0-indexed axis (i.e., their default values are all 1.0 such that (3-crval)/cdelt+crpix - 1 is 2).
4. A "POLCONV" or "POLCCONV" card must be supplied with either "COSMO" or "IAU". If previous conditions met and either card is not supplied, assume data in IAU and convert to COSMO. If previous conditions met and the card value is not "COSMO" or "IAU," assume data in IAU and do no conversion.

I agree there is no clean way of avoiding the need for specific set of cards/values for automatic introspection, but the behavior should be made clear for new users I think.

[zatkins2]

If stakeholders agree on above, I am happy to submit a PR with documentation/more user-supplied arguments, e.g. the desired output convention. However, I think it's worth discussing if the default behavior -- to try to put the data in the COSMO convention -- should be changed to IAU?

[amaurea]

It's not arbitrary which representation the read-in map ends up with. Several hamonic conversion functions break if they do not get data in the convention they expect, which is COSMO. Changing these functions would break lots of existing code.

I am not very keen on adding arguments to all those functions to make them support both conventions. It is much cleaner to convert once on input and then work with a consistent convention than it is to need to worry about which particular convention some particular map somewhere in the code has.

So I'm against changing the default output from COSMO. I'm OK with adding an argument to the reading function to control this, with options like "cosmo", "iau" and "raw". But I don't think it will be useful to most people, as they will need to convert to cosmo to do harmonic operations.

Documenting the behavior is a good idea, though.

If previous conditions met and the card value is not "COSMO" or "IAU," assume data in IAU and do no conversion.

That's not right, is it? If POLCCONV has an unrecognized value, then it's assumed to be COSMO, and therefore doesn't need to be converted.

[zatkins2]

Hey Sigurd, thanks for taking a look. I agree it's not arbitrary, it's just that it's not as clear as it could be. Definitely agree that adding arguments to harmonic operations is not the right approach, but do think allowing more user control/visibility just on reading data from disk would be good. I think with documentation, we can trust users to do the right thing for their application. I'm OK with keeping COSMO as the default, wasn't aware of the depth of the dependencies on this!

In short, specifically, I propose: (1) documenting the behavior, (2) adding an argument with default "cosmo" to read_map, and your other possible options sound good, (3) increasing verbosity of the output and setting verbose=True as default, which would mimic healpy (basically, have it print by default how the convention was interpreted, whether it was changed, and to what).

That's not right, is it? If POLCCONV has an unrecognized value, then it's assumed to be COSMO, and therefore doesn't need to be converted.

Sorry, correct, should have said "assume data in COSMO and do no conversion."

If the above suggestions look OK, I can submit a PR. Thanks!

[zatkins2]

I'd like to add to this issue with a suggestion and question:

Suggestion

write_map should:

1. be similarly documented so that people know what they need to do to get polarization convention info into the header (currently, no convention info, nor whether there are any STOKES axes, are written by default).
2. Even better, there should be some simple wrapper that just takes the desired output convention (by default "cosmo") and does the flipping/header stuff for you, so that it is easier for users to do what they want. This would also have the (good) effect of recording necessary convention info in the header always, even when it is the default, uninteresting, all-COSMO case.
3. Even even better, I think the "system" for recording polarization convention info in the header isn't flexible enough. It can't easily handle if you have some polarization matrix (e.g. the 3x3 cousin of ivar maps, a.k.a. div maps), where only some elements of a given axis need to be flipped. I think the best way would be to have an auxiliary array of data be stored that explicitly records which elements need flipping in an axis, possibly more than 1. Since this would break old data files, there can be an additional keyword added to the header that indicates which "flipping implementation" is in the file, which defaults to the old implementation if it isn't present. read_map would then need to be updated accordingly.

Question (@msyriac)

This issue is pertinent to the SSIA specification that LAT products are saved in IAU but always loaded by pixell into COSMO. From what @amaurea has said, this is really necessary since a lot of cosmology code depends on data being in the COSMO convention. If that's true, I'm curious what the utility is of even saving maps to disk in IAU? Is there any software that can interact with our maps other than pixell?

[msyriac]

The motivation was that this is an official recommendation from IAU, one that astronomers outside the cosmological community seem to adhere to. We use FITS and WCS, so someone could use astropy or similar (e.g. https://docs.astropy.org/en/stable/generated/examples/io/plot_fits-image.html ) to load our maps. They may then potentially incorrectly interpret our polarization maps if we have them save them to disk using COSMO.

[zatkins2]

Ah right, nice. Sorry, had it in my mind that someone might never see an “IAU SO map” if they were stuck with pixell, but yeah clearly that’s not true :)

Anyway, my bigger point still stands: read_map and write_map need a little work to be more explicit and safe (and I can be the one to propose new code! Just want to be sure folks agree)