Summary
Converts a raster to a NumPy array.
Discussion
A Python NumPy array is designed to deal with large arrays. There are many existing Python functions that have been created to process NumPy arrays, the most noted being contained in the SciPy scientific computing package for Python. You may want to convert an ArcGIS raster to a NumPy array to
- Implement one of the many existing Python functions that can be applied to a NumPy array (for example, run filters on the data, perform multidimensional analysis, or utilize optimization routines).
- Develop a custom function by accessing the individual cells within the NumPy array (for example, to implement neighborhood notation, change individual cell values, or run accumulative operators on an entire raster).
If the array definition (the lower left corner and the number of rows and columns) exceeds the extent of the in_raster, the array values will be assigned NoData.
If the lower_left_corner does not coincide with the corner of a cell, it will automatically be snapped to the lower left of the nearest cell corner applying the same rules as the Snap Raster environment setting. This snapping action within the RasterToNumPy function is not to be confused with the Snap Raster environment setting; the function only uses the same interaction. For more information, see the following:
RasterToNumPyArray supports the direct conversion of multiband rasters to an N-dimensional array (ndarray).
- If the input Raster instance is based on a multiband raster, it returns an ndarry, where the length of the first dimension represents the number of bands. The ndarray will have the dimensions (bands, rows, columns).
- If the input Raster instance is based on a single raster or a specific band from a multiband raster, it returns a two-dimensional array with the dimensions (rows, columns).
Syntax
RasterToNumPyArray (in_raster, {lower_left_corner}, {ncols}, {nrows}, {nodata_to_value})| Parameter | Explanation | Data Type |
in_raster | The input raster to convert to a NumPy array. | Raster |
lower_left_corner | The lower left corner within the in_raster from which to extract the processing block to convert to an array. The x- and y-values are in map units. If no value is specified, the origin of the input raster will be used. (The default value is None) | Point |
ncols | The number of columns from the lower_left_corner in the in_raster to convert to the NumPy array. If no value is specified, the number of columns of the input raster will be used. (The default value is None) | Integer |
nrows | The number of rows from the lower_left_corner in the in_raster to convert to the NumPy array. If no value is specified, the number of rows of the input raster will used. (The default value is None) | Integer |
nodata_to_value | The value to assign the in_raster NoData values in the resulting NumPy array. If no value is specified, the NoData value of in_raster will be used. (The default value is None) | Variant |
| Data Type | Explanation |
| NumPyArray | The output NumPy array. |
Code sample
A raster is converted to a NumPy array to calculate the percentage of the cell value for each row of the raster. A new raster is then created.
import arcpy
import numpy
# Get input Raster properties
inRas = arcpy.Raster('C:/data/inRaster')
lowerLeft = arcpy.Point(inRas.extent.XMin,inRas.extent.YMin)
cellSize = ras.meanCellWidth
# Convert Raster to numpy array
arr = arcpy.RasterToNumPyArray(inRas,nodata_to_value=0)
# Calculate percentage of the row for each cell value
arrSum = arr.sum(1)
arrSum.shape = (arr.shape[0],1)
arrPerc = (arr)/arrSum
#Convert Array to raster (keep the origin and cellsize the same as the input)
newRaster = arcpy.NumPyArrayToRaster(arrPerc,lowerLeft,cellSize,
value_to_nodata=0)
newRaster.save("C:/output/fgdb.gdb/PercentRaster")Block process an input multiband raster and calculate cell statistics across the bands. This script converts a multiband raster to a three-dimensional NumPy array and processes the array by dividing it into data blocks. It then calculates the mean of values across the rows of the block, converts the block numpy array to raster, and recombines the bands via mosaicking. A new multiband raster is created.
# Note that, if the input raster is multiband, the data blocks will also be
# multiband, having dimensions (bands, rows, columns). Otherwise, they will
# have dimensions (rows, columns).
import arcpy
import numpy
import os
# Input raster
filein = os.path.join(os.getcwd(),r"input\input.tif")
# Output raster (after processing)
fileout = os.path.join(os.getcwd(),r"output\blockprocessingrdb22.tif")
# Size of processing data block
# where memorysize = datatypeinbytes*nobands*blocksize^2
blocksize = 512
# ----------------------------------------------------------------------------
# Create raster object from file
myRaster = arcpy.Raster(filein)
# Set environmental variables for output
arcpy.env.overwriteOutput = True
arcpy.env.outputCoordinateSystem = filein
arcpy.env.cellSize = filein
# Loop over data blocks
filelist = []
blockno = 0
for x in range(0, myRaster.width, blocksize):
for y in range(0, myRaster.height, blocksize):
# Lower left coordinate of block (in map units)
mx = myRaster.extent.XMin + x * myRaster.meanCellWidth
my = myRaster.extent.YMin + y * myRaster.meanCellHeight
# Upper right coordinate of block (in cells)
lx = min([x + blocksize, myRaster.width])
ly = min([y + blocksize, myRaster.height])
# noting that (x, y) is the lower left coordinate (in cells)
# Extract data block
myData = arcpy.RasterToNumPyArray(myRaster, arcpy.Point(mx, my),
lx-x, ly-y)
# PROCESS DATA BLOCK -----------------------------
# e.g. Calculate mean of each cell of all bands.
myData -= numpy.mean(myData, axis=0, keepdims=True)
# ------------------------------------------------
# Convert data block back to raster
myRasterBlock = arcpy.NumPyArrayToRaster(myData, arcpy.Point(mx, my),
myRaster.meanCellWidth,
myRaster.meanCellHeight)
# Save on disk temporarily as 'filename_#.ext'
filetemp = ('_%i.' % blockno).join(fileout.rsplit('.',1))
myRasterBlock.save(filetemp)
# Maintain a list of saved temporary files
filelist.append(filetemp)
blockno += 1
# Mosaic temporary files
arcpy.Mosaic_management(';'.join(filelist[1:]), filelist[0])
if arcpy.Exists(fileout):
arcpy.Delete_management(fileout)
arcpy.Rename_management(filelist[0], fileout)
# Remove temporary files
for fileitem in filelist:
if arcpy.Exists(fileitem):
arcpy.Delete_management(fileitem)
# Release raster objects from memory
del myRasterBlock
del myRaster
# ----------------------------------------------------------------------------