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.

If it is possible, then 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.

| improve this answer | |
  • 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
  • 3
    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


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

| improve this answer | |
  • 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

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

# 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'

# 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'

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
        outData[i,j] = 0

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

# flush data to disk, set the NoData value and calculate stats

# georeference the image and set the projection

del outData
| improve this answer | |

Accessing ArcObjects from Python? 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.

| improve this answer | |

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.

| improve this answer | |

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



# lets try to use more complex functions

# lets build and expression first

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


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

| improve this answer | |
  • 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 Silva 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.

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.