def running_difference(mc, offset=1, use_offset_for_meta='mean',
                       image_normalize=True):
    """
    Calculate the running difference of a mapcube.

    Parameters
    ----------
    mc : sunpy.map.MapCube
       A sunpy mapcube object

    offset : [ int ]
       Calculate the running difference between map 'i + offset' and image 'i'.

    use_offset_for_meta : {'ahead', 'behind', 'mean'}
       Which meta header to use in layer 'i' in the returned mapcube, either
       from map 'i + offset' (when set to 'ahead') and image 'i' (when set to
       'behind').  When set to 'mean', the ahead meta object is copied, with
       the observation date replaced with the mean of the ahead and behind
       observation dates.

    image_normalize : bool
        If true, return the mapcube with the same image normalization applied
        to all maps in the mapcube.

    Returns
    -------
    sunpy.map.MapCube
       A mapcube containing the running difference of the input mapcube.
       The value normalization function used in plotting the data is changed,
       prettifying movies of resultant mapcube.
    """
    # Create a list containing the data for the new map object
    new_mc = []
    for i in range(0, len(mc.maps) - offset):
        new_data = mc[i + offset].data - mc[i].data
        if use_offset_for_meta == 'ahead':
            new_meta = mc[i + offset].meta
            plot_settings = mc[i + offset].plot_settings
        elif use_offset_for_meta == 'behind':
            new_meta = mc[i].meta
            plot_settings = mc[i].plot_settings
        elif use_offset_for_meta == 'mean':
            new_meta = deepcopy(mc[i + offset].meta)
            new_meta['date_obs'] = _mean_time([parse_time(mc[i + offset].date),
                                               parse_time(mc[i].date)])
            plot_settings = mc[i + offset].plot_settings
        else:
            raise ValueError('The value of the keyword "use_offset_for_meta" has not been recognized.')

        # Update the plot scaling.  The default here attempts to produce decent
        # looking images
        new_map = Map(new_data, new_meta)
        new_map.plot_settings = plot_settings
        new_mc.append(new_map)

    # Create the new mapcube and return
    if image_normalize:
        return movie_normalization(Map(new_mc, cube=True), stretch=LinearStretch())
    else:
        return Map(new_mc, cube=True)

def persistence(mc, func=np.max, image_normalize=True):
    """
    Parameters
    ----------
    mc : sunpy.map.MapCube
       A sunpy mapcube object

    Returns
    -------
    sunpy.map.MapCube
       A mapcube containing the persistence transform of the input mapcube.
       The value normalization function used in plotting the data is changed,
       prettifying movies of resultant mapcube.
    """

    # Get the persistence transform
    new_datacube = persistence_dc(mc.as_array(), func=func)

    # Create a list containing the data for the new map object
    new_mc = []
    for i, m in enumerate(mc):
        new_map = Map(new_datacube[:, :, i], m.meta)
        new_map.plot_settings = deepcopy(m.plot_settings)
        new_mc.append(new_map)

    # Create the new mapcube and return
    if image_normalize:
        return movie_normalization(Map(new_mc, cube=True))
    else:
        return Map(new_mc, cube=True)

def accumulate(mc, accum, normalize=True):
    """
    Parameters
    ----------
    mc : sunpy.map.MapCube
       A sunpy mapcube object

    accum :

    normalize :

    Returns
    -------
    sunpy.map.MapCube
       A summed mapcube in the map layer (time) direction.

    """

    # counter for number of maps.
    j = 0

    # storage for the returned maps
    maps = []
    nmaps = len(mc)

    while j + accum <= nmaps:
        i = 0
        these_map_times = []
        while i < accum:
            this_map = mc[i + j]
            these_map_times.append(parse_time(this_map.date))
            if normalize:
                normalization = this_map.exposure_time
            else:
                normalization = 1.0
            if i == 0:
                # Emission rate
                m = this_map.data / normalization
            else:
                # Emission rate
                m += this_map.data / normalization
            i += 1
        j += accum
        # Make a copy of the meta header and set the exposure time to accum,
        # indicating that 'n' normalized exposures were used.
        new_meta = deepcopy(this_map.meta)
        new_meta['exptime'] = np.float64(accum)

        # Set the observation time to the average of the times used to form
        # the map.
        new_meta['date_obs'] = _mean_time(these_map_times)

        # Create the map list that will be used to make the mapcube
        new_map = Map(m, new_meta)
        new_map.plot_settings = deepcopy(this_map.plot_settings)
        maps.append(new_map)

    # Create the new mapcube and return
    return Map(maps, cube=True)

