Python Map.data方法代码示例

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

# 需要导入模块: from sunpy.map import Map [as 别名]
# 或者: from sunpy.map.Map import data [as 别名]
    plt.axvline(35.0, color='white')
    plt.axhline(5.6, color='white')
    plt.axhline(7.0, color='white')
    #plt.savefig(path.join(outdir, 'tempsolutions_em={}'.format(np.log10(wid)).replace('.', '_')))
    plt.savefig(path.join(outdir, 'tempsolutions_wid={:.3f}'.format(wid).replace('.', '_')))
    plt.close()

    if n_pars == 3:
        emdata = Map(newmap.emission_measure, newmap.meta.copy())
    else:
        if emwlen == 'three':
            total = np.zeros(newmap.shape)
            for w in ['171', '193', '211']:
                emdata = newmap.calculate_em(w, model=True)
                total += emdata.data
            emdata.data = total/3.0
        elif emwlen == 'all':
            total = np.zeros(newmap.shape)
            for w in ['94', '131', '171', '193', '211', '335']:
                emdata = newmap.calculate_em(w, model=True)
                total += emdata.data
            emdata.data = total/6.0
        else:
            emdata = newmap.calculate_em(emwlen, model=True)
    trueem = np.array(list(heights)*n_temps).reshape(n_temps, n_heights)
    diffem = Map((abs(trueem - emdata.data) / trueem) * 100, newmap.meta.copy())
    #print wid, emdata.xrange, emdata.yrange, emdata.scale
    #print wid, diffem.xrange, diffem.yrange, diffem.scale
    #print wid, fitsmap.xrange, fitsmap.yrange, fitsmap.scale
    fig = plt.figure(figsize=(24, 12))
    ax = fig.add_subplot(1, 3, 1)