¿Obtener una imagen ráster como matriz en Python con ArcGIS Desktop?

Cuando comencé a trabajar con Python y ArcGIS 9.3, supuse que habría una manera simple de obtener una imagen ráster en una matriz Python para poder manipularla antes de almacenarla como otra imagen ráster. Sin embargo, parece que no puedo descubrir cómo hacer esto.

Si es posible, ¿cómo?

No creo que esto sea posible con ArcGIS <= 9.3.1

Utilizo la API de GDAL de código abierto para tareas como esta.

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

ArcGIS 10 tiene la capacidad de escribir y leer matrices numPy.

fmark ya respondió la pregunta, pero aquí hay un ejemplo de código OSGEO Python que escribí para leer un ráster (tif) en una matriz NumPy, reclasificar los datos y luego escribirlos en un nuevo archivo tif. Puede leer y escribir cualquier formato compatible con gdal.

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

¿Acceso a ArcObjects desde Python? discute la integración de arcobjects con python.

Quizás el código en este ejemplo podría adaptarse para que pueda llamarse desde python.

No estoy seguro de si hay una manera de pasar una matriz por referencia a Python. Si lo hay, entonces valdría la pena probar IPixelBlock.PixelDatabyRef .

Puede guardar su ráster como una cuadrícula ascii de ESRI y leer / manipular ese archivo con numpy.

Esto proporciona algunos puntos de partida: http://sites.google.com/site/davidpfinlayson2/esriasciigridformat

Pero cuidado: parece que el formato de cuadrícula ascii no siempre sigue las especificaciones, por lo que leerlos correctamente cada vez puede ser un desafío.

No estoy seguro de que pueda manipular el ráster píxel por píxel, pero puede usar los objetos de geoprocesamiento junto con la API de Python.

Puede usar cualquier caja de herramientas para ese tipo de manipulación. Un script de muestra sería:

#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")

Aquí hay un seguimiento de su pregunta . Aún no es posible. No estoy seguro en la versión 10.0.

La forma más fácil sería convertir el ráster a netCDF y luego abrirlo y recorrer la cuadrícula. Hice casi lo mismo para un proyecto que involucra convertir rásteres en datos de características basados ​​en datos asignados a las celdas de ráster. Observé esto por años, y llegué a la conclusión de que recorrer los datos de la cuadrícula sería más fácil con netCDF.