topic Re: reclassify using Numpy, out of memory in Python Questions

reclassify using Numpy, out of memory

RebeccaStrauch__GISP — Sat, 11 Dec 2021 16:03:34 GMT

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 Python python snippets Imagery and Remote Sensing

Re: reclassify using Numpy, out of memory

RebeccaStrauch__GISP — Thu, 19 May 2016 23:20:40 GMT

BTW - I ran the script above on the other study area, and it ran without any issue, i.e. no out of memory error, and I was able to use

test20a = arcpy.NumPyArrayToRaster(arr_rc)

to get it back into a raster, so that makes me believe that is it just an issue with the size of the raster (Col x row), so again, maybe there is an easy way for me to process half at a time, then merge them back....without having overlap/gap issues? That is the approach I will take next if I don't get any better suggestions/options. thanks again for looking.

Re: reclassify using Numpy, out of memory

DanPatterson_Retired — Thu, 19 May 2016 23:49:02 GMT

quick questions

python version (please say 3.4)
numpy version (please say 1.9.x)

comment

how many rasters are you saving (count them) before you actually get to the array reclass???!!!
too many... they are still in memory and references to them don't get deleted until AFTER the array is reclassed

solution

prepare your array for reclassification or whatever, then do
np.save('c:/test/your_data.npy', your_array) # saves in binary format
in_array = np.load('c:/test/your_data.npy')

Homework

>>> np.save('/tmp/123', np.array([[1, 2, 3], [4, 5, 6]]))

>>> a = np.load('/tmp/123.npy')

>>> a

array([[1, 2, 3],

[4, 5, 6]])

Now if memory really sucks, you can read in bits

Mem-map the stored array, and then access the second row

directly from disk:

>>> X = np.load('/tmp/123.npy', mmap_mode='r')

>>> X[1, :]

memmap([4, 5, 6])

I will be up most of the night writing, so get back to me.

And if we can't solve this, you may have to MOVE it to python, since it won't be 'seen' unless they have access to the N.R

Re: reclassify using Numpy, out of memory

RebeccaStrauch__GISP — Thu, 19 May 2016 23:59:07 GMT

quick questions

python version (please say 3.4) Default ArcGIS 10.3.1, so assuming 2.7
numpy version (please say 1.9.x) same, so ??

comment

how many rasters are you saving (count them) before you actually get to the array reclass???!!! only dealing with one raster/array...size specified in OP
too many... they are still in memory and references to them don't get deleted until AFTER the array is reclassed Just one reclass, and it crashes...rebooted to clear memory,

solution
prepare your array for reclassification or whatever, then do
np.save('c:/test/your_data.npy', your_array) # saves in binary format
in_array = np.load('c:/test/your_data.npy')
Homework
>>> np.save('/tmp/123', np.array([[1, 2, 3], [4, 5, 6]]))
>>> a = np.load('/tmp/123.npy')
>>> a
array([[1, 2, 3],
[4, 5, 6]])
Now if memory really sucks, you can read in bits
Mem-map the stored array, and then access the second row
directly from disk:
>>> X = np.load('/tmp/123.npy', mmap_mode='r')
>>> X[1, :]
memmap([4, 5, 6])

I'll give that a try. Thanks for responding Dan.

Re: reclassify using Numpy, out of memory

RebeccaStrauch__GISP — Sat, 11 Dec 2021 16:03:37 GMT

I did the "homework" so see how it should work, but when doing it with the real data, still getting memory issues

.....
arr = arcpy.RasterToNumPyArray(inRas, nodata_to_value=9999999)
np.save(r"C:\___bear2016\_Skwentna2017\outras.npy", arr)   #save in binary format
in_array = np.load(r"C:\___bear2016\_Skwentna2017\outras.npy")
mask = False
mask_val = None
print(arr)
#arr_rc = arr_reclass(arr, bins=bins, new_bins=new_bins)
arr_rc = arr_reclass(in_array, bins=bins, new_bins=new_bins)

