Obtenir une image raster en tant que tableau en Python avec ArcGIS Desktop?

10

Lorsque j'ai commencé à travailler avec Python et ArcGIS 9.3, j'ai supposé qu'il y aurait un moyen simple d'obtenir une image raster dans un tableau Python afin de pouvoir la manipuler avant de la stocker à nouveau en tant qu'autre image raster. Cependant, je n'arrive pas à trouver comment procéder.

Si c'est possible alors, comment?

arcpy raster arcgis-9.3 robintw
la source

Réponses:

6

Je ne pense pas que cela soit possible avec ArcGIS <= 9.3.1

J'utilise l'open source API GDAL pour des tâches telles que cela.

fmark
la source
Génial! J'ai utilisé les programmes utilitaires GDAL dans le passé, mais je n'ai jamais pensé à les utiliser pour ce faire.
robintw
3
Je suis d'accord, le module gdal Python vous permet de lire facilement un raster et de vider les données dans un tableau Numpy. Chris Garrard a un cours sur l'utilisation d'OpenSource Python dans les SIG, il couvre ce sujet. Vous pouvez le trouver sur: gis.usu.edu/~chrisg/python/2008
DavidF
6

fmark a déjà répondu à la question, mais voici un exemple de code OSGEO Python que j'ai écrit pour lire un raster (tif) dans un tableau NumPy, reclasser les données puis les écrire dans un nouveau fichier tif. Vous pouvez lire et écrire n'importe quel format pris en charge par 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
DavidF
la source
1

Vous pouvez enregistrer votre raster en tant que grille ascii ESRI et lire / manipuler ce fichier avec numpy.

Cela fournit quelques points de départ: http://sites.google.com/site/davidpfinlayson2/esriasciigridformat

Mais attention - il semble que le format de la grille ascii ne suive pas toujours les spécifications, donc les lire correctement à chaque fois peut être un défi.

snorris
la source
1

Je ne suis pas sûr que vous puissiez manipuler le raster pixel par pixel, mais vous pouvez utiliser les objets de géotraitement en conjonction avec l'API python.

Vous pouvez utiliser n'importe quelle boîte à outils pour ce type de manipulation. Un exemple de script serait:

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

Voici un suivi de votre question . Toujours pas possible. Pas sûr de la version 10.0.

George Silva
la source
Merci - c'est très utile. Cependant, idéalement, j'aimerais pouvoir parcourir la matrice raster en faisant diverses choses. J'aurais pensé qu'il y aurait eu un moyen dans ArcGIS de le faire, mais peut-être pas!
robintw
robintw, pour ce que j'ai regardé dans la référence, il n'y a aucun moyen d'obtenir un pixel spécifique d'un raster. Je ne sais pas si dans ArcPy (disponible à partir de la v10), vous pouvez récupérer ces cellules individuelles, car elles ont étendu l'API python avec beaucoup de nouvelles fonctionnalités.
George Silva
0

Le moyen le plus simple serait de convertir le raster en netCDF, puis de l'ouvrir et de parcourir la grille. J'ai fait à peu près la même chose pour un projet impliquant la transformation de rasters en données d'entités basées sur les données affectées aux cellules raster. J'ai regardé cela pendant des siècles et je suis arrivé à la conclusion que parcourir les données de la grille serait plus facile avec netCDF.

Poilu
la source