Como obter coordenadas XY e o valor da célula de cada pixel em uma varredura usando Python?

Sou realmente novo em Python e gostaria de saber se existe um método rápido para obter valores de célula de um pixel raster por pixel e as coordenadas (mapear coordenadas XY do centro de cada pixel) usando Python no ArcGIS 10?

Para descrever isso ainda mais, preciso obter o mapa X, o mapa Y e o valor da célula do primeiro pixel, atribuir esses três valores a três variáveis ​​e repetir esta etapa para o restante dos outros pixels (percorrer toda a varredura).

Eu acho que preciso descrever minha pergunta mais. O problema é que eu preciso obter a localização XY de um pixel da primeira varredura e obter os valores das células de várias outras rasters correspondentes a essa localização XY. Esse processo deve percorrer todos os pixels da primeira varredura sem criar nenhum arquivo shapefile de ponto intermediário, pois consumirá muito tempo, pois tenho que lidar com uma varredura com quase 8 bilhões de pixels. Além disso, preciso fazer isso usando Python no ArcGIS 10.

@ JamesS: Muito obrigado pela sua sugestão. Sim, isso funcionaria para uma varredura, mas eu preciso coletar os valores das células para várias outras rasters também. O problema é que, depois de obter as coordenadas X e Y do primeiro pixel da primeira varredura, preciso obter o valor da célula da segunda varredura correspondente àquele local X, Y da primeira varredura, depois à terceira varredura e assim por diante. Então, acho que, ao percorrer o primeiro raster, obter a localização X e Y de um pixel e obter os valores das células da outra raster correspondente a esse local devem ser feitos simultaneamente, mas não tenho certeza. Isso pode ser feito convertendo o primeiro raster em um shapefile de ponto e executando Extrair multivalores para a função de ponto no ArcGIS 10, mas eu '

@ hmfly: Obrigado, Sim, este método (RastertoNumpyarray) funcionará se eu conseguir obter a coordenada de um valor conhecido de linha e coluna da matriz.

@ whuber: Eu não quero fazer nenhum cálculo, tudo o que preciso fazer é escrever coordenadas XY e valores de célula em um arquivo de texto e isso é tudo

Seguindo a ideia do @ Dango, criei e testei (em pequenos rasters com a mesma extensão e tamanho de célula) o seguinte código:

import arcpy, numpy

inRaster = r"C:\tmp\RastersArray.gdb\InRaster"
inRaster2 = r"C:\tmp\RastersArray.gdb\InRaster2"

##Get properties of the input raster
inRasterDesc = arcpy.Describe(inRaster)

#coordinates of the lower left corner
rasXmin = inRasterDesc.Extent.Xmin
rasYmin = inRasterDesc.Extent.Ymin

# Cell size, raster size
rasMeanCellHeight = inRasterDesc.MeanCellHeight
rasMeanCellWidth = inRasterDesc.MeanCellWidth
rasHeight = inRasterDesc.Height
rasWidth = inRasterDesc.Width

##Calculate coordinates basing on raster properties
#create numpy array of coordinates of cell centroids
def rasCentrX(rasHeight, rasWidth):
    coordX = rasXmin + (0.5*rasMeanCellWidth + rasWidth)
    return coordX
inRasterCoordX = numpy.fromfunction(rasCentrX, (rasHeight,rasWidth)) #numpy array of X coord

def rasCentrY(rasHeight, rasWidth):
    coordY = rasYmin + (0.5*rasMeanCellHeight + rasHeight)
    return coordY
inRasterCoordY = numpy.fromfunction(rasCentrY, (rasHeight,rasWidth)) #numpy array of Y coord

#combine arrays of coordinates (although array for Y is before X, dstack produces [X, Y] pairs)
inRasterCoordinates = numpy.dstack((inRasterCoordY,inRasterCoordX))


##Raster conversion to NumPy Array
#create NumPy array from input rasters 
inRasterArrayTopLeft = arcpy.RasterToNumPyArray(inRaster)
inRasterArrayTopLeft2 = arcpy.RasterToNumPyArray(inRaster2)

#flip array upside down - then lower left corner cells has the same index as cells in coordinates array
inRasterArray = numpy.flipud(inRasterArrayTopLeft)
inRasterArray2 = numpy.flipud(inRasterArrayTopLeft2)


# combine coordinates and value
inRasterFullArray = numpy.dstack((inRasterCoordinates, inRasterArray.T))

#add values from second raster
rasterValuesArray = numpy.dstack((inRasterFullArray, inRasterArray2.T))

Com base no código @hmfly, você pode ter acesso aos valores desejados:

(height, width, dim )=rasterValuesArray.shape
for row in range(0,height):
    for col in range(0,width):
        #now you have access to single array of values for one cell location

