Trouble making a map with SOLIS Full-disk Magnetograms

[amunozj]

Description

Hi guys! thank you so much for all the work that you do! I'm using SunPy to load a bunch of different full-disk magnetograms and I have been able to make it work in all cases (HMI, MDI, Mt. Wilson, KPVT, SOT) But not for SOLIS. In all cases, the issues I had to tackle were Header fields that were named differently or were missing, but for this one I'm at a loss.

Expected behavior

I'm making the maps by loading the fits files, select the layer/cube slice and header and assembling the maps using the map function. I expect it to work like the other instruments

Actual behavior

I get the error:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-42-1bfc7962b6e7> in <module>
     11 
     12 # Assemble Sunpy map
---> 13 SOLISmap = sunpy.map.Map(dataSOLIS, SOLIS_fits[0].header)
     14 
     15 # # # Fix instrument rotation

C:\ProgramData\Anaconda3\lib\site-packages\sunpy\map\map_factory.py in __call__(self, composite, sequence, silence_errors, *args, **kwargs)
    297 
    298         if not len(new_maps):
--> 299             raise RuntimeError('No maps loaded')
    300 
    301         # If the list is meant to be a sequence, instantiate a map sequence

RuntimeError: No maps loaded

Steps to Reproduce

Here is what I'm doing for SOLIS:

from astropy.io import fits
import sunpy.map

SOLISfile = 'Data/SOLIS/k4v93110401t161025_FDISK.fts.gz'

SOLIS_fits = fits.open(SOLISfile, cache=False)
SOLIS_fits.verify('fix')
dataSOLIS = SOLIS_fits[0].data[6,:,:]
SOLISmap = sunpy.map.Map(dataSOLIS, SOLIS_fits[0].header)

The file can be donwloaded here:

https://solis.nso.edu/pubkeep/v93/201104/k4v93110401/k4v93110401t161025_FDISK.fts.gz

System Details

• SunPy Version: 1.0.1
• Python Version: 3.7
• OS information: Windows

[amunozj]

A quick update. The problem is with the header. If I use a random header from another instrument there is no problem. The problem is that I'm not sure what is the problem with the header :(

[Cadair]

👋 @amunozj Really glad to hear that SunPy is being useful for you. If you have code that is making those other data sources work, or examples of using them we would love to have those contributed.

The main error here seems to be the fact we enforce the presence of the unit keywords in 1.0. You should have got a warning log line explaining why the header failed to validate.

In [5]: sunpy.map.Map("k4v93110401t161025_FDISK.fts.gz")                                                                                                                                                                           
WARNING: SunpyUserWarning: This file contains more than 2 dimensions. Data will be truncated to the first two dimensions. [sunpy.map.mapbase]
/home/stuart/Git/sunpy/sunpy/map/map_factory.py:292: UserWarning: One of the data, header pairs failed to validate with: Image coordinate units for axis 1 not present in metadata.
Image coordinate units for axis 2 not present in metadata.
See https://docs.sunpy.org/en/stable/code_ref/map.html#fixing-map-metadata` for instructions on how to add missing metadata.

I have got something that looks approximately correct using this:

In [30]: data, header = sunpy.io.read_file("k4v93110401t161025_FDISK.fts.gz")[0]                                                                                                                                                   

In [31]: header['CUNIT1'] = "arcsec"; header['CUNIT2'] = "arcsec"                                                                                                                                                                  

In [32]: sm = sunpy.map.Map(data, header)                                                                                                                                                                                          
WARNING: SunpyUserWarning: This file contains more than 2 dimensions. Data will be truncated to the first two dimensions. [sunpy.map.mapbase]

In [33]: sm.peek()                                                                                                                                                                                                                 
WARNING: SunpyUserWarning: Missing metadata for observer: assuming Earth-based observer.The following sets of keys were checked:
('hgln_obs', 'hglt_obs', 'dsun_obs')
('crln_obs', 'crlt_obs', 'dsun_obs') [sunpy.map.mapbase]

You will note that it's only creating a Map for the first slice of the 3rd (FITS) dimension. If you want the other dimensions you will need to slice the map yourself:

In [32]: sm = sunpy.map.Map(data[1], header)

or if you want all 12 you can do this:

smaps = sunpy.map.Map(list(zip(data, [header]*data.shape[0])))

Let me know if you have any other questions or if I can help you in decoding the header or adding support for the other data sources you are using.

[amunozj]

Thank you so much @Cadair ! I'm an idiot, I had turned off warnings and missed the warning. I had already solved for another instrument. Here is a jupyter notebook with the code I did to use map for all these magnetographs. It is working great!

https://nbviewer.jupyter.org/github/amunozj/Magnetograms/blob/master/Magnetogram_Comparison.ipynb

[Cadair]

Would you be interested in adding Map classes for some of these data? It would be great if the header mapping was in SunPy.

[amunozj]

Unfortunately I'm over-committed for the summer, so it would be some time before I can even start thinking about it. If my commitments lessen I will let you know.