"""
=============
Map Histogram
=============

How to inspect the histogram of the data of a map.
"""
import numpy as np
import matplotlib.pyplot as plt

import astropy.units as u

import sunpy.map
from astropy.coordinates import SkyCoord
from sunpy.data.sample import AIA_171_IMAGE

###############################################################################
# We first create the Map using the sample data and we will create a submap
# of a quiet region.
aia = sunpy.map.Map(AIA_171_IMAGE)
bl = SkyCoord(-400 * u.arcsec, 0 * u.arcsec, frame=aia.coordinate_frame)
tr = SkyCoord(0 * u.arcsec, 400 * u.arcsec, frame=aia.coordinate_frame)
aia_smap = aia.submap(bl, tr)

###############################################################################
# We now create a histogram of the data in this region.
dmin = aia_smap.min()
dmax = aia_smap.max()
num_bins = 50
hist, bins = np.histogram(aia_smap.data, bins=np.linspace(dmin, dmax, num_bins))
width = 0.7 * (bins[1] - bins[0])
x = (bins[:-1] + bins[1:]) / 2

###############################################################################
# Let's plot the histogram as well as some standard values such as mean
# upper, and lower value and the one-sigma range.
plt.figure()
plt.bar(x, hist, align='center', width=width, label='Histogram')
plt.xlabel('Intensity')
plt.axvline(dmin, label='Data min={:.2f}'.format(dmin), color='black')
plt.axvline(dmax, label='Data max={:.2f}'.format(dmax), color='black')
plt.axvline(aia_smap.data.mean(),
            label='mean={:.2f}'.format(aia_smap.data.mean()), color='green')
one_sigma = np.array([aia_smap.data.mean() - aia_smap.data.std(),
                      aia_smap.data.mean() + aia_smap.data.std()])
plt.axvspan(one_sigma[0], one_sigma[1], alpha=0.3, color='green',
            label='mean +/- std = [{0:.2f}, {1:.2f}]'.format(
            one_sigma[0], one_sigma[1]))
plt.axvline(one_sigma[0], color='green')
plt.axvline(one_sigma[1], color='red')
plt.legend()
plt.show()

###############################################################################
# Finally let's overplot what the one-sigma range means on the map
fig = plt.figure()
fig.add_subplot(projection=aia_smap)
aia_smap.plot()
levels = one_sigma / dmax * u.percent * 100
aia_smap.draw_contours(levels=levels, colors=['red', 'green'])
plt.show()

# -*- coding: utf-8 -*-
"""
=====================
Histograming map data
=====================

How to inspect the histogram of the data of a map.
"""
import numpy as np
import matplotlib.pyplot as plt

import astropy.units as u

import sunpy.map
from astropy.coordinates import SkyCoord
from sunpy.data.sample import AIA_171_IMAGE

###############################################################################
# We start with the sample data and create a cutout.
aia = sunpy.map.Map(AIA_171_IMAGE)
bottom_left = SkyCoord(-300 * u.arcsec, 0 * u.arcsec, frame=aia.coordinate_frame)
top_right = SkyCoord(100 * u.arcsec, 400 * u.arcsec, frame=aia.coordinate_frame)
aia_smap = aia.submap(bottom_left, top_right)

###############################################################################
# The image of a `~sunpy.map.GenericMap` is always available in the data attribute.
# Map also provides shortcuts to the image minimum and maximum values.
# Let's create a histogram of the data in this submap.
num_bins = 50
bins = np.linspace(aia_smap.min(), aia_smap.max(), num_bins)
hist, bins = np.histogram(aia_smap.data, bins=bins)
width = 0.7 * (bins[1] - bins[0])
x = (bins[:-1] + bins[1:]) / 2

###############################################################################
# Let's plot the histogram as well as some standard values such as mean
# upper, and lower value and the one-sigma range.
plt.figure()
plt.bar(x, hist, align='center', width=width, label='Histogram')
plt.xlabel('Intensity')
plt.axvline(aia_smap.min(), label='Data min={:.2f}'.format(aia_smap.min()), color='black')
plt.axvline(aia_smap.max(), label='Data max={:.2f}'.format(aia_smap.max()), color='black')
plt.axvline(aia_smap.data.mean(),
            label='mean={:.2f}'.format(aia_smap.data.mean()), color='green')
one_sigma = np.array([aia_smap.data.mean() - aia_smap.data.std(),
                      aia_smap.data.mean() + aia_smap.data.std()])
plt.axvspan(one_sigma[0], one_sigma[1], alpha=0.3, color='green',
            label='mean +/- std = [{0:.2f}, {1:.2f}]'.format(
            one_sigma[0], one_sigma[1]))
plt.axvline(one_sigma[0], color='green')
plt.axvline(one_sigma[1], color='red')
plt.legend()
plt.show()

###############################################################################
# Finally let's overplot the one-sigma contours
fig = plt.figure()
fig.add_subplot(projection=aia_smap)
aia_smap.plot()
levels = one_sigma / aia_smap.max() * u.percent * 100
aia_smap.draw_contours(levels=levels, colors=['blue'])
plt.show()

Update map_data_histogram.py