1

可能的重复:
大卫星图像处理

我正在尝试在双时相 RapidEye 多光谱图像上运行 Mort Canty 的 Python iMAD 实现。它基本上计算了两个图像的典型相关,然后减去它们。我遇到的问题是图像的大小为 5000 x 5000 x 5(条带)像素。当我在整个图像上运行此程序时,我的计算机严重崩溃,我必须将其关闭。

有谁知道什么可以让 python像这样使计算机崩溃?例如,如果我选择每个波段 2999x2999 像素,一切运行正常。

8 GB 内存,I7-2617M 1.5 1.5 ghz,Windows7 64 位。我使用 64 位版本的一切:python(2.7)、numpy、scipy 和 gdal。

先感谢您!

    def covw(dm,w):
    # weighted covariance matrix and means 
    # from (transposed) data array    
       N = size(dm,0) 
       n = size(w)
       sumw = sum(w)
       ws = tile(w,(N,1))
       means = mat(sum(ws*dm,1)/sumw).T
       means = tile(means,(1,n))
       dmc = dm - means
       dmc = multiply(dmc,sqrt(ws))
       covmat = dmc*dmc.T/sumw
       return (covmat,means)

    def main():
    # ------------test---------------------------------------------------------------    
    if len(sys.argv) == 1:        
    (sys.argv).extend(['-p','[0,1,2,3,4]','-  
    d','[0,4999,0,4999]',
'c://users//pythonxy//workspace//1uno.tif','c://users//pythonxy//workspace//2dos.tif'])
    # -------------------------------------------------------------------------------        

options, args = getopt.getopt(sys.argv[1:],'hp:d:')
pos = None
dims = None            
for option, value in options:
    if option == '-h':
        print 'Usage: python %s [-p "bandPositions" -d "spatialDimensions"] 
        filename1   filename2' %sys.argv[0]
        print '       bandPositions and spatialDimensions are quoted lists, 
        e.g., -p "[0,1,3]" -d "[0,400,0,400]"  \n'
        sys.exit(1) 
    elif option == '-p':
        pos = eval(value)
    elif option == '-d':
        dims = eval(value) 
if len(args) != 2:
    print 'Incorrect number of arguments'
    print 'Usage: python %s [-p "bandspositions" -d "spatialdimensions"] 
    filename1 filename2 \n' %sys.argv[0]
    sys.exit(1)                                    
gdal.AllRegister()
fn1 = args[0]
fn2 = args[1]
path = os.path.dirname(fn1)
basename1 = os.path.basename(fn1)
root1, ext = os.path.splitext(basename1)
basename2 = os.path.basename(fn2)
outfn = path+'\\MAD['+basename1+'-'+basename2+']'+ext
inDataset1 = gdal.Open(fn1,GA_ReadOnly)     
inDataset2 = gdal.Open(fn2,GA_ReadOnly)
cols = inDataset1.RasterXSize
rows = inDataset1.RasterYSize    
bands = inDataset1.RasterCount
cols2 = inDataset2.RasterXSize
rows2 = inDataset2.RasterYSize    
bands2 = inDataset2.RasterCount
if (rows != rows2) or (cols != cols2) or (bands != bands2):
    sys.stderr.write("Size mismatch")
    sys.exit(1)
if pos is None:
    pos = range(bands)
else:
    bands = len(pos) 
if dims is None:
    x0 = 0
    y0 = 0
else:
    x0 = dims[0]
    y0 = dims[2]  
    cols = dims[1]-dims[0] + 1  
    rows = dims[3]-dims[2] + 1                       
# initial weights
wt = ones(cols*rows)      
# data array (transposed so observations are columns)
dm = zeros((2*bands,cols*rows),dtype='float32')
k = 0
for b in pos:
    band1 = inDataset1.GetRasterBand(b+1)
    band1 = band1.ReadAsArray(x0,y0,cols,rows).astype(float)
    dm[k,:] = ravel(band1)
    band2 = inDataset2.GetRasterBand(b+1)
    band2 = band2.ReadAsArray(x0,y0,cols,rows).astype(float)        
    dm[bands+k,:] = ravel(band2)
    k += 1
