Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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?

share|improve this question

7 Answers 7

up vote 6 down vote accepted

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

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

share|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

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.

share|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 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.

share|improve this answer

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

share|improve this answer

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.

share|improve this answer

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

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

share|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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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