Расчеты на большом спутниковом снимке вылетает питон и зависает компьютер
Возможный дубликат:
Большая спутниковая обработка изображений
Я пытаюсь запустить реализацию Mort Canty's Python iMAD для битемпоральных изображений RapidEye Multispectral. Который в основном вычисляет каноническую корреляцию для двух изображений и затем вычитает их. У меня проблема в том, что изображения имеют 5000 х 5000 х 5 (полос) пикселей. Когда я запускаю это на всем изображении, мой компьютер ужасно падает, и мне приходится его выключать.
Кто-нибудь есть идеи, что может заставить Python сбить компьютер, как это? Если я, например, выберу 2999x2999 пикселей на полосу, все будет нормально.
Оперативная память 8 ГБ, I7-2617M 1,5 1,5 ГГц, 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 '---------------------------------'
if name == 'main': main ()
1 ответ
Похоже, что эта операция просто занимает больше памяти, чем ваш компьютер может предоставить. Это упрощение, но когда системе не хватает реальной оперативной памяти для использования, она иногда записывает разделы памяти, которые, по-видимому, используются меньше на вашем жестком диске, поэтому она может использовать эту реальную память для чего-то другого. Жесткий диск на много порядков медленнее основной памяти, поэтому, когда вашему программному обеспечению требуются части памяти, которые были записаны на диск, все может стать очень медленным. Когда это происходит в драматическом масштабе, и части вашего программного обеспечения и операционной системы постоянно пытаются использовать фрагменты памяти, которые были выгружены (записаны на диск), ваш жесткий диск может получить большую тренировку, пытаясь искать туда-сюда, написание большого количества материала и чтение большого количества материала, написание большего количества материала и т. д. Система может перестать отвечать на запросы в такой ситуации.
Вы можете увидеть, действительно ли это происходит, наблюдая за мониторами активности вашей системы (я забываю, как они называются в Windows, но я знаю, что они есть; некоторые программы показывают, сколько памяти выделено, используется и т. Д., и рисует хороший график для вас). Наблюдая за ними, запустите вашу программу и посмотрите, как выглядит скорость выделения памяти.
Вероятно, есть несколько способов облегчить использование памяти в этом коде, если в памяти одновременно хранится меньше вещей, но я не вижу, как они выглядят. Вы также можете добавить больше оперативной памяти в вашу систему в надежде, что это исправит это.