I'm trying to extract a small raster from a bigger one, just using python and numpy (I'm not allow to use any other software). This code just reads a from a big raster, extracts an smaller window, convert it to numpy array, convert it back to raster format and then save the small raster (Stepping stone to add more complexity later). RasterToNumpyArray seems to work fine, because I'm able to extract just a small patch from a bigger raster (8G+). So I'm running two nested loops in order to create the moving window. The running windows sweeps the source raster in (5 x 3, columns x rows). The code runs fine with the first inner loop. It generates two small rasters perfectly, but by the next rows run, I've got memory issues (type exceptions.MemoryError)
My thought is that RasterToNumpyArray doesn't loads all the raster in memory (and I don't think it does, otherwise It couln't load an 8G raster into 2G ram). My question is: what is causing the memory issues? How to solve them?
Some specs: Using ArcGIS 10.2, Operating system Microsoft Windows 7 16bit, 16G RAM, I don't know about SWAP RAM.
ps: When using DefineProjection in an empty raster (or more precisely: in a raster with all its cells null), throws an error. That's why I used try/except within the nested loops.
This is a portion of the code:
for rows in range (0,4):
for columns in range (0,5):
arcpy.AddMessage('> Attempt to extract raster [row]-[col]: [' + str(rows)+ ']-['+str(columns)+']...')
try:
small_raster_lower_left = arcpy.Point(large_raster_extent.XMin + (columns * window_size * large_raster_cellsize_X),
large_raster_extent.YMin + (rows * window_size * large_raster_cellsize_Y)
)
arcpy.AddMessage('> Raster to Numpy...')
numpy_array = arcpy.RasterToNumPyArray(large_raster,
small_raster_lower_left,
window_size,
window_size,
nodata_to_value=-99999
)
arcpy.AddMessage('> small raster coordinates. X: {0}, Y: {1}'.format(small_raster_lower_left.X,small_raster_lower_left.Y))
arcpy.AddMessage('> Numpy to raster...')
small_raster = arcpy.NumPyArrayToRaster(numpy_array,
small_raster_lower_left,
large_raster_cellsize_X,
large_raster_cellsize_Y,
value_to_nodata = -99999
)
arcpy.DefineProjection_management(small_raster, large_raster_spatial_ref)
small_raster_name = "small_dem_" + str(columns+rows*10)
arcpy.AddMessage('> Success! Saving raster: ' + small_raster_name + '...')
small_raster.save(small_raster_name)
del numpy_array
del small_raster
except:
arcpy.AddMessage('> Fault! Trying next window...')
del numpy_array
del small_raster
Two small rasters succesfully extracted form a bigger one, in color.
1 Answer 1
If your Arc license includes spatial analyst you could just use 'extract by rectangle' with a couple of iterators set up to re-define the extent object for each new extraction window. Without SA, 'Data management - clip' should be able to be used in a similar way.
del) large objects instead of allowing Python to remove them when unreferenced. This can help prevent hyper-inefficient heap utilization, where blocks of memory aren't made available until after more space is used, and the "collected" space isn't large enough for the next requested block, so even more space is used, fragmenting the heap until no space is available. You should also consider trying the 64-bit background processing module.