NumPy memory error on large rasters

Question

I have multiple large rasters. Their 5 percentile and 95 percentile are needed to be calculated. After loading the rasters to the ArcMap, I am using the following codes -

import numpy
import arcpy

rasters = ["raster1.tif", "raster2.tif", "raster3.tif", "raster4.tif"]

for i in rasters:
    array = arcpy.RasterToNumPyArray(i,nodata_to_value=10000000000)
    marray = numpy.ma.masked_values(array,10000000000)
    del array
    percentile5 = numpy.percentile(marray.compressed(),(5))
    percentile95 = numpy.percentile(marray.compressed(),(95))
    del marray
    print percentile5 , percentile95 

For the first two rasters, it provides the percentile values. For the third raster, it is showing memory errors as provided below.

Traceback (most recent call last):
  File "E:\test\Stat\Inlinear data.py", line 89, in <module>
    marray = numpy.ma.masked_values(array,10000000000)
  File "C:\Python27\ArcGIS10.5\lib\site-packages\numpy\ma\core.py", line 2200, in masked_values
    condition = umath.less_equal(mabs(xnew - value), atol + rtol * mabs(value))
MemoryError

Failed to execute (IntelinAir).

I am using ArcMap 10.5 in desktop with 32 GB Ram.

Is there any way to get rid of this memory issue?

Answer 1

A MemoryError means that you have exhausted the memory available. If your ArcGIS is 32-bit, then that limit is 4GB. You will see this error with a large array with float64:

#  Python 2.7.13 (v2.7.13:a06454b1afa1, Dec 17 2016, 20:42:59) [MSC v.1500 32 bit (Intel)] on win32
import numpy as np
ar = np.empty((16963, 11367), np.float32)
del ar
ar = np.empty((16963, 11367), np.float64)
# MemoryError

As mentioned by others, the only solution for Esri tools is to install 64-bit Background Geoprocessing, which is a separate installation on top of ArcGIS Desktop. After it's installed, you may find (e.g.) c:\Python27\ArcGISx6410.5\python.exe to process with the unmodified script in the question.

---

An open source solution is to use 64-bit Python with other tools, like GDAL or rasterio, which can use the full amount of RAM on your computer.

# Python 2.7.15 | packaged by conda-forge | (default, Jul 27 2018, 15:07:58) [MSC v.1500 64 bit (AMD64)]
import numpy as np
import rasterio

rasters = ["raster1.tif", "raster2.tif", "raster3.tif", "raster4.tif"]

for rname in rasters:
    with rasterio.open(rname, 'r') as src:
        mar = src.read(masked=True)
    percentiles = np.percentile(mar.compressed(), [5, 95])
    print('{0}: {1}'.format(rname, list(percentiles)))

A slightly longer recipie can be done with osgeo.gdal too.

Answer 2

Slightly off topic, but anyway.

Brute force is not always the best approach, because there will be one that won't fit into memory. This will do the job for any size raster:

import arcpy, os, traceback, sys
from arcpy import env
env.overwriteOutput = True
from arcpy.sa import *

def getCells(aRaster):
    arcpy.MakeRasterLayer_management(aRaster, "RASTERLAYER")
    tbl=arcpy.da.TableToNumPyArray("RASTERLAYER",("value","count"))
    return tbl[0][1]

inRaster=arcpy.Raster("DEM_2006")
minZ = float(arcpy.GetRasterProperties_management(inRaster, "MINIMUM")[0])
maxZ = float(arcpy.GetRasterProperties_management(inRaster, "MAXIMUM")[0])
area=Con(~IsNull(inRaster),1)
target=getCells(area)/20
while True:
    mid=(maxZ+minZ)/2
    if (maxZ-minZ)<0.01:break
    below=Con(inRaster<mid,1)
    below=getCells(below)
    arcpy.AddMessage("Current value %s" % mid)
    if below>target:maxZ=mid
    else: minZ=mid
arcpy.AddMessage('5 percentile value =%s' %mid)