How to get raster image as array in Python with ArcGIS

http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/RasterToNumPyArray/000v0000012z000000/

ArcGIS 10 has the ability to write and read numPy arrays.

fmark already answered the question, but here is some example OSGEO Python code that I wrote to read a raster (tif) into a NumPy array, reclass the data and then write it out to a new tif file. You can read and write any gdal supported format.

"""
Example of raster reclassification using OpenSource Geo Python

"""
import numpy, sys
from osgeo import gdal
from osgeo.gdalconst import *


# register all of the GDAL drivers
gdal.AllRegister()

# open the image
inDs = gdal.Open("c:/workshop/examples/raster_reclass/data/cropland_40.tif")
if inDs is None:
  print 'Could not open image file'
  sys.exit(1)

# read in the crop data and get info about it
band1 = inDs.GetRasterBand(1)
rows = inDs.RasterYSize
cols = inDs.RasterXSize

cropData = band1.ReadAsArray(0,0,cols,rows)

listAg = [1,5,6,22,23,24,41,42,28,37]
listNotAg = [111,195,141,181,121,122,190,62]

# create the output image
driver = inDs.GetDriver()
#print driver
outDs = driver.Create("c:/workshop/examples/raster_reclass/output/reclass_40.tif", cols, rows, 1, GDT_Int32)
if outDs is None:
  print 'Could not create reclass_40.tif'
  sys.exit(1)

outBand = outDs.GetRasterBand(1)
outData = numpy.zeros((rows,cols), numpy.int16)


for i in range(0, rows):
  for j in range(0, cols):

    if cropData[i,j] in listAg:
        outData[i,j] = 100
    elif cropData[i,j] in listNotAg:
        outData[i,j] = -100
    else:
        outData[i,j] = 0


# write the data
outBand.WriteArray(outData, 0, 0)

# flush data to disk, set the NoData value and calculate stats
outBand.FlushCache()
outBand.SetNoDataValue(-99)

# georeference the image and set the projection
outDs.SetGeoTransform(inDs.GetGeoTransform())
outDs.SetProjection(inDs.GetProjection())

del outData

This thread discusses integration of arcobjects with python.

Perhaps the code in this sample could be adapted so that it could be called from python.

I'm not sure if there's a way to pass an array byref back to python. If there is, then IPixelBlock.PixelDatabyRef would be worth trying.

You can save your raster as an ESRI ascii grid and read/manipulate that file with numpy.

This provides some starting points: http://sites.google.com/site/davidpfinlayson2/esriasciigridformat

But watch out - it seems the ascii grid format doesn't always follow spec, so reading them properly every time can be a challenge.

I'm not sure you can manipulate the raster pixel by pixel, but you can use the geoprocessing objects in conjuntion with python API.

You can use any toolbox for that kind of manipulation. A sample script would be:

#import arcgisscripting

gp = arcgisscripting.create(9.3)

gp.AddToolbox("SA") # addint spatial analyst toolbox

rasterA = @"C:\rasterA.tif"
rasterB = @"C:\rasterB.tif"

rasterC = @"C:\rasterC.tif" # this raster does not yet exist
rasterD = @"C:\rasterD.tif" # this raster does not yet exist

gp.Minus_SA(rasterA,rasterB,rasterC)

gp.Times_SA(rasterA,rasterB,rasterD)

# lets try to use more complex functions

# lets build and expression first

expression1 = "slope( " + rasterC + ")"
expression2 = "(" + rasterC " + " rasterD + ") - " + rasterA 

gp.SingleOutputMapAlgebra_SA(expression1,@"C:\result_exp1.tif")
gp.SingleOutputMapAlgebra_SA(expression2,@"C:\result_exp2.tif")

Here is a follow up on your question. Still not possible. Not sure on version 10.0.

The easiest way would be to convert the raster to netCDF and then open that and walk through the grid. I did very much the same thing for a project involving turning rasters into feature data based on data assigned to the raster cells. I looked at this for ages, and come to the conclusion walking the grid data would be easier from netCDF.

Hope it helps.