Infelizmente, existe um 'mas' - o código é adequado para matrizes NumPy que podem ser manipuladas pela memória do sistema. Para o meu sistema (8 GB), a maior matriz era de cerca de 9000.9000.

Como minha experiência não me permite fornecer mais ajuda, você pode considerar algumas sugestões sobre como lidar com matrizes grandes: /programming/1053928/python-numpy-very-large-matrices

arcpy.RasterToNumPyArrayO método permite especificar o subconjunto da varredura convertida em matriz NumPy ( página de ajuda do ArcGIS10 ) o que pode ser útil ao agrupar grandes conjuntos de dados em submatrizes.

Se você deseja apenas obter os valores de pixel (linha, coluna), pode escrever um script arcpy como este:

import arcpy
raster = arcpy.Raster("yourfilepath")
array = arcpy.RasterToNumPyArray(raster)
(height, width)=array.shape
for row in range(0,height):
    for col in range(0,width):
        print str(row)+","+str(col)+":"+str(array.item(row,col))

Mas, se você deseja obter a coordenada do pixel, o NumPyArray não pode ajudá-lo. Você pode converter a varredura em ponto pela ferramenta RasterToPoint e, em seguida, obter a coordenada arquivada pelo Shape.

O método mais simples para gerar coordenadas e valores de células em um arquivo de texto no ArcGIS 10 é a função de amostra , sem necessidade de código e, especialmente, sem necessidade de repetir cada célula. Na calculadora raster ArcGIS <= 9.3x , costumava ser tão simples quantooutfile.csv = sample(someraster) produzir um arquivo de texto com todos os valores e coordenadas de células (não nulos) (no formato z, x, y). No ArcGIS 10, parece que o argumento "in_location_data" agora é obrigatório, então você precisa usar a sintaxe Sample(someraster, someraster, outcsvfile).

Editar: Você também pode especificar vários rasters: Sample([someraster, anotherraster, etc], someraster, outcsvfile) . Se isso funcionaria em 8 bilhões de células, não tenho idéia ...

Edit: Note, eu não testei isso no ArcGIS 10, mas usei a função de amostra por anos em <= 9.3 (e Workstation).

Edit: Agora testei no ArcGIS 10 e ele não será gerado em um arquivo de texto. A ferramenta altera a extensão do arquivo para ".dbf" automaticamente. No entanto ... o código python a seguir funciona como instruções de álgebra de mapa SOMA e MOMA ainda são suportadas no ArcGIS 10:

import arcgisscripting
gp=arcgisscripting.create()
gp.multioutputmapalgebra(r'%s=sample(%s)' % (outputcsv,inputraster))

Uma maneira de fazer isso seria usar a ferramenta Raster_To_Point seguida pela ferramenta Add_XY_Coordinates . Você terminará com um shapefile em que cada linha na tabela de atributos representa um pixel da sua varredura com colunas para X_Coord , Y_Coord e Cell_Value . Em seguida, você pode percorrer esta tabela usando um cursor (ou exportá-la para algo como o Excel, se preferir).

Se você tiver apenas uma varredura para processar, provavelmente não vale a pena criar scripts - basta usar as ferramentas do ArcToolbox. Se você precisar fazer isso para muitos rasters, tente algo como isto:

[ Nota: Eu não tenho o ArcGIS 10 e não estou familiarizado com o ArcPy, portanto, este é apenas um esboço muito aproximado. Não foi testado e quase certamente precisará de ajustes para que funcione.]

import arcpy, os
from arcpy import env

# User input
ras_fold = r'path/to/my/data'           # The folder containing the rasters
out_fold = r'path/to/output/shapefiles' # The folder in which to create the shapefiles

# Set the workspace
env.workspace = ras_fold

# Get a list of raster datasets in the raster folder
raster_list = arcpy.ListRasters("*", "All")

# Loop over the rasters
for raster in raster_list:
    # Get the name of the raster dataset without the file extension
    dataset_name = os.path.splitext(raster)[0]

    # Build a path for the output shapefile
    shp_path = os.path.join(out_fold, '%s.shp' % dataset_name)

    # Convert the raster to a point shapefile
    arcpy.RasterToPoint_conversion(raster, shp_path, "VALUE")

    # Add columns to the shapefile containing the X and Y co-ordinates
    arcpy.AddXY_management(shp_path)

Você pode fazer um loop sobre as tabelas de atributos shapefile usando um Cursor de Pesquisa ou (possivelmente mais simples) usando o dbfpy . Isso permitirá que você leia os dados da sua varredura (agora armazenados em uma tabela .dbf do shapefile) em variáveis ​​python.

from dbfpy import dbf

# Path to shapefile .dbf
dbf_path = r'path\to\my\dbf_file.dbf'

# Open the dbf file
db = dbf.Dbf(dbf_path)

