Python Map.submap方法代码示例

# 需要导入模块: from sunpy.map import Map [as 别名]
# 或者: from sunpy.map.Map import submap [as 别名]
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)

# 需要导入模块: from sunpy.map import Map [as 别名]
# 或者: from sunpy.map.Map import submap [as 别名]
    def calculate_em(self, wlen='171', dz=100, model=False):
        """
        Calculate an approximation of the coronal EmissionMeasure using a given
        TemperatureMap object and a particular AIA channel.
    
        Parameters
        ----------
        tmap : CoronaTemps.temperature.TemperatureMap
            A TemperatureMap instance containing coronal temperature data
        wlen : {'94' | '131' | '171' | '193' | '211' | '335'}
            AIA wavelength used to approximate the emission measure. '171', '193'
            and '211' are most likely to provide reliable results. Use of other
            channels is not recommended.
        """
        # Load the appropriate temperature response function
        tresp = read('/imaps/holly/home/ajl7/CoronaTemps/aia_tresp')
        resp = tresp['resp{}'.format(wlen)]
    
        # Get some information from the TemperatureMap and set up filenames, etc
        tempdata = self.data.copy()
        tempdata[np.isnan(tempdata)] = 0.0
        date = sunpy.time.parse_time(self.date)
        if not model:
            data_dir = self.data_dir
            fits_dir = path.join(data_dir, '{:%Y/%m/%d}/{}'.format(date, wlen))
            filename = path.join(fits_dir,
                                 '*{0:%Y?%m?%d}?{0:%H?%M}*fits'.format(date))
            if wlen == '94': filename = filename.replace('94', '094')
    
            # Load and appropriately process AIA data
            filelist = glob.glob(filename)
            if filelist == []:
                print 'AIA data not found :('
                return
            aiamap = Map(filename)
            aiamap.data /= aiamap.exposure_time
            aiamap = aiaprep(aiamap)
            aiamap = aiamap.submap(self.xrange, self.yrange)
        else:
            fname = '/imaps/holly/home/ajl7/CoronaTemps/data/synthetic/{}/model.fits'.format(wlen)
            if wlen == '94': fname = fname.replace('94', '094')
            aiamap = Map(fname)

        # Create new Map and put EM values in it
        emmap = Map(self.data.copy(), self.meta.copy())
        indices = np.round((tempdata - 4.0) / 0.05).astype(int)
        indices[indices < 0] = 0
        indices[indices > 100] = 100
        #print emmap.shape, indices.shape, tempdata.shape, aiamap.shape, resp.shape
        emmap.data = np.log10(aiamap.data / resp[indices])
        #emmap.data = aiamap.data / resp[indices]

        emmapcubehelix = _cm.cubehelix(s=2.8, r=-0.7, h=1.4, gamma=1.0)
        cm.register_cmap(name='emhelix', data=emmapcubehelix)
        emmap.cmap = cm.get_cmap('emhelix')
    
        return emmap

