1

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?

  • 1
    32 GiB of ram is meaningless to ArcGIS, it's still a 32bit process and thereby subject to 4 GiB maximum RAM allocation. What is the memory issue? if there's an error message I'd like to see it. If it takes 3 rasters before it fails try running them one at a time as a process and let Windows settle down between each one. – Michael Stimson Aug 7 '18 at 4:56
  • @MichaelStimson I have added the error message. – Sourav Aug 7 '18 at 5:19
  • 1
    What are the dimensions (bands, rows, cols) and datatype of the rasters? – user2856 Aug 7 '18 at 5:26
  • What happens when you run each image in one process? Does it still give the memory error on the same raster? i.e. is it the iteration or is it the data that's likely at fault here. – Michael Stimson Aug 7 '18 at 5:30
  • 1
    Please don't add additional info (i.e. the dimensions (bands, rows, cols) and datatype) as a comment, edit your question. – user2856 Aug 7 '18 at 5:57
4

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.

  • 1
    Whilst generally I agree with you, installing 64bit GDAL would solve this problem it may be too big of an idiom shift for the OP. They may be on a course to learn arcpy or in a workplace that only allows Esri products. As there is no version (yet) of ArcGIS Desktop that is 64bit a viable Esri alternative is ArcGIS Pro 64bit which offers the same function pro.arcgis.com/en/pro-app/arcpy/functions/… and can access 32 GiB of RAM. – Michael Stimson Aug 7 '18 at 6:11
  • While ArcMap will never be 64-bit, there has been a 64-bit arcpy option since at least 10.2 (probably 10.1), via the 64-bit Background Geoprocessing package. This is the option with the least learning curve. It requires an install, but no further licensing. – Vince Aug 7 '18 at 11:09
  • 1
    Thanks for the answer. I just ran the script in ArcGIS Pro which comes with 64-bit python 3.5. Ran the whole program without any issue! Python 3 dealt with the matrix excellently. – Sourav Aug 8 '18 at 1:25
1

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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.