How to get raster image as array in Python with ArcGIS

When starting to work with Python and ArcGIS 9.3, I assumed there would be a simple way to get a raster image into a Python array so that I can manipulate it before storing it back as another raster image. However, I can't seem to find out how to do this. Is it possible? If so, how?

-

I don't think this is possible with ArcGIS <= 9.3.1

I use the open source GDAL API for tasks such as this.

-
Great! I've used the GDAL utility programs in the past, but never thought about using them to do this. – robintw Jul 24 '10 at 11:51
I agree, the gdal Python module allows you to easily read a raster and dump the data to a Numpy array. Chris Garrard has a course on using OpenSource Python in GIS, it covers this subject. You can find it at: gis.usu.edu/~chrisg/python/2008 – DavidF Sep 1 '10 at 13:46

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)

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.

-
Thanks - that's very helpful. However, ideally I'd like to be able to iterate across the raster array doing various things to it. I would have thought there would have been a way in ArcGIS to do this, but maybe not! – robintw Jul 23 '10 at 17:13
robintw, for what i've looked in the reference, there is no way of getting a specific pixel of a raster. I'm not sure if in ArcPy (available from v10) you can fetch these individual cells, since they extended the python API with lots of new funcionality. – George Jul 23 '10 at 17:15

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.

-

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)

band1 = inDs.GetRasterBand(1)
rows = inDs.RasterYSize
cols = inDs.RasterXSize

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
``````
-
Thanks for the sample code - that's very helpful. – robintw Sep 7 '10 at 9:31

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.

-

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.

-

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

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

-
Thanks! Sadly, I'm on ArcGIS 9.3 at the moment - do you know if there is any way to do this in 9.3? – robintw Jul 23 '10 at 18:22