RasterToNumPyArray gives an error

[GraziaZulian]

I'm having problems by using RasterToNumPyArray.

With the smal test data base it works, but if I try with a huge dataset (columns and rows: 67000, 58000, cellsize: 100 m) it gives me this error:

Traceback (most recent call last):
  File "C:\Python26\ArcGIS10.0\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 309, in RunScript
    debugger.run(codeObject, __main__.__dict__, start_stepping=0)
  File "C:\Python26\ArcGIS10.0\Lib\site-packages\pythonwin\pywin\debugger\__init__.py", line 60, in run
    _GetCurrentDebugger().run(cmd, globals,locals, start_stepping)
  File "C:\Python26\ArcGIS10.0\Lib\site-packages\pythonwin\pywin\debugger\debugger.py", line 624, in run
    exec cmd in globals, locals
  File "D:\ES_Map\recreationnew\RecNumpy2.py", line 31, in <module>
    in1 = arcpy.RasterToNumPyArray(arcpy.env.workspace + "\\eu_coastna") #arcpy.RasterToNumPyArray(input1)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\__init__.py", line 1113, in RasterToNumPyArray
    return _RasterToNumPyArray(*args, **kwargs)
RuntimeError: ERROR 999998: Unexpected Error.

could you help me please.
here is a simple example of my code

[Luke_Pinner]

That's because your raster is too big to fit in memory. A 67000*58000 (1 band) 8 bit raster is 4GB, 16 bit is 7GB and 32 bit is 14GB...

If you must use numpy, you will have to clip out smaller sections/tiles of your raster, convert to numpy array, process each tile separately and then convert back to raster and mosaic.

Looking at the commented out section of your code, you don't need to use numpy at all, all those operations can be done in spatial analyst map algebra.

[GraziaZulian]

Thanks,
I know that with the map algebra a can complete the same tasks, but I'm looking for a way to create less intermediate datasets and make my process less time consuming (I need 1 week to complete it). I'm working at continental scale with 25 m resolution data (The code I put was just a test).
Have you some suggestions?
Is there a way, for example to work with array by row?

[Luke_Pinner]

You can specify what part of the raster to extract to a numpy array - RasterToNumPyArray(in_raster, lower_left_corner, ncols, nrows).
http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/NumPyArrayToRaster/000v00000130000000/

So you could loop through and process a few rows at a time.

[GraziaZulian]

Thanks, but I can't figure out how to rebuild the output dataset.
Can you give me an example?

[Luke_Pinner]

Sorry, I'm not at work so don't have access to ArcGIS to test an example. I think the mosaic tool will do what you need though:
http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/Mosaic_To_New_Raster/00170000009800000... or http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//001700000097000000

[GraziaZulian]

Hi,
I need to go back to office as well to try.
I'll let you know tomorrow,
the question is:
Do I loop for rows and cols?
inside the loop I calculate

and then? do I create n rasters and than make the mosaic?
then I could delete the temporary rasters.... this could let me save space.

Could you give a look to the code that I'll post tomorrow?

[Luke_Pinner]

Do I loop for rows and cols?

Depends on the type of raster you are using. Some rasters are tiled internally and it is much more efficient to access by tiles or blocks. Geodatabase rasters have 128x128 pixel tiles (by default), GRIDS are 256x4 (col then row), tiffs can be tiled or untiled (safest to process one row at a time).

Here is some code to process by tile (the commented out stuff is arcpy related that I can't test as I'm not at back work until January).

ncols=67000
nrows=58000
cellsize=2
xmin=0
ymax=58000
## ^^ Get all this info from your original raster

tile_size=(128,128) #GDB rasters=(128,128), tif=(ncols,1), GRID=(256,4)
ntiles = 4 #no. of tiles to process at a time

#arcpy.env.Extent=in_raster
#arcpy.env.outputCoordinateSystem=in_raster
#arcpy.env.cellSize=in_raster
#new_raster=arcpy.CreateRasterDataset(out_path, out_name, cellsize, pixel_type, 
#                                     spatial_reference)
for row in range(0,nrows,tile_size[1]):
    for col in range(0,ncols,tile_size[0]*ntiles):
        colend=min([ncols,col+tile_size[0]*ntiles-1])
        rowend=min([nrows,row+tile_size[1]-1])
        xstart=xmin+(cellsize*col)
        ystart=ymax-(cellsize*(rowend+1))
        print 'Processing cols',col,':',colend
        print 'Processing rows',row,':',rowend
        print 'Extracting lower left corner:',xstart,ystart
        #array=RasterToNumPyArray (in_raster, arcpy.Point(xmin,ystart), ncols, step)
        #Do your processing
        #out_raster=NumPyArrayToRaster(array, arcpy.Point(xmin,ystart), cellsize, cellsize)
        #arcpy.Mosaic_management(out_raster,new_raster)