# Loop over the records
for rec in db:
    cell_no = rec['POINTID'] # Numbered from top left, running left to right along each row
    cell_x = rec['POINT_X']
    cell_y = rec['POINT_Y']
    cell_val = rec['GRID_CODE']

    # Print values
    print cell_no, cell_x, cell_y, cell_val

Talvez você possa criar um arquivo mundial para a varredura, encobri-la em uma matriz numpy. então, se você fizer um loop sobre a matriz, obterá os valores das células e, se você atualizar constantemente x, y do arquivo mundial, também terá as coordenadas para cada valor de célula. espero que seja útil.

O código de Marcin funcionou bem, exceto que um problema nas funções rasCentrX e rasCentrY estava fazendo com que as coordenadas de saída aparecessem em uma resolução diferente (como observou Grazia). Minha correção foi mudar

coordX = rasXmin + (0.5*rasMeanCellWidth + rasWidth)

para

coordX = rasXmin + ((0.5 + rasWidth) * rasMeanCellWidth)

e

  coordY = rasYmin + (0.5*rasMeanCellHeight + rasHeight)

para

  coordY = rasYmin + ((0.5 + rasHeight) * rasMeanCellHeight)

Usei o código para converter uma grade ESRI em um arquivo CSV. Isso foi obtido removendo a referência ao inRaster2 e usando um csv.writer para gerar as coordenadas e os valores:

out = csv.writer(open(outputCSV,"wb"), delimiter=',', quoting=csv.QUOTE_NONNUMERIC)
out.writerow(['X','Y','Value'])
(height, width, dim )=inRasterFullArray.shape
for row in range(0,height):
    for col in range(0,width):
        out.writerow(inRasterFullArray[row,col])

Também não achei a transposição necessária

inRasterFullArray = numpy.dstack((inRasterCoordinates, inRasterArray.T))

tão convertido que para

inRasterFullArray = numpy.dstack((inRasterCoordinates, inRasterArray))

Feio, mas altamente eficaz:

1. Crie um novo recurso de ponto com 4 pontos fora dos cantos da varredura em questão. Verifique se o mesmo sistema de coordenadas da varredura em questão.
2. Adicione os campos duplos 'xcor' e 'ycor'
3. Calcular geometria para obter coordenadas para esses campos
4. Analista espacial-> Interpolação-> Tendência -> Regressão linear
5. Configurações do ambiente: ajuste o tamanho da célula e a varredura para a mesma que a varredura em questão
6. Executar separadamente para 'xcor' e 'ycor'
7. Saem os avaliadores com coordenadas como valores de célula, usados ​​como entrada para scripts.

Uma solução simples usando pacotes python de código aberto:

import fiona
import rasterio
from pprint import pprint


def raster_point_coords(raster, points):

    # initialize dict to hold data
    pt_data = {}

    with fiona.open(points, 'r') as src:
        for feature in src:
            # create dict entry for each feature
            pt_data[feature['id']] = feature

    with rasterio.open(raster, 'r') as src:
        # read raster into numpy array
        arr = src.read()
        # rasterio always reads into 3d array, this is 2d, so reshape
        arr = arr.reshape(arr.shape[1], arr.shape[2])
        # get affine, i.e. data needed to work between 'image' and 'raster' coords
        a = src.affine

    for key, val in pt_data.items():
        # get coordinates
        x, y = val['geometry']['coordinates'][0], val['geometry']['coordinates'][1]
        # use affine to convert to row, column
        col, row = ~a * (x, y)
        # remember numpy array is indexed array[row, column] ie. y, x
        val['raster_value'] = arr[int(row), int(col)]

    pprint(pt_data) 

if __name__ == '__main__':
    # my Landsat raster
    ras = '/data01/images/sandbox/LT05_040028_B1.tif'
    # my shapefile with two points which overlap raster area
    pts = '/data01/images/sandbox/points.shp'
    # call function
    raster_point_coords(ras, pts)

Fiona é útil, pois você pode abrir um shapefile, percorrer os recursos e (como eu) anexá-los a um dictobjeto. De fato, a featureprópria Fiona é como umdict , por isso é fácil acessar propriedades. Se meus pontos tivessem algum atributo, eles apareceriam neste ditado junto com as coordenadas, id, etc.

O Rasterio é útil porque é fácil de ler no raster como uma matriz numpy, um tipo de dados leve e rápido. Também temos acesso a uma dictdas propriedades de varredura, incluindo o affine, que é todos os dados que precisamos para converter as coordenadas x, y de varredura em coordenadas de linha, col. Veja a excelente explicação de @ perrygeo aqui .

Terminamos com um pt_datatipo dictque possui dados para cada ponto e o extraído raster_value. Poderíamos reescrever facilmente o shapefile com os dados extraídos, se quiséssemos.