Rasterbild als Array in Python mit ArcGIS Desktop abrufen?

10

Als ich anfing, mit Python und ArcGIS 9.3 zu arbeiten, ging ich davon aus, dass es eine einfache Möglichkeit gibt, ein Rasterbild in ein Python-Array zu übertragen, damit ich es bearbeiten kann, bevor ich es wieder als anderes Rasterbild speichere. Ich kann jedoch nicht herausfinden, wie das geht.

Wenn es möglich ist, wie dann?

arcpy raster arcgis-9.3 robintw
quelle

Antworten:

6

Ich denke nicht, dass dies mit ArcGIS <= 9.3.1 möglich ist

Ich verwende die Open-Source- GDAL-API für solche Aufgaben.

fmark
quelle
Groß! Ich habe in der Vergangenheit die GDAL-Hilfsprogramme verwendet, aber nie daran gedacht, sie dazu zu verwenden.
Robintw
3
Ich stimme zu, mit dem gdal Python-Modul können Sie einfach ein Raster lesen und die Daten in einem Numpy-Array speichern. Chris Garrard hat einen Kurs zur Verwendung von OpenSource Python in GIS, der dieses Thema behandelt. Sie finden es unter: gis.usu.edu/~chrisg/python/2008
DavidF
6

fmark hat die Frage bereits beantwortet, aber hier ist ein Beispiel für einen OSGEO-Python-Code, den ich geschrieben habe, um ein Raster (tif) in ein NumPy-Array zu lesen, die Daten neu zu klassifizieren und sie dann in eine neue tif-Datei zu schreiben. Sie können jedes von gdal unterstützte Format lesen und schreiben.

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

Zugriff auf ArcObjects über Python? diskutiert die Integration von Arcobjects mit Python.

Möglicherweise könnte der Code in diesem Beispiel so angepasst werden, dass er von Python aus aufgerufen werden kann.

Ich bin nicht sicher, ob es eine Möglichkeit gibt, ein Array byref an Python zurückzugeben. Wenn ja , dann wäre IPixelBlock.PixelDatabyRef einen Versuch wert.

Kirk Kuykendall
quelle
1

Sie können Ihr Raster als ESRI-ASCII-Raster speichern und diese Datei mit numpy lesen / bearbeiten.

Dies bietet einige Ausgangspunkte: http://sites.google.com/site/davidpfinlayson2/esriasciigridformat

Aber aufgepasst - es scheint, dass das ASCII-Rasterformat nicht immer den Spezifikationen entspricht, daher kann es eine Herausforderung sein, sie jedes Mal richtig zu lesen.

Snorris
quelle
1

Ich bin nicht sicher, ob Sie das Raster Pixel für Pixel bearbeiten können, aber Sie können die Geoverarbeitungsobjekte in Verbindung mit der Python-API verwenden.

Sie können jede Toolbox für diese Art der Manipulation verwenden. Ein Beispielskript wäre:

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

Hier ist eine Antwort auf Ihre Frage . Immer noch nicht möglich. Nicht sicher in Version 10.0.

George Silva
quelle
Danke - das ist sehr hilfreich. Idealerweise möchte ich jedoch in der Lage sein, über das Raster-Array zu iterieren und dabei verschiedene Dinge zu tun. Ich hätte gedacht, dass es in ArcGIS eine Möglichkeit gegeben hätte, dies zu tun, aber vielleicht auch nicht!
Robintw
robintw, für das, was ich in der Referenz gesehen habe, gibt es keine Möglichkeit, ein bestimmtes Pixel eines Rasters zu erhalten. Ich bin mir nicht sicher, ob Sie in ArcPy (verfügbar ab Version 10) diese einzelnen Zellen abrufen können, da sie die Python-API um viele neue Funktionen erweitert haben.
George Silva
0

Am einfachsten wäre es, das Raster in netCDF zu konvertieren, es dann zu öffnen und durch das Raster zu gehen. Ich habe das Gleiche für ein Projekt getan, bei dem Raster in Feature-Daten umgewandelt wurden, die auf Daten basieren, die den Rasterzellen zugewiesen wurden. Ich habe mir das ewig angesehen und bin zu dem Schluss gekommen, dass das Gehen der Rasterdaten von netCDF einfacher wäre.

Behaart
quelle