def test2():
    aia = Map(sunpy.AIA_171_IMAGE)
    fig = plt.figure()
    ax = plt.subplot(111)
    aia.plot()
    plt.colorbar()
    aia.draw_limb()
    plt.show()

def create_tempmap(date, n_params=1, data_dir=home+'SDO_data/',
                   maps_dir=home+'temperature_maps/'):
    wlens = ['94', '131', '171', '193', '211', '335']
    t0 = 5.6
    images = []
    #imdates = {}
    
    print 'Finding data for {}.'.format(date.date())
    # Loop through wavelengths
    for wl, wlen in enumerate(wlens):
        #print 'Finding {}A data...'.format(wlen),
        fits_dir = data_dir + '{}/{:%Y/%m/%d}/'.format(wlen, date)
        filename = fits_dir + 'aia*{0}*{1:%Y?%m?%d}?{1:%H?%M}*lev1?fits'.format(wlen, date)
        temp_im = Map(filename)
        # Download data if not enough found
        client = vso.VSOClient()
        if temp_im == []:
            print 'File not found. Downloading from VSO...'
            # Wavelength value for query needs to be an astropy Quantity
            wquant = u.Quantity(value=int(wlen), unit='Angstrom')
            qr = client.query(vso.attrs.Time(date,# - dt.timedelta(seconds=6),
                                             date + dt.timedelta(seconds=12)),#6)),
                              vso.attrs.Wave(wquant, wquant),
                              vso.attrs.Instrument('aia'),
                              vso.attrs.Provider('JSOC'))
            res = client.get(qr, path=fits_dir+'{file}', site='NSO').wait()
            temp_im = Map(res)
        if temp_im == []:
            print 'Downloading failed.'
            print res, len(qr), qr
            return np.zeros((512, 512)), None, None
        if isinstance(temp_im, list):
            temp_im = temp_im[0]
        # TODO: save out level 1.5 data so it can be loaded quickly.
        temp_im = aiaprep(temp_im)
        temp_im.data = temp_im.data / temp_im.exposure_time # Can probably increase speed a bit by making this * (1.0/exp_time)
        images.append(temp_im)
        #imdates[wlen] = temp_im.date
    
    normim = images[2].data.copy()
    # Normalise images to 171A
    print 'Normalising images'
    for i in range(len(wlens)):
        images[i].data = images[i].data / normim
    
    # Produce temperature map
    if n_params == 1:
        tempmap = find_temp(images, t0)#, force_temp_scan=True)
    else:
        #tempmap = find_temp_3params(images, t0)
        pass

    return tempmap

 def save(self):
     date = sunpy.time.parse_time(self.date)
     if not path.exists(self.maps_dir):
         makedirs(self.maps_dir)
     fname = path.join(self.maps_dir,
                       '{:%Y-%m-%dT%H_%M_%S}.fits'.format(date))
     alldata = np.zeros((self.shape[0], self.shape[1], self.n_params+1))
     alldata[..., 0] = self.data
     if self.n_params != 1:
         fname = fname.replace('.fits', '_full.fits')
         alldata[..., 1] = self.dem_width
         alldata[..., 2] = self.emission_measure
     alldata[..., -1] = self.goodness_of_fit
     outmap = Map(alldata, self.meta.copy())
     outmap.save(fname, clobber=True)

    def get_properties(cls, header):
        """Parses SXT image header"""
        properties = Map.get_properties(header)

        # 2012/12/19 - the SXT headers do not have a value of the distance from
        # the spacecraft to the center of the Sun.  The FITS keyword 'DSUN_OBS'
        # appears to refer to the observed diameter of the Sun.  Until such
        # time as that is calculated and properly included in the file, we will
        # use simple trigonometry to calculate the distance of the center of
        # the Sun from the spacecraft.  Note that the small angle approximation
        # is used, and the solar radius stored in SXT FITS files is in arcseconds.
        properties["dsun"] = constants.au
        yohkoh_solar_r = header.get("solar_r", None)
        if yohkoh_solar_r == None:
            properties["dsun"] = constants.au
        else:
            properties["dsun"] = constants.radius / (np.deg2rad(yohkoh_solar_r / 3600.0))

        wavelnth = header.get("wavelnth")
        if wavelnth == "Al.1":
            wavelnth = "Al01"
        if wavelnth.lower() == "open":
            wavelnth = "white light"

        properties.update(
            {
                "detector": "SXT",
                "instrument": "SXT",
                "observatory": "Yohkoh",
                "name": "SXT %s" % wavelnth,
                "nickname": "SXT",
                "cmap": cm.get_cmap(name="yohkohsxt" + wavelnth[0:2].lower()),
            }
        )
        return properties

    def get_properties(cls, header):
        """Parses SXT image header"""
        properties = Map.get_properties(header)
        
        # 2012/11/07 - the SXT headers do not have a value of the distance from
        # the spacecraft to the center of the Sun.  The FITS keyword 'DSUN_OBS'
        # appears to refer to the observed diameter of the Sun.  Until such 
        # time as that is calculated and properly included in the file, we will 
        # use the value of 1 AU as a standard.
        properties['dsun']= constants.au
        
        wavelnth = header.get('wavelnth')
        if wavelnth == 'Al.1':
            wavelnth = 'Al01'
        if wavelnth.lower() == 'open':
            wavelnth = 'white light'

        properties.update({
            "detector": "SXT",
            "instrument": "SXT",
            "observatory": "Yohkoh",
            "name": "SXT %s" % wavelnth,
            "nickname": "SXT",
            "cmap": cm.get_cmap(name='yohkohsxt' + wavelnth[0:2].lower())
        })
        return properties 

