I have a large .tif raster with integer values. I wish to generate Jenks classification breaks using code and then use those breaks to reclassify my .tif file.

I came across a block of code that I am using to define my Jenks function. Its source is: https://gist.github.com/drewda/1299198

So far I have imported my raster, converted it to a (233,679) 2D array and then converts to a (158207,) 1D array.

Then I attempt to run the function on the array but get the error:

Traceback (most recent call last):
  File "C:\Users\agates44\Desktop\FinalProject\script1.py", line 80, in <module>
    new = getJenksBreaks(oneDarr,9)
  File "C:\Users\agates44\Desktop\FinalProject\script1.py", line 40, in getJenksBreaks

My code is below

# Define the folder path
folderPath = r'C:\Users\agates44\Desktop\FinalProject' 

# Import modules
import arcpy 
import numpy as np

# Set workspace
from arcpy import env
env.workspace = folderPath
# Allow overwrite
env.overwriteOutput = True

# Get Input Raster properties
inRas = arcpy.Raster('smallorange.tif')
lowerLeft = arcpy.Point(inRas.extent.XMin,inRas.extent.YMin)
cellSize = inRas.meanCellWidth

# Convert Raster to numpy array
arr = arcpy.RasterToNumPyArray(inRas,nodata_to_value=0)
print arr.shape
print arr.dtype
oneDarr = arr.ravel()
print oneDarr.shape
print oneDarr.dtype

def getJenksBreaks(dataList, numClass):
  mat1 = []
  for i in range(0,len(dataList)+1):
    temp = []
    for j in range(0,numClass+1):
  mat2 = []
  for i in range(0,len(dataList)+1):
    temp = []
    for j in range(0,numClass+1):
  for i in range(1,numClass+1):
    mat1[1][i] = 1
    mat2[1][i] = 0
    for j in range(2,len(dataList)+1):
      mat2[j][i] = float('inf')
  v = 0.0
  for l in range(2,len(dataList)+1):
    s1 = 0.0
    s2 = 0.0
    w = 0.0
    for m in range(1,l+1):
      i3 = l - m + 1
      val = float(dataList[i3-1])
      s2 += val * val
      s1 += val
      w += 1
      v = s2 - (s1 * s1) / w
      i4 = i3 - 1
      if i4 != 0:
        for j in range(2,numClass+1):
          if mat2[l][j] >= (v + mat2[i4][j - 1]):
            mat1[l][j] = i3
            mat2[l][j] = v + mat2[i4][j - 1]
    mat1[l][1] = 1
    mat2[l][1] = v
  k = len(dataList)
  kclass = []
  for i in range(0,numClass+1):
  kclass[numClass] = float(dataList[len(dataList) - 1])
  countNum = numClass
  while countNum >= 2:#print "rank = " + str(mat1[k][countNum])
    id = int((mat1[k][countNum]) - 2)
    #print "val = " + str(dataList[id])
    kclass[countNum - 1] = dataList[id]
    k = int((mat1[k][countNum] - 1))
    countNum -= 1
  return kclass

new = getJenksBreaks(oneDarr,9)
print new

### used after running getJenksBreaks()
##def classify(value, breaks):
##  for i in range(1, len(breaks)):
##    if value < breaks[i]:
##      return i
##  return len(breaks) - 1

How can I figure out which line of my code is using all the memory and what is a workaround I can use?

****** EDIT The function ran (although it took forever) on the raster size listed above. But that was only a subset of my much larger dataset.

(2400, 2337) 2D array (5608800,) 1D array

This is when I get the memory error.

  • 1
    What does task manager say about your memory usage? 233 x 679 isn't a large raster, are you sure that your raster is only 233 rows and 679 columns? Be aware that ArcGIS Desktop is only a 32bit application and can only address 2**32 bytes of memory for each instance of Arc - program, script and data. May 29, 2018 at 2:51
  • @MichaelStimson I apologize, I was wrong. The raster that is generating the memory error is a much larger raster. I was able to generate break values for the smaller subset of the array, which was my original post: [0, 35, 150, 350, 516, 696, 896, 3337, 4132, 428760.0]. Yet I need to run the function on the larger raster, which I fixed in the edit below. May 29, 2018 at 3:28
  • @MichaelStimson My task manager shows my memory reaching 87% capacity when I run it on the larger raster, and then I get the MemoryError. May 29, 2018 at 3:35
  • 87% isn't a valuable figure, what is important is the physical load in mebibytes, Esri can only address 4 GiB of memory and it sounds like you're filling that up. Considering you're reading the whole raster into memory and performing a purely python function I would suggest you consider either ArcGIS Pro (64 bit) or GDAL 64 bit.. my preference is for GDAL. Read gis.stackexchange.com/questions/32995/… about using GDAL to read a raster into a numpy array. May 29, 2018 at 3:45
  • How long did the 233, 679 (158207) array take to run? I roughly guesstimated ~5 hours based on how long it took to run on a 10,000 element array.
    – user2856
    May 29, 2018 at 5:16

1 Answer 1


A 2400*2337 (assuming single band 32bit) raster is only ~20mb. I'm not sure why you get a MemoryError. I can create a dummy 32bit 2400*2337 ndarray in numpy and run getJenksBreaks on it "fine"* in 32bit Python 2.7 from ArcGIS Desktop 10.3.

*By "fine", I mean, it will run, but very slowly and I don't see the memory use climb, but I only let it run for a few minutes.

However... You will not be able to run the pure python Jenks classification on your entire raster. You will need to sample your raster

Jenks classification in pure python is horribly slow. You should use a native (i.e. C/C++) jenks library, such as jenkspy.

And you need to sample your data, don't try to run it on the entire raster. You can use numpy.random.choice(data, size=sample_size, replace=False) to sample it. I have tested and found fairly similar class breaks with 25,000 - 100,000 samples (from a 32 bit flow accumulation raster).

Consider the following plots of the time to calculate Jenks breaks against the number of elements. It was generated (in Excel, one day I'll learn how to use matplotlib...) from the results of the code that follows using 100 - 500,000 samples of the dummy data The solid line shows the pure python Jenks calculation and the dotted line is using jenkspy. Even the jenkspy calculation is quite slow once you get over 10,000 values. I haven't tried running it on the whole 5,608,800, I'm sure it would take a very, very, very long time...

enter image description here enter image description here

import time

import jenkspy
import numpy as np

def getJenksBreaks(dataList, numClass):
    blah blah blah, same code as yours...

nbands, nrows, ncols = 1, 2400, 2337
zmin, zmax = 0, 12345

data = np.random.randint(zmin, zmax, (nbands, nrows, ncols)).astype(np.float32).ravel()

samples = (100, 500, 1000, 5000, 10000, 50000, 100000, 500000)
for sample in samples:
    sampled = np.random.choice(data, size=sample, replace=False)

    breaks = jenkspy.jenks_breaks(sampled, nb_class=5)
    tj = time.time() - t

    if sample <= 100000:  # Don't bother trying more than 100K
        t=time.time()    # pixels with the pure python version
        breaks = getJenksBreaks(sampled, 5)
        print('%s\t%s\t%s'%(sample, tj, time.time() - t))
        print('%s\t%s'%(sample, tj)
  • it is a single band 32bit flow accumulation raster, 21.40mb. Running on ArcMap Desktop 10.5 May 29, 2018 at 19:49
  • @AustinGates Doesn't matter. You will not be able to run the pure python Jenks code on your entire large raster. Just sample it.
    – user2856
    May 29, 2018 at 23:52

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.