So, I think my best bet work be to find a "clean" and automated (since this will be done in different study areas) to split, reclass, con, and then convert to poly. I'll go ahead and move this to the python space to get more input. (I like the design of you numpy space...very organized!)

I'm not sure how much more I'll be able to work on this until Monday, but will try if I find anything that looks promising.

Re: reclassify using Numpy, out of memory

DanPatterson_Retired — Sat, 11 Dec 2021 16:03:40 GMT

It works... here is the script and how I created a floating point array in the range 0-100.0 and reclassed to 8 bins of 1 to 8.

# coding: utf-8
"""
Script:  reclass_array_demo.py
Author:  Dan.Patterson@carleton.ca
Modified: 2016-05-19
Purpose: For Rebecca... array has no, nodata to simplify
References:
    RasterToNumPyArray and NumPyArrayToRaster can be used to and from
    raster formats of different types.
"""
import numpy as np
from textwrap import dedent
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

if __name__=="__main__":
    """For Rebecca  Uncomment the #z and #np.save lines to create the array, then you can load it blah blah"""
    r = 7000
    c = 14000
    #z = np.abs(100 * np.random.random_sample((r,c)) -100)  # make an array with rand 0-100
    #np.save("f:/temp/z_7k_14k.npy",z)
    z = np.load("f:/temp/z_7k_14k.npy")
    bins = [0,5,10,15,20,25,30,60,100]
    new_bins = [1, 2, 3, 4, 5, 6, 7, 8]
    mask = False
    mask_val = None
    a_rc = arr_reclass(z, bins=bins, new_bins=new_bins)
    np.save("f:/temp/reclass_z_7k_14k.npy",a_rc)

So it works, but it too about 2 seconds to complete, so I may have to see if I can vectorize it some more to keep within in my one coffee sip requirement

And here are the results which are as expected for the classes given and the sampling framework

>>> a_rc

array([[ 7., 6., 8., ..., 8., 7., 7.],

[ 7., 7., 7., ..., 8., 7., 8.],

[ 7., 7., 8., ..., 8., 7., 7.],

...,

[ 8., 8., 7., ..., 8., 4., 6.],

[ 8., 7., 7., ..., 7., 7., 8.],

[ 2., 1., 8., ..., 8., 3., 8.]])

>>> a_rc.shape

(7000, 14000)

>>> np.histogram(a_rc,bins=[1, 2, 3, 4, 5, 6, 7, 8])

(array([ 4903001, 4898616, 4900343, 4899584, 4898414, 4901199, 68598843]), array([1, 2, 3, 4, 5, 6, 7, 8]))

>>>

array size on disk

Re: reclassify using Numpy, out of memory

DanPatterson_Retired — Fri, 20 May 2016 00:33:59 GMT

See my comments about saving and loading the array from disk. Arcpy is notorious for stealing and not giving back. Once you have your raster converted to an array, you don't need arcmap any more

Re: reclassify using Numpy, out of memory

RebeccaStrauch__GISP — Fri, 20 May 2016 01:07:36 GMT

Thanks for that Dan, but even creating your random array, saving etc (and ~765 MB...2x the size of mine), I'm getting a memory error, although different than the "out_of_memory" I was getting before

Runtime error
Traceback (most recent call last):
File "<string>", line 51, in <module>
File "<string>", line 33, in arr_reclass
MemoryError

My desktop is older, but semi-beefy...better than what most get. I may try to run your script on my home machine to see if it works.

I still think my best bet will be to split the array "in half" or thirds, process each separately, then merge them back together. Just thinking out loud, but my guess that this would be a simple thing to do, once it is in an array. Query the number of items, based on the size, decide how many parts, write each part to a file, load one at a time (overwriting the memory space) to do the reclass, write it back to disk, process the next. Once done, stitch them back together and get it back into a raster for processing. what do you think??

Also, is there aeasy what to flush the memory if using this in a tool? This will be part of my Toolbox (and maybe addin), run by various users and areas, so making it work...even if in chunks..is important. I've tried using the function clearINMEM() by James Crandall, but my guess this isn't actually using the "IN_MEMORY" space ??

