Mendapatkan gambar raster sebagai array di Python dengan ArcGIS Desktop?

10

Ketika mulai bekerja dengan Python dan ArcGIS 9.3, saya berasumsi akan ada cara sederhana untuk mendapatkan gambar raster ke dalam array Python sehingga saya dapat memanipulasinya sebelum menyimpannya kembali sebagai gambar raster lainnya. Namun, sepertinya saya tidak tahu bagaimana melakukan ini.

Jika memungkinkan, lalu bagaimana caranya?

arcpy raster arcgis-9.3 Robintw
sumber

Jawaban:

6

Saya rasa ini tidak mungkin dengan ArcGIS <= 9.3.1

Saya menggunakan GDAL API open source untuk tugas-tugas seperti ini.

fmark
sumber
Bagus! Saya telah menggunakan program utilitas GDAL di masa lalu, tetapi tidak pernah berpikir untuk menggunakannya untuk melakukan ini.
robintw
3
Saya setuju, modul Python gdal memungkinkan Anda untuk dengan mudah membaca raster dan membuang data ke array Numpy. Chris Garrard memiliki kursus tentang penggunaan OpenSource Python di GIS, ini mencakup topik ini. Anda dapat menemukannya di: gis.usu.edu/~chrisg/python/2008
DavidF
6

fmark sudah menjawab pertanyaan itu, tetapi di sini ada beberapa contoh kode Python OSGEO yang saya tulis untuk membaca raster (tif) ke dalam array NumPy, mengklasifikasikan ulang data dan kemudian menuliskannya ke file tif baru. Anda dapat membaca dan menulis format apa pun yang didukung 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
sumber
2

Mengakses ArcObjects dari Python? membahas integrasi objek busur dengan python.

Mungkin kode dalam sampel ini dapat diadaptasi sehingga dapat dipanggil dari python.

Saya tidak yakin apakah ada cara untuk melewatkan array byref kembali ke python. Jika ada, maka IPixelBlock.PixelDatabyRef akan layak untuk dicoba.

Kirk Kuykendall
sumber
1

Anda dapat menyimpan raster Anda sebagai kisi ascii ESRI dan membaca / memanipulasi file itu dengan numpy.

Ini memberikan beberapa poin awal: http://sites.google.com/site/davidpfinlayson2/esriasciigridformat

Tapi hati-hati - sepertinya format grid ascii tidak selalu mengikuti spec, jadi membacanya dengan benar setiap waktu bisa menjadi tantangan.

Snorris
sumber
1

Saya tidak yakin Anda dapat memanipulasi raster pixel demi pixel, tetapi Anda dapat menggunakan objek geoprocessing bersamaan dengan python API.

Anda dapat menggunakan kotak alat apa pun untuk manipulasi semacam itu. Contoh skrip adalah:

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

Ini adalah tindak lanjut dari pertanyaan Anda . Masih tidak mungkin. Tidak yakin pada versi 10.0.

George Silva
sumber
Terima kasih - itu sangat membantu. Namun, idealnya saya ingin dapat beralih di array raster melakukan berbagai hal untuk itu. Saya akan berpikir akan ada cara di ArcGIS untuk melakukan ini, tetapi mungkin tidak!
robintw
robintw, untuk apa yang telah saya lihat dalam referensi, tidak ada cara untuk mendapatkan pixel raster tertentu. Saya tidak yakin apakah di ArcPy (tersedia dari v10) Anda dapat mengambil sel-sel individual ini, karena mereka memperpanjang API python dengan banyak fungsi baru.
George Silva
0

Cara termudah adalah dengan mengubah raster ke netCDF dan kemudian membukanya dan berjalan melalui grid. Saya melakukan banyak hal yang sama untuk sebuah proyek yang melibatkan mengubah raster menjadi data fitur berdasarkan data yang ditugaskan ke sel raster. Saya melihat ini sejak lama, dan sampai pada kesimpulan berjalan data grid akan lebih mudah dari netCDF.

Berbulu
sumber