reclassify using Numpy, out of memory

18836
37
Jump to solution
05-19-2016 03:41 PM
RebeccaStrauch__GISP
MVP Emeritus

EDIT: although I marked Dan's answer re: making the array an integer array, which was the last piece I needed, I created a summary at the end https://community.esri.com/message/610325#comment-610325  to put all the pieces that worked together.  I suggest starting there if running into this problem.

Bottom line: trying to get reclass elevation raster (in m), that include “water” and NoData, and Denali foothills, to 3-4 class (based on feet---keeping the converted layer, not  in memory), then convert to polygons.

For the study are I used for development (and one other), the raster Reclassify and arcpy.RasterToPolygon_conversion this worked relatively fast. (max few minutes)

I’m now running it on an area about 2x larger (~25m cells, 6657 col 13035 row), with a wider range of input values.  When using my script/tool or testing in the python window, I run out of memory after about 4 hours.

So I’m diving into Numpy and https://community.esri.com/people/Dan_Patterson ‘s blog https://community.esri.com/blogs/dan_patterson/2016/03/12/reclassify-raster-data-simply

After a clean reboot, I’m trying my test script in the python window, but I’m still getting an out of error message.  I’m not sure if I have something in error in my script, or it is just too large for the 32bit-python to handle.  So, I’m open to suggestions…..maybe processing half at a time??  I have tried but a full raster with no nodata, and one that has nodata outside the side area boundary.

