큰 위성 이미지에 대한 계산이 Python과 충돌하고 컴퓨터가 정지됩니다.
-
11-12-2019 - |
문제
가능한 중복:
대형위성영상처리
저는 양측 RapidEye Multispectral 이미지에서 Mort Canty의 Python iMAD 구현을 실행하려고 합니다.기본적으로 두 이미지의 표준 상관관계를 계산한 다음 이를 뺍니다.내가 겪고 있는 문제는 이미지가 5000 x 5000 x 5(밴드) 픽셀이라는 것입니다.내 컴퓨터 전체 이미지에서 이것을 실행하면 충돌 정말 꺼야 해요.
파이썬으로 무엇을 만들 수 있는지 아는 사람 있나요? 충돌 컴퓨터가 이렇게 돼?예를 들어 밴드당 2999x2999픽셀을 선택하면 모든 것이 잘 실행됩니다.
8GB RAM, I7-2617M 1.5 1.5GHz, 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 '---------------------------------'
만약에 이름 == '기본':기본()
해결책
이 작업은 컴퓨터가 제공할 수 있는 것보다 더 많은 메모리를 사용하는 것 같습니다.이는 지나치게 단순화한 것이지만, 시스템에서 사용할 실제 RAM이 부족한 경우에는 덜 사용되는 것처럼 보이는 메모리 섹션을 하드 디스크에 기록하여 해당 실제 메모리를 다른 용도로 사용할 수 있습니다.하드 디스크는 주 메모리보다 훨씬 느리므로 소프트웨어에 디스크에 기록된 메모리 일부가 필요한 경우 모든 것이 매우 느려질 수 있습니다.이러한 일이 극적인 규모로 발생하고 소프트웨어 및 운영 체제의 일부가 교체된(디스크에 기록된) 메모리 조각을 지속적으로 사용하려고 하면 하드 드라이브가 앞뒤로 탐색하려고 큰 작업을 수행할 수 있습니다. 많은 것을 쓰고, 많은 것을 읽고, 더 많은 것을 쓰세요.이와 같은 상황에서는 시스템이 응답하지 않을 수 있습니다.
시스템의 활동 모니터를 보면 이것이 실제로 무슨 일이 일어나고 있는지 확인할 수 있습니다(Windows에서 뭐라고 부르는지는 잊어버렸지만 거기에 있다는 것은 알고 있습니다.할당된 메모리 양, 사용 중인 메모리 양 등을 보여주고 멋진 그래프를 그리는 일부 소프트웨어).이를 보면서 프로그램을 시작하고 메모리 할당 비율이 어떤지 살펴보세요.
한 번에 메모리에 보관되는 항목의 수가 적다면 이 코드에서 메모리 사용량을 줄일 수 있는 몇 가지 방법이 있을 수 있지만 그 내용을 직접 확인할 수는 없습니다.이 문제가 해결되기를 바라면서 시스템에 RAM을 더 추가할 수도 있습니다.