Loading an HMI daily synoptic image into a Map

In this example we load the Daily Synoptic Maps produced by the HMI team. This data is an interesting demonstration of SunPy’s Map class as it is not in the more common Helioprojective coordinate system, it is in Heliographic Carrington coordinates and in a non-trivial Cylindrical Equal Area projection.

import numpy as np
import matplotlib.pyplot as plt

from astropy.utils.data import download_file

import sunpy.map

Let’s download the file and read it into a Map

filename = download_file(
    'http://jsoc.stanford.edu/data/hmi/synoptic/hmi.Synoptic_Mr.2191.fits', cache=True)
syn_map = sunpy.map.Map(filename)

There are a couple of oddities with this file, firstly the value of ‘CUNIT2’:

print(syn_map.meta['CUNIT2'])

Out:

Sine Latitude

That is not a unit! What this is telling us is that the latitude coordinate is actually the sine of latitude. According to the Thompson (2006) paper, CUNIT2 should be in degrees and CDELT2 should be multiplied by 180/pi. Also the value of CDELT1 has the wrong sign so let’s fix this.

syn_map.meta['CUNIT2'] = 'degree'
syn_map.meta['CDELT2'] = 180 / np.pi * syn_map.meta['CDELT2']
syn_map.meta['CDELT1'] *= -1

Let’s also fix the plot settings

syn_map.plot_settings['cmap'] = 'hmimag'
syn_map.plot_settings['norm'] = plt.Normalize(-1500, 1500)

Let’s plot the results

fig = plt.figure(figsize=(12, 5))
axes = plt.subplot(projection=syn_map)
im = syn_map.plot()

# Set up the Sine Latitude Grid
x = axes.coords[0]
y = axes.coords[1]

x.set_coord_type('longitude', coord_wrap=360.)
y.set_coord_type('latitude')

x.set_major_formatter('dd')
y.set_major_formatter('dd')

x.set_axislabel("Carrington Longitude [deg]")
y.set_axislabel("Latitude [deg]")

x.set_ticks(color='black', exclude_overlapping=True)
y.set_ticks(color='black', exclude_overlapping=True)

axes.coords.grid(color='black', alpha=0.6, linestyle='dotted', linewidth=0.5)

cb = plt.colorbar(im, fraction=0.019, pad=0.1)
cb.set_label("LOS Magnetic Field [Gauss]")

# Another horrible hack to make the ticks draw on the RHS
axes.set_ylim((1, syn_map.data.shape[0] - 1))
plt.title("{} {}-{}".format(syn_map.meta['content'], syn_map.meta['CAR_ROT'],
                            syn_map.meta['CAR_ROT'] + 1))
plt.show()

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/sunpy-1.0.3.dev1+gece6362cd-py3.7-linux-x86_64.egg/sunpy/map/mapbase.py:415: SunpyUserWarning: Missing metadata for observation time: setting observation time to current time.
  SunpyUserWarning)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/sunpy-1.0.3.dev1+gece6362cd-py3.7-linux-x86_64.egg/sunpy/map/mapbase.py:601: 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')
  warnings.warn(warning_message, SunpyUserWarning)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/sunpy-1.0.3.dev1+gece6362cd-py3.7-linux-x86_64.egg/sunpy/map/mapbase.py:601: 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')
  warnings.warn(warning_message, SunpyUserWarning)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/astropy/visualization/wcsaxes/coordinate_helpers.py:355: AstropyDeprecationWarning: exclude_overlapping= should be passed to set_ticklabel instead of set_ticks
  AstropyDeprecationWarning)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/astropy/visualization/wcsaxes/coordinate_helpers.py:355: AstropyDeprecationWarning: exclude_overlapping= should be passed to set_ticklabel instead of set_ticks
  AstropyDeprecationWarning)