Test script I am using (parts borrowed from Dan's demo scripts)

import time
import os
import arcpy
from arcpy import env
from arcpy.sa import *
import numpy as np
from textwrap import dedent

arcpy.CheckOutExtension("spatial")

def arr_reclass(a, bins=[], new_bins=[], mask=False, mask_val=None):
  """a    - integer or floating point array to be reclassed using
     bins - sequential list/array of the lower limits of each class
        include one value higher to cover the upper range.
     mask - whether the raster contains nodata values or values to
        be masked with mask_val
     array dimensions will be squeezed 
  """ 
  a_rc = np.zeros_like(a)
  if (len(bins) < 2): # or (len(new_bins <2)):
    print("Bins = {} new = {} won't work".format(bins,new_bins))
    return a
  if len(new_bins) < 2:
    new_bins = np.arange(1,len(bins)+2)   
  new_classes = zip(bins[:-1],bins[1:],new_bins)
  for rc in new_classes:
    q1 = (a >= rc[0])
    q2 = (a < rc[1])
    z = np.where(q1 & q2, rc[2],0)
    a_rc = a_rc + z
  return a_rc

theWorkspace = r'C:\___bear2016\_Skwentna2017\Skwentna_prep.gdb'
arcpy.env.workspace = theWorkspace
arcpy.env.overwriteOutput = True
wsPath, wsGDB = os.path.split(theWorkspace)
print("Project gdb (vector, i.e. for points, lines, polygons): {0}".format(theWorkspace))

inStudy = r'C:\___bear2016\_Skwentna2017\Skwentna_prep.gdb\Study16b'
stdyPath, stdyName = os.path.split(inStudy)
print("Study boundary: {0}".format(stdyName))
spatialref = arcpy.Describe(inStudy).spatialReference

rasterIn = r'C:\___bear2016\_Skwentna2017\Skwentna_raster.gdb\elevExtent'
#rasterIn = r'C:\___bear2016\_Skwentna2017\Skwentna_raster.gdb\elevClip'
inRas = arcpy.Raster(rasterIn)

lowerLeft = arcpy.Point(inRas.extent.XMin, inRas.extent.YMin)
cellSize = inRas.meanCellWidth

rasterGDB, rasterFile = os.path.split(rasterIn)
elevUnits = "meters"
elevCutoffFt = "6000"
flatElevBreak = "2000"
elevBreaks = flatElevBreak.split("; ")
elevBreaks.sort(key=float)  
flatElevBreak = ";".join(elevBreaks)
flatElevBreak
'2000'
hydroBuff = r'C:\___bear2016\_Skwentna2017\Skwentna_prep.gdb\hydroBuff250'
finalTransZones = arcpy.os.path.join(theWorkspace, "TransAreas")
outRiparian = arcpy.os.path.join(theWorkspace, "TransRiparian")
pElevClass = arcpy.os.path.join(theWorkspace, "elevReclassPoly")  
outDEMtmp = arcpy.os.path.join(theWorkspace, "elevReclass")  

if elevUnits == "meters":   #  3.280839895013123
  toElevFeet = 3.28084
else:      # note, IF had to go from feet to meters, factor would be 0.3048
  toElevFeet = 1
elevFeet = (inRas * toElevFeet)
saveElevFeet = arcpy.os.path.join(rasterGDB, "elevFt")
elevFeet.save(saveElevFeet)
elevMin = elevFeet.minimum  - 1
elevMax = elevFeet.maximum  + 1
noDataReclass = 99999

elevSplit = flatElevBreak.split(";")
lstCnt = 1
#reclassArg = []
bins = []
new_bins = []
for reclassValue in elevSplit:
  if lstCnt == 1:
    #reclassArg.append([int(elevMin), int(reclassValue), int(reclassValue)])
    bins.append(int(elevMin))
    bins.append(int(-1))
    bins.append(int(reclassValue))
    new_bins.append(int(-1))
    new_bins.append(int(reclassValue))
  else:
    #reclassArg.append([int(prevValue), int(reclassValue), int(reclassValue)])
    bins.append(int(reclassValue))    
    new_bins.append(int(0))    
    new_bins.append(int(reclassValue))    
  lstCnt += 1
  prevValue = reclassValue
#reclassArg.append([int(prevValue), int(elevCutoffFt), int(elevCutoffFt)])
#reclassArg.append([int(elevCutoffFt), int(elevMax), 9])

bins.append(int(elevCutoffFt))
bins.append(int(elevMax))
new_bins.append(int(elevCutoffFt))
new_bins.append(9)  
new_classes = zip(bins[:-1],bins[1:],new_bins)

print("   ...Old Ranges:\n   {0}".format(bins))
print("   ...New values:\n   {0}".format(new_bins))
print("   ...New classes:\n   {0}".format(new_classes))


arr = arcpy.RasterToNumPyArray(inRas, nodata_to_value=9999999)
mask = False
mask_val = None
print(arr)
arr_rc = arr_reclass(arr, bins=bins, new_bins=new_bins)

The error I am getting is on the last line, however, looking at the resources, there is plenty of RAM available, so I'm assuming it is just the Arcmap/Python limit.

Output from script (in ArcMap python window)

Project gdb (vector, i.e. for points, lines, polygons): C:\___bear2016\_Skwentna2017\Skwentna_prep.gdb

Study boundary: Study16b

...Old Ranges:

   [-107504, -1, 2000, 6000, 20036]

   ...New values:

   [-1, 2000, 6000, 9]

   ...New classes:

   [(-107504, -1, -1), (-1, 2000, 2000), (2000, 6000, 6000), (6000, 20036, 9)]

[[  205.16717529   205.25009155   205.29598999 ...,  1037.9387207

   1037.86242676  1037.83959961]

[  205.22860718   205.3127594    205.37818909 ...,  1037.93701172

   1037.84533691  1037.76965332]

[  205.29447937   205.38237      205.46986389 ...,  1037.91918945

   1037.75964355  1037.76477051]

...,

[ 1481.97961426  1492.3326416   1504.91308594 ...,   510.62741089

    504.97958374   501.19992065]

[ 1469.56079102  1486.19445801  1505.37487793 ...,   506.6071167

    501.30090332   497.77194214]

[ 1453.34509277  1471.94152832  1493.07250977 ...,   501.66564941

    497.09884644   494.20864868]]

Runtime error

Traceback (most recent call last):

  File "<string>", line 115, in <module>

  File "<string>", line 30, in arr_reclass

MemoryError

>>>

If I need to split it due to size, some guidance for that would be good.  The one nice thing about this numpy process, I get the error immediately instead of 4-hours in....

Thanks for any help!

ps - Dan....is this question seen by more than the 13 members/followers?  Nice "place" but want to make sure others have a chance to find/reply

tagging a few more areas  Pythonpython snippetsImagery and Remote Sensing

37 Replies
TimothyHales
Esri Notable Contributor

The majority represents those not signed into the community. With the Re: Double log in to geonet​ issue, many viewers are not be authenticated correctly. There is also an Jive bug with views being double counted. Also views are counted each time you come to this thread, read it in an activity stream, or interact within your inbox. With the 30+ comments of Rebecca Strauch, GISP​ and Dan Patterson​ I'm sure there are 100+ views between the two of you.

0 Kudos
RebeccaStrauch__GISP
MVP Emeritus

I probably have 100+ myself, between double logins, multiple machines, and the iPad wanting to refresh the tab many times after I've changed focus.   I don't read too much into the stats, but is always interesting to see numbers...even if only relative and not completely accurate.

0 Kudos
DanPatterson_Retired
MVP Emeritus

I just thought the programmers were trolling for tips

0 Kudos
DanPatterson_Retired
MVP Emeritus

NOW!  lets try something completely different which may simplify a few steps, like the float to in stuff...

I had this one buried

>>> # and now for something completely different....
>>> r = 100
>>> c = 120
>>> z = np.abs(100 * np.random.random_sample((r,c)) -100)  # make an array with rand 0-100
>>> np.save("f:/test/Rebecca/z_7_14k.npy",z)
>>> # ----- compare ---- select
>>> condlist = [z<20, z<40, z<60, z<80, z>=80]
>>> choicelist = [20,40,60,80,100]
>>> z_rc = np.select(condlist, choicelist)
>>>
>>> np.set_printoptions(edgeitems=5,linewidth=100,precision=1,suppress=True,threshold=200)
>>>
>>> z
array([[ 71.8,  39.9,  58.4,  8.1,  95.6, ...,  49. ,  56.7,  16. ,  64.8,  16.4],
      [ 92. ,  44.8,  5.2,  22.8,  83.2, ...,  57.4,  81. ,  39.5,  99.7,  45.9],
      [ 68.2,  43.3,  72. ,  57.7,  30.6, ...,  80.7,  82.6,  30.5,  89.6,  7. ],
      [ 44.6,  86.8,  59.7,  1.4,  83.1, ...,  41.8,  54.8,  64.7,  82.4,  7. ],
      [ 91. ,  84.6,  33.5,  23.6,  53.1, ...,  63.9,  70.9,  28.4,  99.4,  19.6],
      ...,
      [ 58.6,  92.7,  2. ,  28.6,  78.8, ...,  8.9,  59.4,  71.5,  42.2,  47.8],
      [  6.7,  80.6,  6.2,  81.5,  18.1, ...,  10.9,  71.7,  43.3,  38.6,  7.5],
      [ 90.7,  26.8,  89. ,  46.9,  64.4, ...,  6.6,  16.6,  54.3,  70.6,  46.9],
      [ 74.7,  85.1,  54.5,  24. ,  75.9, ...,  16.2,  89.5,  99.7,  83.3,  0.3],
      [ 60.6,  0.7,  52.6,  23.4,  42.7, ...,  66.9,  15.5,  88.1,  57.5,  55.5]])
>>> z_rc
array([[ 80,  40,  60,  20, 100, ...,  60,  60,  20,  80,  20],
      [100,  60,  20,  40, 100, ...,  60, 100,  40, 100,  60],
      [ 80,  60,  80,  60,  40, ..., 100, 100,  40, 100,  20],
      [ 60, 100,  60,  20, 100, ...,  60,  60,  80, 100,  20],
      [100, 100,  40,  40,  60, ...,  80,  80,  40, 100,  20],
      ...,
      [ 60, 100,  20,  40,  80, ...,  20,  60,  80,  60,  60],
      [ 20, 100,  20, 100,  20, ...,  20,  80,  60,  40,  20],
      [100,  40, 100,  60,  80, ...,  20,  20,  60,  80,  60],
      [ 80, 100,  60,  40,  80, ...,  20, 100, 100, 100,  20],
      [ 80,  20,  60,  40,  60, ...,  80,  20, 100,  60,  60]])
>>>

You can try it yourself... obviously, I kept the reclass classes as the upper value of the reclass range to make it easy to see...

ie

71.8 is < 80 (but > 60)

39.9  is < 40 (but > 20)

58.4 is < 60 (but > 40)

ad nauseum

NOTE  your choicelist is what you want them reclassed to 'could be any thing'

>>> choicelist = ['I','am','so','confused','Dan']
>>> z_rc2 = np.select(condlist, choicelist)
>>> z_rc2
array([['confused', 'am', 'so', 'I', 'Dan', ..., 'so', 'so', 'I', 'confused', 'I'],
       ['Dan', 'so', 'I', 'am', 'Dan', ..., 'so', 'Dan', 'am', 'Dan', 'so'],
       ['confused', 'so', 'confused', 'so', 'am', ..., 'Dan', 'Dan', 'am', 'Dan', 'I'],
       ['so', 'Dan', 'so', 'I', 'Dan', ..., 'so', 'so', 'confused', 'Dan', 'I'],
       ['Dan', 'Dan', 'am', 'am', 'so', ..., 'confused', 'confused', 'am', 'Dan', 'I'],
       ..., 
       ['so', 'Dan', 'I', 'am', 'confused', ..., 'I', 'so', 'confused', 'so', 'so'],
       ['I', 'Dan', 'I', 'Dan', 'I', ..., 'I', 'confused', 'so', 'am', 'I'],
       ['Dan', 'am', 'Dan', 'so', 'confused', ..., 'I', 'I', 'so', 'confused', 'so'],
       ['confused', 'Dan', 'so', 'am', 'confused', ..., 'I', 'Dan', 'Dan', 'Dan', 'I'],
       ['confused', 'I', 'so', 'am', 'so', ..., 'confused', 'I', 'Dan', 'so', 'so']], 
      dtype='<U11')
DanPatterson_Retired
MVP Emeritus

I just ran a comparison with the conventional load and a memmap version.

It made no difference on outputs with the np.select method.

Since I don't have memory problems (in the gis sense ) both worked without incident

The arrays are largish but not huge (7000 x 14000)

>>> z = np.load("f:/temp/z_7k_14k.npy")
>>> condlist = [z<20, z<40, z<60, z<80, z>=80]
>>> choicelist = [20,40,60,80,100]
>>> z_rc2 = np.select(condlist, choicelist)
>>> z_rc2.shape, z_rc2.ndim
((7000, 14000), 2)
>>> z_rc2
array([[ 40,  40, 100,  40,  20, ...,  20,  60,  80,  40,  60],
      [ 60,  40,  60,  20,  60, ...,  40,  40, 100,  60, 100],
      [ 40,  60,  80,  80,  40, ..., 100,  20, 100,  40,  40],
      [ 20,  40,  40,  60,  20, ...,  80,  60,  20,  40,  20],
      [ 40,  80,  80,  60,  60, ...,  40, 100, 100,  60,  80],
      ...,
      [ 80,  80,  80,  20,  40, ..., 100,  60,  20,  80,  20],
      [ 40,  20,  20, 100,  40, ..., 100,  20,  40,  60,  20],
      [ 80, 100,  40,  80,  40, ...,  40,  40,  80,  20,  40],
      [ 80,  40,  40, 100,  20, ..., 100,  60,  40,  60, 100],
      [ 20,  20,  80,  60,  80, ...,  20,  40, 100,  20,  80]])
>>>
>>>
>>> zm = np.memmap("f:/temp/z_7k_14k.npy")
>>> zm_rc = np.select(condlist,choicelist)
>>> zm_rc.shape,zm_rc.ndim
((7000, 14000), 2)
>>> zm_rc
array([[ 40,  40, 100,  40,  20, ...,  20,  60,  80,  40,  60],
      [ 60,  40,  60,  20,  60, ...,  40,  40, 100,  60, 100],
      [ 40,  60,  80,  80,  40, ..., 100,  20, 100,  40,  40],
      [ 20,  40,  40,  60,  20, ...,  80,  60,  20,  40,  20],
      [ 40,  80,  80,  60,  60, ...,  40, 100, 100,  60,  80],
      ...,
      [ 80,  80,  80,  20,  40, ..., 100,  60,  20,  80,  20],
      [ 40,  20,  20, 100,  40, ..., 100,  20,  40,  60,  20],
      [ 80, 100,  40,  80,  40, ...,  40,  40,  80,  20,  40],
      [ 80,  40,  40, 100,  20, ..., 100,  60,  40,  60, 100],
      [ 20,  20,  80,  60,  80, ...,  20,  40, 100,  20,  80]])
>>>
>>> zm_rc == z_rc2
array([[ True,  True,  True,  True,  True, ...,  True,  True,  True,  True,  True],
      [ True,  True,  True,  True,  True, ...,  True,  True,  True,  True,  True],
      [ True,  True,  True,  True,  True, ...,  True,  True,  True,  True,  True],
      [ True,  True,  True,  True,  True, ...,  True,  True,  True,  True,  True],
      [ True,  True,  True,  True,  True, ...,  True,  True,  True,  True,  True],
      ...,
      [ True,  True,  True,  True,  True, ...,  True,  True,  True,  True,  True],
      [ True,  True,  True,  True,  True, ...,  True,  True,  True,  True,  True],
      [ True,  True,  True,  True,  True, ...,  True,  True,  True,  True,  True],
      [ True,  True,  True,  True,  True, ...,  True,  True,  True,  True,  True],
      [ True,  True,  True,  True,  True, ...,  True,  True,  True,  True,  True]], dtype=bool)
>>> np.allclose(zm_rc,z_rc2)
>>>
>>> np.sum( (zm_rc == z_rc2).astype('int') ) # returns the row*col product
98000000
>>> np.sum( (~(zm_rc == z_rc2)).astype('int') ) # returns the equality sum... notice the '~' to negate and the extra brackets
0
>>>

>>>

RebeccaStrauch__GISP
MVP Emeritus

Summary (final solution and comments):

Since this  is a long thread with many detailed suggestions/demos, I wanted to provide a summary that’s a bit more simplified (I hope).  Dan gets the points for the solution, the last key point being I needed an integer array (which is the one I marked), but other comments helped in the final solution.

The final script that worked used a combination of many of the comments/suggestions offered by Dan, and a couple items that were necessary, even though it was suggested to do otherwise.  Just to reminder from initial post, using a pure arcpy process worked for previous projects (array size 6262 x 7352  vs. this project’s array 13035 x 6657….about 25m cellsize for both, not that that makes a difference in processing).  Software: ArcGIS 10.3.1 with default python install (Dan mentioned that python 3.4 and Pro might work better, but my choice is 10.3.1 at this time).

Key things in the final (and leaving any one of these off didn’t seem to work):

  • Using numpy.
  • Converted initial raster to an integer NumPyArray  (this was ok  for my needs)
    • arr = (arcpy.RasterToNumPyArray(inRas)).astype(int)
  • Split the numpy array into smaller arrays (I split into thirds), then vstack them back together
  • Use the LowerLeft in Cellsize information, grabbed from the initial raster, when creating the new raster (and save it).  Do NOT apply the projection until after converting to polygon, even though it should be ok, seems to have an issue “sometimes”
    • newRas = arcpy.NumPyArrayToRaster(arr_stack, lowerLeft, cellSize) 
    • newRas.save(outDEMtmp) 
  • Reset the temporary arrays to [] as soon as they are not needed to keep the ArcMap memory clear. Although Dan mentioned that it should be “del” in many cases, looking at other threads and repeated testing (including running the script multiple times in same session), shows that only “anArray = []” seems to work.  I added a “del” line at the end for cleanup, but doesn’t seem to make a difference.
  • Use the arr_reclass() function provided in the demo by Dan in his blog (I modified slightly for my needs)  https://community.esri.com/blogs/dan_patterson/2016/03/12/reclassify-raster-data-simply
  • NOTE: although irrelevant to the testing and solution to this issue, I meant to use my raster that I had creating with values in feet, vs the inRas.  This wasn’t critical for testing, but would have been in my final.  Correct raster works too….and looks much better.

The script below works as needed, and I can run it multiple times (three for sure) without running a memory issue now.,,,but on the 4th attempt, get same error, so something is still holding on to memory.  Just fyi.

Hope this will help others.

import time
import os
import arcpy
from arcpy import env
from arcpy.sa import *
import numpy as np

arcpy.CheckOutExtension("spatial")

# Function by https://community.esri.com/people/Dan_Patterson with slight modifications
def arr_reclass(a, bins=[], new_bins=[], mask=False, mask_val=None):
  """a    - integer point array to be reclassed using
     bins - sequential list/array of the lower limits of each class
        include one value higher to cover the upper range.
     array dimensions will be squeezed 
  """ 
  a_rc = np.zeros_like(a)
  if (len(bins) < 2): # or (len(new_bins <2)):
    print("Bins = {} new = {} won't work".format(bins,new_bins))
    return a
  if len(new_bins) < 2:
    new_bins = np.arange(1,len(bins)+2)   
  new_classes = zip(bins[:-1],bins[1:],new_bins)
  for rc in new_classes:
    q1 = (a >= rc[0])
    q2 = (a < rc[1])
    z = np.where(q1 & q2, rc[2],0)
    a_rc = a_rc + z
  return a_rc

theWorkspace = r'C:\___bear2016\_Skwentna2017\Skwentna_prep.gdb'
arcpy.env.workspace = theWorkspace
arcpy.env.overwriteOutput = True
wsPath, wsGDB = os.path.split(theWorkspace)
print("Project gdb (vector, i.e. for points, lines, polygons): {0}".format(theWorkspace))

inStudy = r'C:\___bear2016\_Skwentna2017\Skwentna_prep.gdb\Study16b'
stdyPath, stdyName = os.path.split(inStudy)
print("Study boundary: {0}".format(stdyName))
spatialref = arcpy.Describe(inStudy).spatialReference

rasterIn = r'C:\___bear2016\_Skwentna2017\Skwentna_raster.gdb\elevExtent'
inRas = arcpy.Raster(rasterIn)
rasterGDB, rasterFile = os.path.split(rasterIn)
# grab the lowerleft and cellsize info, to be use when creating new raster
lowerLeft = arcpy.Point(inRas.extent.XMin, inRas.extent.YMin)
cellSize = inRas.meanCellWidth
noDataReclass = 999999      # used when creating reclass "bins"

#  My variables...typically assigned by tool input
elevUnits = "meters"
elevCutoffFt = "6000"
flatElevBreak = "2000"
elevBreaks = flatElevBreak.split("; ")
elevBreaks.sort(key=float)  
flatElevBreak = ";".join(elevBreaks)
flatElevBreak
'2000'

# final output
pElevClass = arcpy.os.path.join(theWorkspace, "elevReclassPoly")  
outDEMtmp = arcpy.os.path.join(theWorkspace, "elevReclass")  

# for elev values in meters, creates vursion in feet for processing
if elevUnits == "meters":   #  3.280839895013123
  toElevFeet = 3.28084
else:      # note, IF had to go from feet to meters, factor would be 0.3048
  toElevFeet = 1
elevFeet = (inRas * toElevFeet)
saveElevFeet = arcpy.os.path.join(rasterGDB, "elevFt")
elevFeet.save(saveElevFeet)
elevMin = elevFeet.minimum  - 1
elevMax = elevFeet.maximum  + 1

# Creates the from (bin) and to (new_bin) values for reclass function
#   uses the input elev breaks, min/max, and nodata variable
elevSplit = flatElevBreak.split(";")
lstCnt = 1
bins = []
new_bins = []
for reclassValue in elevSplit:
  if lstCnt == 1:
    bins.append(int(elevMin))
    bins.append(int(-1))
    bins.append(int(reclassValue))
    new_bins.append(int(9))
    new_bins.append(int(reclassValue))
  else:
    bins.append(int(reclassValue))    
    new_bins.append(int(9))    
    new_bins.append(int(reclassValue))    
  lstCnt += 1
  prevValue = reclassValue

bins.append(int(elevCutoffFt))
bins.append(int(elevMax))
bins.append(int(noDataReclass))
new_bins.append(int(elevCutoffFt))
new_bins.append(9)  
new_bins.append(int(noDataReclass))
new_classes = zip(bins[:-1],bins[1:],new_bins)

print("   ...Old Ranges:\n   {0}".format(bins))
print("   ...New values:\n   {0}".format(new_bins))
print("   ...New classes:\n   {0}".format(new_classes))

# Creats integer NumPy array
arr = (arcpy.RasterToNumPyArray(saveElevFeet)).astype(int)   #, nodata_to_value=9999999)
print(arr)
print(arr.ndim)
print(arr.shape)
#(13035, 6657)

# saves the array to disc...can they free up arr memory, or use in different process, if needed
np.save(r"C:\___bear2016\_Skwentna2017\outras.npy", arr)   #save in binary format

# This works to free up the memory in this script....similar lines later in script...
arr = []

# loads the array back in from file...could have just reused "arr", but as precaution during testing...
arrLoad = np.load(r"C:\___bear2016\_Skwentna2017\outras.npy")

# count number of elements in array, and determin record numbers to split into thirds
arrCnt = len(arrLoad)
arrSplit0 = int(arrCnt/3)
arrSplit1 = (int(arrCnt/3) * 2)
# now reclass in slices and stack back  
arr_rc0 = arr_reclass(arrLoad[:arrSplit0,:], bins=bins, new_bins=new_bins)  
arr_rc1 = arr_reclass(arrLoad[arrSplit0:arrSplit1,:], bins=bins, new_bins=new_bins)  
arr_rc2 = arr_reclass(arrLoad[arrSplit1:,:], bins=bins, new_bins=new_bins)  
# clear arrLoad
arrLoad = []
# merge the thirds into new array
arr_stack = np.vstack((arr_rc0,arr_rc1, arr_rc2))  
#print(arr_stack.ndim)
#print(arr_stack.shape)

# clear the 3 arrays
arr_rc0 = []
arr_rc1 = []
arr_rc2 = []

#print(lowerLeft)
# create new raster from merged array, using inital lowerlet and cellsize
newRas = arcpy.NumPyArrayToRaster(arr_stack, lowerLeft, cellSize)  
newRas.save(outDEMtmp)  
arr_stack = []

# This was my goal...get reclassed raster
arcpy.RasterToPolygon_conversion(outDEMtmp, pElevClass, "SIMPLIFY", "Value")

# assigning initial spatialref, since numpy doesn't reserve or care
arcpy.DefineProjection_management(pElevClass, spatialref)
arcpy.DefineProjection_management(outDEMtmp, spatialref)

# added this for cleanup....but makes no difference whether there or not.
del arr, arr_rc0, arr_rc1, arr_rc2, arr_stack, arrLoad
DanPatterson_Retired
MVP Emeritus

Great!  glad you perservered in true numpythion fashion... it looks good...Is it safe now to delete the 48 GB of test data???

0 Kudos
RebeccaStrauch__GISP
MVP Emeritus

lol.  Yes, you bet Dan.  Thanks again for all your help over the weekend.  It's nice to get this resolved.  Now the fun part of merging into my tool/script to see if it still works.  If not, may need to break the tool in two.

...still which I could get the memory to truly clear....but chance of me worrying about running this 4x in a row.....I can work with that. 

0 Kudos