def submap(mc, range_a, range_b, **kwargs):
    """
    Parameters
    ----------
    mc : sunpy.map.MapCube
       A sunpy mapcube object

    range_a : list


    range_b : list

    Returns
    -------
    sunpy.map.MapCube
       A mapcube containing maps that have had the map submap
       method applied to each layer.
    """
    nmc = len(mc)
    if (len(range_a) == nmc) and (len(range_b) == nmc):
        ra = range_a
        rb = range_b
    elif (len(range_a) == 1) and (len(range_b) == 1):
        ra = [range_a for i in range(0, nmc)]
        rb = [range_b for i in range(0, nmc)]
    else:
        raise ValueError('Both input ranges must be either of size 1 or size '
                         'equal to the number of maps in the mapcube')

    # Storage for the returned maps
    maps = []
    for im, m in enumerate(mc):
        maps.append(Map.submap(m, ra[im], rb[im], **kwargs))
    # Create the new mapcube and return
    return Map(maps, cube=True)

    def add_map(self, input_, zorder=None, alpha=1, levels=False):
        """Adds a map to the CompositeMap
        
        Parameters
        ----------
        input_ : {sunpy.map, string}
            Map instance or filepath to map to be added
        zorder : int
            The index to use when determining where the map should lie along
            the z-axis; maps with higher z-orders appear above maps with lower
            z-orders.
        alpha : float
            Opacity at which the map should be displayed. An alpha value of 0
            results in a fully transparent image while an alpha value of 1
            results in a fully opaque image. Values between result in semi-
            transparent images.

        """
        if zorder is None:
            zorder = max([m.zorder for m in self._maps]) + 10
        
        m = Map.read(input_)
        m.zorder = zorder
        m.alpha = alpha
        m.levels = levels
        
        self._maps.append(m)

    def __new__(cls, *args, **kwargs):
        """Creates a new Map instance"""
        
        maps = []
        data = []
        headers = []
    
        # convert input to maps
        for item in args:
            if isinstance(item, Map):
                maps.append(item)
            else:
                maps.append(Map.read(item))

        # sort data
        sortby = kwargs.get("sortby", "date")
        if hasattr(cls, '_sort_by_%s' % sortby):
            maps.sort(key=getattr(cls, '_sort_by_%s' % sortby)())

        # create data cube
        for map_ in maps:
            data.append(np.array(map_))
            headers.append(map_._original_header)

        obj = np.asarray(data).view(cls)
        obj._headers = headers

        return obj

    def get_properties(cls, header):
        """Parses XRT image header"""
        properties = Map.get_properties(header)
        # XRT uses DATE_OBS, not date-obs.
        properties["date"] = parse_time(header.get('date_obs', None))

        #TODO: proper exception handling here - report to the user that there is
        # an unexpected value
        fw1 = header.get('EC_FW1_')
        if not(fw1.lower() in [x.lower() for x in cls.filter_wheel1_measurements]):
            pass
        fw2 = header.get('EC_FW2_')
        if not(fw2.lower() in [x.lower() for x in cls.filter_wheel2_measurements]):
            pass

        # All images get the same color table - IDL Red temperature (loadct, 3)
        properties.update({
            "detector": "XRT",
            "instrument": "XRT",
            "observatory": "Hinode",
            "name": "XRT %s-%s " % (fw1.replace('_', ' '),
                                       fw2.replace('_', ' ')),
            "nickname": "XRT",
            "cmap": cm.get_cmap(name='hinodexrt')
        })
        return properties

