Получение растрового изображения в виде массива в Python с ArcGIS Desktop?

10

Когда я начинал работать с Python и ArcGIS 9.3, я предполагал, что будет простой способ вставить растровое изображение в массив Python, чтобы я мог манипулировать им, прежде чем сохранить его как другое растровое изображение. Тем не менее, я не могу понять, как это сделать.

Если это возможно, то как?

arcpy raster arcgis-9.3 robintw
источник

Ответы:

6

Я не думаю, что это возможно с ArcGIS <= 9.3.1

Я использую GDAL API с открытым исходным кодом для таких задач, как эта.

fmark
источник
Большой! В прошлом я использовал служебные программы GDAL, но никогда не думал об их использовании для этого.
robintw
3
Я согласен, модуль Python gdal позволяет легко читать растр и выгружать данные в массив Numpy. Крис Гаррард имеет курс по использованию OpenSource Python в ГИС, он охватывает эту тему. Вы можете найти его по адресу: gis.usu.edu/~chrisg/python/2008
DavidF
6

fmark уже ответил на вопрос, но вот несколько примеров кода Python для OSGEO, который я написал, чтобы прочитать растр (tif) в массив NumPy, переклассифицировать данные и записать их в новый файл tif. Вы можете читать и писать любой формат, поддерживаемый 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
источник
2

Доступ к ArcObjects из Python? обсуждается интеграция arcobjects с python.

Возможно, код в этом примере можно адаптировать так, чтобы его можно было вызывать из python.

Я не уверен, есть ли способ передать массив byref обратно в python. Если есть, то стоит попробовать IPixelBlock.PixelDatabyRef .

Кирк Куйкендалл
источник
1

Вы можете сохранить свой растр как сетку ESRI ascii и читать / манипулировать этим файлом с помощью numpy.

Это обеспечивает некоторые отправные точки: http://sites.google.com/site/davidpfinlayson2/esriasciigridformat

Но будьте осторожны - кажется, что формат сетки ascii не всегда соответствует спецификациям, поэтому правильное их чтение каждый раз может быть проблемой.

snorris
источник
1

Я не уверен, что вы можете манипулировать растровым пикселем за пикселем, но вы можете использовать объекты геообработки в сочетании с Python API.

Вы можете использовать любой набор инструментов для такого рода манипуляций. Пример сценария будет:

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

Вот продолжение вашего вопроса . Все еще невозможно. Не уверен в версии 10.0.

Джордж Сильва
источник
Спасибо - это очень полезно. Однако в идеале я хотел бы иметь возможность выполнять итерацию по растровому массиву, выполняя различные действия с ним. Я бы подумал, что в ArcGIS был бы способ сделать это, но, возможно, нет!
robintw
robintw, для того, что я посмотрел в ссылке, нет никакого способа получить определенный пиксель растра. Я не уверен, что в ArcPy (доступна с v10) вы можете получить эти отдельные ячейки, так как они расширили Python API с множеством новых функций.
Джордж Сильва
0

Самый простой способ - преобразовать растр в netCDF, а затем открыть его и пройтись по сетке. Я сделал то же самое для проекта, включающего превращение растров в данные объектов, основанные на данных, присвоенных растровым ячейкам. Я смотрел на это целую вечность и пришел к выводу, что ходить по сетке данных будет проще от netCDF.

волосатый
источник