print '========================='
print '       iMAD'
print '========================='
print 'time1: '+fn1
print 'time2: '+fn2   
print 'Delta    [canonical correlations]'   
# iteration of MAD        
delta = 1.0
oldrho = zeros(bands)
iter = 0
while (delta > 0.001) and (iter < 50):    
#     weighted covariance matrices and means 
    sigma,means = covw(dm,wt)          
    s11 = mat(sigma[0:bands,0:bands])
    s22 = mat(sigma[bands:,bands:]) 
    s12 = mat(sigma[0:bands,bands:])
    s21 = mat(sigma[bands:,0:bands])
#     solution of generalized eigenproblems
    s22i = mat(linalg.inv(s22))
    lama,a = linalg.eig(s12*s22i*s21,s11) 
    s11i = mat(linalg.inv(s11))    
    lamb,b = linalg.eig(s21*s11i*s12,s22) 
#     sort a   
    idx = argsort(lama)
    a = a[:,idx]
#     sort b         
    idx = argsort(lamb)
    b = b[:,idx]           
#     canonical correlations        
    rho = sqrt(real(lamb[idx]))             
#     normalize dispersions   
    a = mat(a)
    tmp1 = a.T*s11*a
    tmp2 = 1./sqrt(diag(tmp1))
    tmp3 = tile(tmp2,(bands,1))
    a = multiply(a,tmp3)
    b = mat(b) 
    tmp1 = b.T*s22*b
    tmp2 = 1./sqrt(diag(tmp1))
    tmp3 = tile(tmp2,(bands,1))
    b = multiply(b,tmp3)
#     assure positive correlation
    tmp = diag(a.T*s12*b)
    b = b*diag(tmp/abs(tmp))
#     canonical and MAD variates
    U = a.T*mat(dm[0:bands,:]-means[0:bands,:])    
    V = b.T*mat(dm[bands:,:]-means[bands:,:])           
    MAD = U-V  
#     new weights        
    var_mad = tile(mat(2*(1-rho)).T,(1,rows*cols))    
    chisqr = sum(multiply(MAD,MAD)/var_mad,0)
    wt = 1-stats.chi2.cdf(chisqr,[bands])
#     continue iteration         
    delta = sum(abs(rho-oldrho))
    oldrho = rho
    print delta
    iter += 1   
# write results to disk
driver = inDataset1.GetDriver()    
outDataset = driver.Create(outfn,cols,rows,bands+1,GDT_Float32)
projection = inDataset1.GetProjection()
geotransform = inDataset1.GetGeoTransform()
if geotransform is not None:
    gt = list(geotransform)
    gt[0] = gt[0] + x0*gt[1]
    gt[3] = gt[3] + y0*gt[5]
    outDataset.SetGeoTransform(tuple(gt))
if projection is not None:
    outDataset.SetProjection(projection)        
for k in range(bands):        
    outBand = outDataset.GetRasterBand(k+1)
    outBand.WriteArray(resize(MAD[k,:],(rows,cols)),0,0) 
    outBand.FlushCache()
outBand = outDataset.GetRasterBand(bands+1)    
outBand.WriteArray(resize(chisqr,(rows,cols)),0,0) 
outBand.FlushCache()    
outDataset = None
inDataset1 = None
inDataset2 = None  
print 'result written to: '+outfn
print '---------------------------------'     

如果名称== '主要': main()

4

1 回答 1

1

听起来这个操作占用的内存比您的计算机所能提供的更多。这是一种过度简化,但是当系统用完实际可用的 RAM 时,它有时会将似乎使用较少的内存部分写入硬盘,因此它可以将实际内存用于其他用途。硬盘比主存慢很多数量级,所以当您的软件需要已写入磁盘的部分内存时,一切都会变得非常慢。当这种情况以惊人的规模发生时,并且您的软件和操作系统的某些部分不断尝试使用已被换出(写入磁盘)的内存,您的硬盘驱动器可能会在尝试来回寻找时进行重大锻炼,写很多东西,读很多东西,写更多东西,等等。

您可以通过观察系统的活动监视器来了解这是否真的发生了(我忘记了它们在 Windows 上的名称,但我知道它们在那里;一些软件可以显示分配的内存量、正在使用的内存量等,并为您绘制一个漂亮的图表)。在观看这些的同时,启动您的程序并观察内存分配率的样子。

可能有一些方法可以减轻这段代码中的内存使用,如果一次在内存中保存的东西更少,但我不知道它们是什么。您还可以向系统添加更多 RAM,希望它能解决此问题。

于 2012-05-18T19:56:12.590 回答