ArcPy & NumPy： MemoryError 处理 Landsat 8 图像

Question

我正试图使用ArcPy（经验丰富）和NumPy（新手）比较NDVI计算，但在下面脚本中指示的行遇到了MemoryError。

由RasterToNumPyArray创建的阵列的形状和数据类型是(8191, 8101)和uint16。这是标准陆地卫星8号图像的大小。

这是一个不合理的规模期望处理吗？除法是否促进了数据类型容纳小数（diff和sum成功完成）？如果是，我能以某种方式强制输出整数吗？

import cProfile, numpy
def arcpyNDVI():
    NIRras = arcpy.Raster("LC80460222013184LGN00_B5.TIF")
    REDras = arcpy.Raster("LC80460222013184LGN00_B4.TIF")
    NDVIap = (NIRras - arcpy.sa.Float(REDras))/ (NIRras + REDras)
    NDVIap.save(r'C:/junk/ap_ras.tif')
def numpyNDVI():
    NIRras = arcpy.Raster("LC80460222013184LGN00_B5.TIF")
    NIRll = arcpy.Point(NIRras.extent.XMin,NIRras.extent.YMin)
    NIRcs = NIRras.meanCellWidth
    REDras = arcpy.Raster("LC80460222013184LGN00_B4.TIF") 
    arcpy.env.outputCoordinateSystem = NIRras.spatialReference
    NIRnp = arcpy.RasterToNumPyArray(NIRras)
    REDnp = arcpy.RasterToNumPyArray(REDras)
    diff = NIRnp - REDnp
    sum = NIRnp + REDnp
    NDVInp = diff / sum # MEMORY ERROR HERE
    NDVInpRas = arcpy.NumPyArrayToRaster(NDVInp,NIRll,NIRcs,NIRcs)
    NDVInpRas.save(r'C:/junk/np_ras.tif')
cProfile.runctx('arcpyNDVI()',None,locals())
cProfile.runctx('numpyNDVI()',None,locals())

update1：上面的脚本有时会运行到完成（有时会给出MemoryError），但它输出uint16，而我需要小数。更改线路后：NDVInp = diff / sum到：NDVInp = numpy.true_divide(diff,sum)，我一直得到：

Runtime error 
Traceback (most recent call last):
  File "<string>", line 24, in <module>
  File "C:AnacondaLibcProfile.py", line 49, in runctx
    prof = prof.runctx(statement, globals, locals)
  File "C:AnacondaLibcProfile.py", line 140, in runctx
    exec cmd in globals, locals
  File "<string>", line 1, in <module>
  File "<string>", line 17, in numpyNDVI
MemoryError

Answer 1

答案似乎是始终使用del语句来释放内存，并在最早的时刻删除变量。如果有人想评论这是否是numpy脚本（或Python）的常见/最佳实践，我会洗耳恭听。这是我的工作脚本：

import numpy
NIRras = arcpy.Raster("LC80460222013184LGN00_B5.TIF")
NIRll = arcpy.Point(NIRras.extent.XMin,NIRras.extent.YMin)
NIRcs = NIRras.meanCellWidth
REDras = arcpy.Raster("LC80460222013184LGN00_B4.TIF") 
arcpy.env.outputCoordinateSystem = NIRras.spatialReference
NIRnp = arcpy.RasterToNumPyArray(NIRras)
del NIRras
REDnp = arcpy.RasterToNumPyArray(REDras)
del REDras
diff = NIRnp - REDnp
sum = NIRnp + REDnp
del NIRnp, REDnp
NDVInp = numpy.true_divide(diff,sum) # MEMORY ERROR HERE
NDVInpRas = arcpy.NumPyArrayToRaster(NDVInp,NIRll,NIRcs,NIRcs)
del NDVInp, NIRll, NIRcs
NDVInpRas.save(r'C:/junk/np_ras.tif')
del NDVInpRas