# 需要导入模块: from sunpy.map import Map [as 别名]
# 或者: from sunpy.map.Map import submap [as 别名]
    def compare(self, display_wlen='171', context_wlen=None, extra_maps=[]):
        #        temp_args=None, temp_kwargs=None,
        #        wlen_args=None, wlen_kwargs=None,
        #        ctxt_args=None, ctxt_kwargs=None,
        #        extr_args=None, extr_kwargs=None):
        valid_wlens = ['94', '131', '171', '195', '211', '335', '304', 'hmi']
        if display_wlen.lower() not in valid_wlens:
            print "Display wavelength provided invalid or None."
            output = self.plot()#*temp_args, **temp_kwargs)
            return output
        save_output = True
        data_dir = self.data_dir
        maps_dir = self.maps_dir
        
        date = self.date
        nmaps = 2 + len(extra_maps)
        if context_wlen:
            nrows = 2
        else:
            nrows = 1
        
        fig = plt.figure(figsize=(24, 14))
        
        fig.add_subplot(nrows, nmaps, nmaps, axisbg='k')
        self.plot()#*temp_args, **temp_kwargs)
        plt.colorbar(orientation='horizontal')
        
        displaymap = Map(data_dir+'{0}/{1:%Y/%m/%d}/aia*{0}*t{1:%H?%M}*lev1?fits'\
            .format(display_wlen, date))
        if isinstance(displaymap, list):
            displaymap = displaymap[0]
        displaymap = aiaprep(displaymap)
        displaymap /= displaymap.exposure_time
        
        fig.add_subplot(nrows, nmaps, 1, axisbg='k')
        displaymap.plot()#*wlen_args, **wlen_kwargs)
        plt.colorbar(orientation='horizontal')
        
        if context_wlen and self.region != None:
            context_plot = fig.add_subplot(nrows, 1, nrows)
            contextmap = Map(data_dir+'{0}/{1:%Y/%m/%d}/aia*{0}*t{1:%H?%M}*lev1?fits'.format(context_wlen, date))
            if isinstance(contextmap, list):
                contextmap = contextmap[0]
            x, y = self.region_coordinate['x'], self.region_coordinate['y']
            contextmap = contextmap.submap([-1000, 1000], [y-300, y+300])
            # Need to figure out how to get 'subimsize' from self. Use the default 150'' for now
            #rect = patches.Rectangle([x-subdx, y-subdx], subimsize[0], subimsize[1], color='white', fill=False)
            rect = patches.Rectangle([x-150, y-150], 300, 300, color='white',
                                     fill=False)
            contextmap.plot()#*ctxt_args, **ctxt_kwargs)
            context_plot.add_artist(rect)
        
        for m, thismap in extra_maps:
            fig.add_subplot(nrows, nmaps, 3+m)
            thismap.plot()#*extr_args, **extr_kwargs)
        
        if save_output:
            error = sys('touch '+maps_dir+'maps/{:%Y/%m/%d/} > shelloutput.txt'.format(date))
            if error != 0:
                sys('{0}{1:%Y}; {0}{1:%Y/%m}; {0}{1:%Y/%m/%d} > shelloutput.txt'\
                        .format('mkdir '+maps_dir+'maps/', date))
            filename = maps_dir+'maps/{:%Y/%m/%d/%Y-%m-%dT%H:%M:%S}_with{}'.\
                    format(date, display_wlen)
            plt.savefig(filename)
            if self.region != None:
                reg_dir = maps_dir + 'maps/region_maps'
                reg_dir = reg_dir + '/{}/'. format(self.region)
                error = sys('touch ' + reg_dir + ' > shelloutput.txt')
                if error != 0:
                    sys('mkdir ' + reg_dir + ' > shelloutput.txt')
                plt.savefig(reg_dir+'{:%Y-%m-%dT%H:%M:%S}'.format(date))
            plt.close()
        else:
            plt.show()

        return

# 需要导入模块: from sunpy.map import Map [as 别名]
# 或者: from sunpy.map.Map import submap [as 别名]
import sunpy.data.sample as s
import numpy as np
from astropy.units import Quantity
import scipy.signal

aia = Map(s.AIA_171_IMAGE)

fig, axs = plt.subplots(1,1)
aia.plot(axes=axs[0])
aia.draw_grid()

pixel_size = 3.4
number_of_pixels = 160

center = np.array([461, 303])
smap = aia.submap(Quantity([center[0] - pixel_size * number_of_pixels/2., center[0] + pixel_size * number_of_pixels/2.], 'arcsec'),
                  Quantity([center[1] - pixel_size * number_of_pixels/2., center[1] + pixel_size * number_of_pixels/2.], 'arcsec'))
smap.plot()

#plate scale for this image is 1.3 arcsec
#plate scale for AIA image is 1.4 arcsec which is very close
psf = np.fromfile('onaxis_psf.dat', dtype=np.float32).reshape(53,53)
psf /= psf.sum()

r = scipy.signal.convolve2d(smap.data, psf, mode='same')

smap.data = r
smap.plot()

new_dim = int(smap.dimensions.x.value * smap.scale.x.value / pixel_size)

rmap = smap.resample(Quantity([new_dim, new_dim], 'pix'), method='spline')