def add_noise(params, wave_maps, verbose=False):
    """
    Adds simulated noise to a list of maps
    """
    wave_maps_noise = []
    for current_wave_map in wave_maps:
        if verbose:
            print("  * Adding noise to map at " + str(current_wave_map.date))

        noise = noise_random(params, current_wave_map.data.shape)
        struct = noise_structure(params, current_wave_map.data.shape)

        noisy_wave_map = Map(current_wave_map.data + noise + struct,
                                       current_wave_map.meta)
        noisy_wave_map.plot_settings = deepcopy(current_wave_map.plot_settings)
        wave_maps_noise.append(noisy_wave_map)

    return Map(wave_maps_noise, cube=True)

def clean(params, wave_maps, verbose=False):
    """
    Cleans a list of maps
    """
    wave_maps_clean = []
    for current_wave_map in wave_maps:
        if verbose:
            print("  * Cleaning map at "+str(current_wave_map.date))

        data = np.asarray(current_wave_map.data)
        if params.get("clean_nans"):
            data[np.isnan(data)] = 0.
                
        cleaned_wave_map = Map(data, current_wave_map.meta)
        # cleaned_wave_map.name = current_wave_map.name
        cleaned_wave_map.meta['date-obs'] = current_wave_map.date
        cleaned_wave_map.plot_settings = deepcopy(current_wave_map.plot_settings)
        wave_maps_clean.append(cleaned_wave_map)

    return Map(wave_maps_clean, cube=True)

    def __getitem__(self, key):
        """Overiding indexing operation"""
        if self.ndim is 3 and isinstance(key, int):
            data = np.ndarray.__getitem__(self, key)
            header = self._headers[key]
            for cls in Map.__subclasses__():
                if cls.is_datasource_for(header):
                    return cls(data, header)

        else:
            return np.ndarray.__getitem__(self, key)

 def get_properties(cls, header):
     """Parses EUVI image header"""
     properties = Map.get_properties(header)
     
     properties.update({
         "date": parse_time(header.get('date-obs',header.get('date_obs'))),
         "detector": "EUVI",
         "instrument": "SECCHI",
         "observatory": header.get('obsrvtry'),
         "cmap": cm.get_cmap('sohoeit%d' % header.get('wavelnth')),
         "nickname": "EUVI-" + header.get('obsrvtry')[-1]
     })
     return properties

 def get_properties(cls, header):
     """Parses SWAP image header"""
     properties = Map.get_properties(header)
     
     properties.update({
         "detector": "SWAP",
         "instrument": "SWAP",
         "observatory": "PROBA2",
         "name": "SWAP %s" % header.get('wavelnth'),
         "nickname": "SWAP",
         "cmap": cm.get_cmap(name='sdoaia171')
     })
     return properties

 def get_properties(cls, header):
     """Parses AIA image header"""
     properties = Map.get_properties(header)
     
     properties.update({
         "detector": "AIA",
         "instrument": "AIA",
         "observatory": "SDO",
         "nickname": "AIA",
         "cmap": cm.get_cmap('sdoaia%d' % header.get('wavelnth')),
         "processing_level": header.get('LVL_NUM')            
     })
     return properties

 def get_properties(cls, header):
     """Parses LASCO image header"""
     properties = Map.get_properties(header)
     
     datestr = "%sT%s" % (header.get('date_obs'), header.get('time_obs'))
     
     properties.update({
         "date": parse_time(datestr),
         "measurement": "white-light",
         "name": "LASCO %s" % header.get('detector'),
         "nickname": "LASCO-%s" % header.get('detector'),
         "cmap": cm.get_cmap('soholasco%s' % properties['detector'][1])
     })
     return properties

    def get_properties(cls, header):
        """Parses AIA image header"""
        properties = Map.get_properties(header)

        properties.update(
            {
                "detector": "AIA",
                "instrument": "AIA",
                "observatory": "SDO",
                "nickname": "AIA",
                "cmap": cm.get_cmap("sdoaia%d" % header.get("wavelnth")),
            }
        )
        return properties