(time for me to head home...)

Re: reclassify using Numpy, out of memory

Luke_Pinner — Fri, 20 May 2016 01:25:22 GMT

Rebecca, the RasterToNumpyArray function supports extracting smaller chunks from the raster using the {lower_left_corner}, {ncols}, {nrows} parameters:

RasterToNumPyArray (in_raster, {lower_left_corner}, {ncols}, {nrows}, {nodata_to_value})

So you could loop through and reclass a smaller array each time.

Re: reclassify using Numpy, out of memory

DanPatterson_Retired — Fri, 20 May 2016 01:33:55 GMT

Luke that is true, but you can keep arcpy out of it using mem-map which is essentially the same.

Mem-map supports 2 GB of memory, so I am not sure if it because I am using python 3.4.x and numpy 1.9.x for these tests (just using the ArcGIS Pro 1.2 install versions... nothing fancy this time)

I might be worth a while, since you have the array saved to disk, you can tile it, horizontally just to make stacking easier when you are done. I will post an example later

From the help topic in np.load.....

Mem-map the stored array, and then access the second row

directly from disk:

>>> X = np.load('/tmp/123.npy', mmap_mode='r')

>>> X[1, :]

memmap([4, 5, 6])

Re: reclassify using Numpy, out of memory

DanPatterson_Retired — Sat, 11 Dec 2021 16:03:42 GMT

It may be split time...

Now... you can create the array, save it to disk and load 'pieces', reclassify, save to disk, repeat then reassemble

I have no clue why you are running out of memory. The arrays aren't that big, the processing within numpy barely makes a glitch time wise or memory wise.

I would suggest you monitor your memory useage. How much ram do you have anyway?

Anyway, the below is a verbose description with results so you can see the process conceptually.

Later...

replace the below with, keeping the reclass function....

