Frage

Mögliches Duplikat:
Verarbeitung großer Satellitenbilder

Ich versuche, die Python-iMAD-Implementierung von Mort Canty auf bitemporalen RapidEye-Multispektralbildern auszuführen.Dabei wird im Grunde die kanonische Korrelation für die beiden Bilder berechnet und dann subtrahiert.Das Problem, das ich habe, ist, dass die Bilder 5000 x 5000 x 5 (Bänder) Pixel haben.Wenn ich dies auf dem gesamten Bild meines Computers ausführe stürzt ab schrecklich und ich muss es ausschalten.

Hat jemand eine Idee, was Python machen kann? Absturz der Computer so?Wenn ich zum Beispiel 2999 x 2999 Pixel pro Band wähle, läuft alles einwandfrei.

8 GB RAM, I7-2617M 1,5 1,5 GHz, Windows7 64 Bit.Ich verwende die 64-Bit-Version von allem:Python(2.7), Numpy, Scipy und Gdal.

Vielen Dank im Voraus!

    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 '---------------------------------'     

Wenn Name == 'hauptsächlich':hauptsächlich()

War es hilfreich?

Lösung

Es hört sich so an, als würde dieser Vorgang einfach mehr Speicher beanspruchen, als Ihr Computer bereitstellen kann.Dies ist eine zu starke Vereinfachung, aber wenn dem System der tatsächlich zu verwendende RAM ausgeht, schreibt es manchmal Speicherabschnitte, die scheinbar weniger genutzt werden, auf Ihre Festplatte, sodass es diesen echten Speicher für etwas anderes verwenden kann.Die Festplatte ist um viele Größenordnungen langsamer als der Hauptspeicher. Wenn Ihre Software also Teile des Speichers benötigt, die auf die Festplatte geschrieben wurden, kann alles sehr langsam werden.Wenn dies in dramatischem Ausmaß geschieht und Teile Ihrer Software und Ihres Betriebssystems ständig versuchen, Teile des Speichers zu nutzen, die ausgelagert (auf die Festplatte geschrieben) wurden, kann Ihre Festplatte bei dem Versuch, hin und her zu suchen, eine große Belastung bekommen. viel schreiben und viel lesen und noch mehr schreiben usw.In einer solchen Situation kann es vorkommen, dass das System nicht mehr reagiert.

Sie könnten sehen, ob das wirklich der Fall ist, indem Sie die Aktivitätsmonitore Ihres Systems beobachten (ich habe vergessen, wie sie unter Windows heißen, aber ich weiß, dass sie da sind;eine Software, die Ihnen zeigt, wie viel Speicher zugewiesen ist, verwendet wird usw. und ein schönes Diagramm für Sie erstellt).Während Sie sich diese ansehen, starten Sie Ihr Programm und beobachten Sie, wie die Speicherbelegungsrate aussieht.

Es gibt wahrscheinlich einige Möglichkeiten, die Speichernutzung in diesem Code zu verringern, wenn weniger Dinge gleichzeitig im Speicher gehalten werden, aber ich sehe nicht ohne weiteres, um welche es sich dabei handelt.Sie könnten Ihrem System auch mehr RAM hinzufügen, in der Hoffnung, das Problem dadurch zu beheben.

Lizenziert unter: CC-BY-SA mit Zuschreibung
Nicht verbunden mit StackOverflow
scroll top