Otrzymujesz obraz rastrowy jako tablicę w Pythonie dzięki ArcGIS Desktop?

10

Zaczynając pracę z Python i ArcGIS 9.3, założyłem, że istnieje prosty sposób na przeniesienie obrazu rastrowego do tablicy Python, aby móc nim manipulować przed zapisaniem go jako kolejnego obrazu rastrowego. Jednak nie mogę się dowiedzieć, jak to zrobić.

Jeśli to możliwe, to jak?

arcpy raster arcgis-9.3

— robintw
źródło

6

Nie sądzę, że jest to możliwe dzięki ArcGIS <= 9.3.1

Korzystam z GDAL API typu open source do takich zadań.

— fmark
źródło

Świetny! W przeszłości korzystałem z programów narzędziowych GDAL, ale nigdy nie myślałem o ich użyciu.

— robintw

3

Zgadzam się, moduł gdal Python pozwala na łatwy odczyt rastra i zrzucenie danych do tablicy Numpy. Chris Garrard ma kurs używania OpenSource Python w GIS, który obejmuje ten temat. Można go znaleźć na stronie: gis.usu.edu/~chrisg/python/2008

— DavidF

10

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

ArcGIS 10 ma możliwość zapisywania i odczytywania tablic numPy.

— George Silva
źródło

Dzięki! Niestety, obecnie korzystam z ArcGIS 9.3 - czy wiesz, czy można to zrobić w 9.3?

— robintw

6

fmark już odpowiedział na pytanie, ale oto przykładowy kod Python OSGEO, który napisałem, aby odczytać raster (tif) do tablicy NumPy, przeklasyfikować dane, a następnie zapisać je w nowym pliku tif. Możesz czytać i pisać w dowolnym formacie obsługiwanym przez 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
źródło

2

Uzyskujesz dostęp do ArcObjects z Python? omawia integrację obiektów arkobowych z pythonem.

Być może kod w tym przykładzie mógłby zostać dostosowany, aby można go było wywoływać z Pythona.

Nie jestem pewien, czy istnieje sposób na przekazanie tablicy byref z powrotem do Pythona. Jeśli tak, to warto spróbować IPixelBlock.PixelDatabyRef .

— Kirk Kuykendall
źródło

1

Możesz zapisać swój raster jako siatkę ESCI ascii i czytać / manipulować tym plikiem za pomocą numpy.

Zapewnia to kilka punktów wyjścia: http://sites.google.com/site/davidpfinlayson2/esriasciigridformat

Ale uważaj - wydaje się, że format siatki ascii nie zawsze jest zgodny ze specyfikacją, więc prawidłowe czytanie za każdym razem może być wyzwaniem.

— Snorris
źródło

1

Nie jestem pewien, czy możesz manipulować rastrem piksel po pikselu, ale możesz używać obiektów geoprzetwarzania w połączeniu z API Pythona.

Do tego rodzaju manipulacji można użyć dowolnego zestawu narzędzi. Przykładowy skrypt to:

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

Oto odpowiedź na twoje pytanie . Wciąż niemożliwe. Nie jestem pewien w wersji 10.0.

— George Silva
źródło

Dzięki - to bardzo pomocne. Jednak idealnie chciałbym móc iterować po tablicy rastrowej, wykonując różne czynności. Wydawało mi się, że w ArcGIS byłby na to sposób, ale może nie!

— robintw

robintw, dla tego, czego szukałem w referencji, nie ma sposobu na uzyskanie określonego piksela rastra. Nie jestem pewien, czy w ArcPy (dostępny od wersji 10) możesz pobrać te pojedyncze komórki, ponieważ rozszerzyły one API Pythona o wiele nowych funkcji.

— George Silva,

0

Najprostszym sposobem byłoby przekonwertowanie rastra na netCDF, a następnie otwarcie go i przejście przez siatkę. Zrobiłem to samo w przypadku projektu polegającego na przekształceniu rastrów w dane funkcji oparte na danych przypisanych do komórek rastrowych. Patrzyłem na to od wieków i doszedłem do wniosku, że chodzenie po siatce danych byłoby łatwiejsze od netCDF.

— Włochaty
źródło