topic Re: RasterToNumPyArray gives an error in Python Questions

RasterToNumPyArray gives an error

GraziaZulian — Fri, 16 Dec 2011 06:30:26 GMT

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

Re: RasterToNumPyArray gives an error

Luke_Pinner — Fri, 16 Dec 2011 20:07:54 GMT

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.

Re: RasterToNumPyArray gives an error

GraziaZulian — Sat, 17 Dec 2011 05:43:21 GMT

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?

Re: RasterToNumPyArray gives an error

Luke_Pinner — Sat, 17 Dec 2011 09:53:34 GMT

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.

Re: RasterToNumPyArray gives an error

GraziaZulian — Sat, 17 Dec 2011 17:02:15 GMT

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

Re: RasterToNumPyArray gives an error

Luke_Pinner — Sat, 17 Dec 2011 20:14:10 GMT

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/001700000098000000/ or http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//001700000097000000

Re: RasterToNumPyArray gives an error

GraziaZulian — Sun, 18 Dec 2011 05:44:55 GMT

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?

Re: RasterToNumPyArray gives an error

Luke_Pinner — Sat, 11 Dec 2021 21:52:09 GMT

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)

Re: RasterToNumPyArray gives an error

GraziaZulian — Mon, 19 Dec 2011 08:32:16 GMT

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

Re: RasterToNumPyArray gives an error

OlgaChesnokova — Thu, 20 Jun 2013 09:44:49 GMT

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

Re: RasterToNumPyArray gives an error

DavidWilton — Sat, 11 Dec 2021 21:52:12 GMT

I ran into this issue so I wrote an iterator class to do this. In my case I wanted to return the X,Y,Z values but you could easily modify to return the numpy array.

It can be used simply like this:

for block in RasterBlockIterator(raster, no_data_val=self.no_data_val):
   for point in block:
      x = point[0]
      y = point[1] 
      z = point[2]‍‍‍‍‍‍‍‍‍‍

The iterator class. If you want the numpy array just return the arr on line 62. By default it works in blocks of 1000*1000

class RasterBlockIterator:
    """Iterator that extracts tiles of data from a Esri Raster."""

    def __init__(self, raster: arcpy.Raster, tile_size=(1000, 1000), no_data_val=np.nan):
        """
        :param raster: Raster to iterate
        :param tile_size: size of each block/tile to extrac to numpy array
        :param no_data_val: value in the raster which should be treated as no data
        """
        self.tile_size = tile_size
        self.raster = raster
        self.no_data_val = no_data_val
        self.cell_size = raster.meanCellHeight
        self.pos = 0

        # work out the size of all the block to create. They won't be extracted at this point
        t_rows = raster.height
        t_cols = raster.width
        x_min = raster.extent.XMin
        y_min = raster.extent.YMin
        self.tiles = []
        for row in range(0, t_rows, tile_size[1]):
            for col in range(0, t_cols, tile_size[0]):
                colend = min([t_cols, col + tile_size[0]])
                rowend = min([t_rows, row + tile_size[1]])
                xstart = x_min + (col * self.cell_size)
                ystart = y_min + (row * self.cell_size)
                self.tiles.append(namedtuple('tile', 'colstart colend rowstart rowend xstart ystart')(
                    colstart=col,
                    colend=colend,
                    rowstart=row,
                    rowend=rowend,
                    xstart=xstart,
                    ystart=ystart
                ))

    def __iter__(self):
        return self

    def __len__(self):
        return len(self.tiles)

    def __next__(self) -> List[Tuple[float, float, float]]:

        try:
            tile = self.tiles[self.pos]
        except IndexError:
            raise StopIteration

        lower_l = arcpy.Point(tile.xstart, tile.ystart)
        ncols = tile.colend - tile.colstart
        rnrows = tile.rowend - tile.rowstart
        arr = arcpy.RasterToNumPyArray(self.raster, nodata_to_value=np.nan,
                                       lower_left_corner=lower_l,
                                       ncols=ncols, nrows=rnrows)
        # Loop over rows and extract X,Y,Z
        output_data = [
            (j * self.cell_size + tile.xstart, i * self.cell_size + tile.ystart, value)
            for i, row in enumerate(arr)
            for j, value in enumerate(row)
            if not self.no_data_val == value and not np.isnan(value)
        ]

        self.pos += 1
        return output_data‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

```