if __name__=="__main__":
    """run demo"""
    r = 8
    c = 10
    #  uncomment the 2 lines below to create the array only once, then recomment
    #  for subsequent testing
    #z = np.abs(100 * np.random.random_sample((r,c)) - 100)  # make rand value in 0-100 
    #np.save("z_8x10.npy",z) 
    z = np.load("z_8x10.npy")
    bins = [0, 10, 20, 30, 40, 60, 70, 80, 90, 100]
    new_bins = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
    z_rc = arr_reclass(z, bins=bins, new_bins=new_bins)
    z_hist = np.histogram(z_rc,bins=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
    result = np.array(list(zip(z_hist[1],z_hist[0])))
    print("input array\n{}\nreclassed array\n{}".format(z, z_rc))
    print("histogram....\n(class, count)\n{}".format(result))
    # now reclass in slices and stack back
    z_rc0 = arr_reclass(z[:4,:], bins=bins, new_bins=new_bins)
    z_rc1 = arr_reclass(z[4:,:], bins=bins, new_bins=new_bins)
    z_stack = np.vstack((z_rc0,z_rc1))
    frmt="""
    Now split the array into some bits, do the reclass and reformulate
    If you want to use memmap (see np.load), it just reads the bits from disk
    If you want, you can reuse the names to clear memory...this was not done
    in this case.  you can save intermediate results to disk as well.
    ....
    first split and reclass...
    {}
    second...
    {}
    all equal? np.allclose(z_stack, z_rc) = {}
    """
    print(dedent(frmt).format(z_rc0,z_rc1, np.allclose(z_stack, z_rc)))

The results

input array
[[ 57.04  17.38  74.4  18.04  25.    99.25  18.09  9.79  45.35  67.94]
[ 28.73  74.36  32.13  60.14  0.79  14.74  60.81  20.31  53.26  89.08]
[ 26.97  70.29  96.59  43.5  97.16  96.33  18.    82.99  35.37  40.31]
[  0.12  52.13  26.93  19.59  48.93  20.31  4.3  24.43  29.33  20.46]
[ 85.84  12.67  94.1  43.57  76.68  22.15  93.13  66.64  61.31  49.3 ]
[ 55.35  25.44  23.02  62.98  12.44  83.66  42.53  92.32  48.41  24.  ]
[ 70.27  45.64  3.    42.58  73.27  66.55  1.68  67.29  48.39  2.47]
[ 96.38  4.19  78.12  82.1  13.83  73.13  52.77  28.02  57.7    5.44]]
reclassed array
[[ 5.  2.  7.  2.  3.  9.  2.  1.  5.  6.]
[ 3.  7.  4.  6.  1.  2.  6.  3.  5.  8.]
[ 3.  7.  9.  5.  9.  9.  2.  8.  4.  5.]
[ 1.  5.  3.  2.  5.  3.  1.  3.  3.  3.]
[ 8.  2.  9.  5.  7.  3.  9.  6.  6.  5.]
[ 5.  3.  3.  6.  2.  8.  5.  9.  5.  3.]
[ 7.  5.  1.  5.  7.  6.  1.  6.  5.  1.]
[ 9.  1.  7.  8.  2.  7.  5.  3.  5.  1.]]
histogram....
(class, count)
[[ 1  9]
[ 2  9]
[ 3 14]
[ 4  2]
[ 5 17]
[ 6  8]
[ 7  8]
[ 8  5]
[ 9  8]]


Now split the array into some bits, do the reclass and reformulate
If you want to use memmap (see np.load), it just reads the bits from disk
If you want, you can reuse the names to clear memory...this was not done
in this case.  you can save intermediate results to disk as well.
....
first split and reclass...
[[ 5.  2.  7.  2.  3.  9.  2.  1.  5.  6.]
[ 3.  7.  4.  6.  1.  2.  6.  3.  5.  8.]
[ 3.  7.  9.  5.  9.  9.  2.  8.  4.  5.]
[ 1.  5.  3.  2.  5.  3.  1.  3.  3.  3.]]
second...
[[ 8.  2.  9.  5.  7.  3.  9.  6.  6.  5.]
[ 5.  3.  3.  6.  2.  8.  5.  9.  5.  3.]
[ 7.  5.  1.  5.  7.  6.  1.  6.  5.  1.]
[ 9.  1.  7.  8.  2.  7.  5.  3.  5.  1.]]
all equal? np.allclose(z_stack, z_rc) = True

Re: reclassify using Numpy, out of memory

Luke_Pinner — Fri, 20 May 2016 03:23:37 GMT

I tend to read rasters directly to arrays in memory as chunks both because of the 2GB memmap limit and because it avoids the doubling up of file system IO (reading it in, saving it back out to npy, reading it back in). But I'm usually dealing with large stacks of timeseries data. Horses for courses

Re: reclassify using Numpy, out of memory

DanPatterson_Retired — Fri, 20 May 2016 04:07:28 GMT

Luke are you doing data prep and processing outside of Arc* ? or are you 'tiling' or equivalent when you mention chunks? In this supposed world of big data, there seems little obvious interest in the data analysis side in the gis field. Most complaints are the 1000 record limits for AGOL

Re: reclassify using Numpy, out of memory

RebeccaStrauch__GISP — Fri, 20 May 2016 05:13:31 GMT

Dan, I'll take a look at this tomorrow. I did test your other script, i.e. creating the random array, on my home win10 machine. Again, 10.3.1 gave the memory error. I then ran it in Pro (using the file written before the other crashed) and it worked, so that does point to the version of Python and/or numpy being an issue. However, I do need to make this work in ArcMap, not Pro, so the split option seems like a better solution.

Luke, I think the speed of the numpy or similar is the way for me to go, especially if I can quickly split the array and remap. With the arcpy remap taking +4 hours of running before it ran into issues with memory, at least the numpy version let me know in a matter of second (< min). I had looked at chunking it up with the options available in RasterToNumpyArray, and that would probably work, if I figure out the logic to find the lowerleft corners of each split...and that might still be an option. I was playing areound with the SpitRaster command, but my first attempt was a disaster since I messed up the attributes. Anyway, I appreciate the effort, and in the long run I'll probably use bits of all the suggestions.

Re: reclassify using Numpy, out of memory

DanPatterson_Retired — Fri, 20 May 2016 10:22:36 GMT

yes it is python version... split works too

Re: reclassify using Numpy, out of memory

RebeccaStrauch__GISP — Sat, 11 Dec 2021 16:03:45 GMT

The splitting of the array has solved the memory issues I believe. FYI - for others that may see this, to make it a little more generic/robust for different size raster, I'm splitting it into 3 raster, based on len(arr)

...
arr = np.load(r"C:\___bear2016\_Skwentna2017\outras.npy")
# count number of elements so it can be split
arrCnt = len(arr)
arrSplit0 = int(arrCnt/3)
arrSplit1 = (int(arrCnt/3) * 2)
# now reclass in slices and stack back  
arr_rc0 = arr_reclass(arr[:arrSplit0,:], bins=bins, new_bins=new_bins)  
arr_rc1 = arr_reclass(arr[arrSplit0:arrSplit1,:], bins=bins, new_bins=new_bins)  
arr_rc2 = arr_reclass(arr[arrSplit1:,:], bins=bins, new_bins=new_bins)  
arr_stack = np.vstack((arr_rc0,arr_rc1, arr_rc2))  
arr = arcpy.NumPyArrayToRaster(arr_stack) #, value_to_nodata=noDataReclass)

I am having an issue yet trying to save it back to the raster, but Dan, unless you or someone can see quickly what the issue is, I'll keep searching/testing an post another thread if necessary. I'm not sure if it is a numpy issue or the arcpy issues on

...closing...will continure editing...

EDIT: last response was getting to slow, so closed and editing to complete response....

error message I'm getting

Runtime error

Traceback (most recent call last):
File "<string>", line 127, in <module>
File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\__init__.py", line 2292, in NumPyArrayToRaster
return _NumPyArrayToRaster(*args, **kwargs)
TypeError: Cannot create raster for numpy array.

I have tried a few varieties...but now I may have to look at the no data. I had set a nodata value, which should have handled it....but, need to look again..

Re: reclassify using Numpy, out of memory

DanPatterson_Retired — Sat, 11 Dec 2021 16:03:48 GMT

quickly offhand, don't use arr as the variable name for output since it 'exists' try arr_out ... I will check when on an arc machine

BACK

I just assembled a test using 3 arrays (floating point) and stacked them

Dan demos in python 3.4
Date...
>>> a_rc
array([[ 7.,  6.,  8., ...,  8.,  7.,  7.],
       [ 7.,  7.,  7., ...,  8.,  7.,  8.],
       [ 7.,  7.,  8., ...,  8.,  7.,  7.],
       ...,
       [ 8.,  8.,  7., ...,  8.,  4.,  6.],
       [ 8.,  7.,  7., ...,  7.,  7.,  8.],
       [ 2.,  1.,  8., ...,  8.,  3.,  8.]])
>>> t = np.vstack((a_rc,a_rc,a_rc))
>>> t
array([[ 7.,  6.,  8., ...,  8.,  7.,  7.],
       [ 7.,  7.,  7., ...,  8.,  7.,  8.],
       [ 7.,  7.,  8., ...,  8.,  7.,  7.],
       ...,
       [ 8.,  8.,  7., ...,  8.,  4.,  6.],
       [ 8.,  7.,  7., ...,  7.,  7.,  8.],
       [ 2.,  1.,  8., ...,  8.,  3.,  8.]])
>>>
>>> import arcpy
>>> tt = arcpy.NumPyArrayToRaster(t)
>>> tt.save("f:/test/tt")
>>> t.ndim
2
>>> t.shape
(21000, 14000)
>>>

The results

Do note the size of the input arrays and the output raster. I am just praying that the inability to save was due to the naming convention and that you don't have to save each raster, then

Mosaic To New Raster—Help | ArcGIS for Desktop . If it is a memory issue you will have to label the LL corner of each array, save it so that they can be merged inside of arcmap into a final array.

Rebecca... time to move to python 3.4 (at least)

Re: reclassify using Numpy, out of memory

DanPatterson_Retired — Sun, 22 May 2016 00:27:50 GMT

Just did some digging, python prior to 2.5 had the 2GB limit (even for 64 bit systems). Numpy does not, there were some pretty ridiculously large files being read and used

see hpaulj's comment http://stackoverflow.com/questions/37365064/memory-mapping-slows-down-over-time-alternatives

Plus there are alternatives

http://stackoverflow.com/questions/26214696/creating-very-large-numpy-arrays-in-small-chunks-pytables-vs-numpy-memmap

h5py, pytables and Dask (haven't checked these all out yet, the first is best for simple single content arrays, the 2nd for tabular arrays, Dask is from Continuum, ) Array — dask 0.9.0 documentation

Re: reclassify using Numpy, out of memory

RebeccaStrauch__GISP — Sat, 11 Dec 2021 16:03:50 GMT

Dan, I don't think the "arr" name itself the issue (since that made no difference), but I'm using a different array names just in case.

What I am doing now is redefining each array as "= []" as soon as I don't need it anymore, which seems to free up the space. From a number of other discussions I've read thru, doing a "del array_variable" only deletes the variable and does not free up the memory. (just an fyi for anyone reading this). So that help.

The next issues (not discovered or resolved in this order) were that I had to specify the lowerleft and the cell-size that I saved at the very beginning of my script from the original raster....I also had to use DefineProjection since that isn't an option in the NumPyArrayToRaster (snippet of these lines below)

np.save(r"C:\___bear2016\_Skwentna2017\outras.npy", arr)   #save in binary format
arr = []
arrLoad = np.load(r"C:\___bear2016\_Skwentna2017\outras.npy")
# count number of elements so it can be split
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)  
arr_stack = np.vstack((arr_rc0,arr_rc1, arr_rc2))  
print(arr_stack.ndim)
print(arr_stack.shape)
arr_rc0 = []
arr_rc1 = []
arr_rc2 = []
print(lowerLeft)
newRas = arcpy.NumPyArrayToRaster(arr_stack, lowerLeft, cellSize)  
newRas.save(outDEMtmp)  
newRas = []
arcpy.DefineProjection_management(outDEMtmp, spatialref)
arcpy.RasterToPolygon_conversion(outDEMtmp, pElevClass, "SIMPLIFY", "Value")

So, now the issue is that the RasterToPolygon does not see the outDEMtmp as a valid raster,

even though it appears to be fine when displaying and/or looking at the properties.

the error in the python window (which I tried first)

>>> arcpy.RasterToPolygon_conversion(outDEMtmp, pElevClass, "SIMPLIFY", "Value")
Runtime error Traceback (most recent call last): File "<string>", line 1, in <module> File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\conversion.py", line 337, in RasterToPolygon raise e ExecuteError: ERROR 000864: Input raster: The input is not within the defined domain. ERROR 000367: Invalid GP data type
>>>

Neither of the error seems to be relevant. A restart didn't help....so on to researching this. ...maybe I'll copy to a new raster, if that works. Any chance the NumPyToRaster doesn't work quite right??

Re: reclassify using Numpy, out of memory

DanPatterson_Retired — Sun, 22 May 2016 02:34:28 GMT

del delete the reference and cleans it up... if it is only being used by python and not any other software maybe like arcmap or a geodatabase being used in arcmap if the source or destination is one...

No...I ran tests today and got a valid rasters and huge ones as you can see. The issue isn't with numpy, NumPyToRaster produces a valid raster. The only thing you did different is that you defined a coordinate system... don't, just convert it to a polygon, save your work, then see if you can define the projection. If you are doing this in a geodatabase... move it outside of a geodatabase...make an esri grid in a folder, make a shapefile in a folder, I want to rule out this whole geodatabase thing and its clutches and clucking on everything that is done. You can always copy stuff back into one, I just want to work with native operating stuff first.

Now... as for the results. Are you planning to mix the results with other data? if so, then you will need the projected featureclass/shapefile, it is not needed for the raster.

keep me posted.