[GraziaZulian]

Thanks, I'm working on it now. I'll let you know

[OlgaChesnokova]

Hello,

I am trying to replace the holes on the Aster DEM with the values from the Ferranti DEM, and I get the same error as described in the first post. Therefore, I've decided to proceed with tiles. Now I receive the following error:

File "C:\Python27\ArcGIS10.1\Lib\site-packages\Pythonwin\pywin\framework\scriptutils.py", line 320, in RunScript
    debugger.run(codeObject, __main__.__dict__, start_stepping=1)
  File "C:\Python27\ArcGIS10.1\Lib\site-packages\Pythonwin\pywin\debugger\__init__.py", line 60, in run
    _GetCurrentDebugger().run(cmd, globals,locals, start_stepping)
  File "C:\Python27\ArcGIS10.1\Lib\site-packages\Pythonwin\pywin\debugger\debugger.py", line 655, in run
    exec cmd in globals, locals
  File "D:\tests\aster_dem\aster_dem.py", line 53, in <module>
    new_raster = arcpy.CreateRasterDataset_management(arcpy.env.workspace, "aster_fix", cellsize, pixel_type, spatialReference)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\management.py", line 11265, in CreateRasterDataset
    raise e
ExecuteError: ERROR 999999: Error executing function.
Failed to execute (CreateRasterDataset).

I am quite new to the python and arcpy, so the probability that I have some stupid mistakes inside is quite high. Any kind of help would be really appreciated!

code:
import os, sys
import string, math
import numpy
import arcpy
from arcpy.sa import *
# os.chdir(r"D:\tests\aster_dem")
arcpy.env.workspace = r"D:\tests\aster_dem"
arcpy.env.overwriteOutput = True

# input data
ferranti = arcpy.Raster(arcpy.env.workspace + "\\dhm25_ferranti_wgs84.tif")  
aster = arcpy.Raster(arcpy.env.workspace + "\\e84_aster_switzerland_clip.tif")

# spatial analyst extension check
arcpy.CheckOutExtension("Spatial")

# difference
difference = aster - ferranti
difference.save(arcpy.env.workspace + "\\difference.tif")
threshold = -300

# data description
spatialReference=arcpy.Describe(ferranti).spatialReference
inRasterDesc = arcpy.Describe(ferranti)

# coordinates of the lower left corner
xmin = inRasterDesc.Extent.XMin
xmax = inRasterDesc.Extent.XMax
ymin = inRasterDesc.Extent.YMin
ymax = inRasterDesc.Extent.YMax

# cell size, raster size
rasMeanCellHeight = inRasterDesc.MeanCellHeight
rasMeanCellWidth = inRasterDesc.MeanCellWidth
nrows = inRasterDesc.Height
ncols = inRasterDesc.Width

cellsize = rasMeanCellHeight

# tiles
tile_size=(128,128)
ntiles = 2
pixel_type = "16_BIT_UNSIGNED"
RasterType = 'TIF'

new_raster = arcpy.CreateRasterDataset_management(arcpy.env.workspace, "aster_fix", cellsize, pixel_type, spatialReference)

for row in range(0,nrows,tile_size[1]):
    for col in range(0,ncols,tile_size[0]*ntiles):
        colend = min([ncols,col+tile_size[0]*ntiles-1])
        rowend = min([nrows,row+tile_size[1]-1])
        xstart = xmin + (cellsize*col)
        ystart = ymax - (cellsize*(rowend+1))
        print 'Processing cols',col,':',colend
        print 'Processing rows',row,':',rowend
        print 'Extracting lower left corner:',xstart,ystart
        ferranti_arr = arcpy.RasterToNumPyArray(ferranti, arcpy.Point(xmin,ystart), ncols, nrows)
        aster_arr = arcpy.RasterToNumPyArray(aster, arcpy.Point(xmin,ystart), ncols, nrows)
        difference_arr = arcpy.RasterToNumPyArray(difference, arcpy.Point(xmin,ystart), ncols, nrows)

        if difference_arr[row,col] < threshold:
            aster_arr[row,col] = ferranti_arr[row,col]
       
        out_raster = NumPyArrayToRaster(aster_arr, arcpy.Point(xmin,ystart), cellsize, cellsize)
        arcpy.Mosaic_management(out_raster,new_raster)

Thank you in advance!

Regards,
Olga