Comment obtenir les coordonnées XY et la valeur de cellule de chaque pixel dans un raster en utilisant Python?

Je suis vraiment nouveau sur Python et j'aimerais savoir s'il existe une méthode rapide pour obtenir les valeurs de cellule d'un raster pixel par pixel et les coordonnées (mapper les coordonnées XY du centre de chaque pixel) en utilisant Python dans ArcGIS 10?

Pour décrire cela plus en détail, je dois obtenir la carte X, la carte Y et la valeur de cellule du premier pixel et affecter ces trois valeurs à trois variables et répéter cette étape pour le reste des autres pixels (boucle à travers le raster entier).

Je pense que je dois décrire ma question davantage. Le problème est que je dois obtenir l'emplacement XY d'un pixel du premier raster et obtenir les valeurs de cellule de plusieurs autres rasters correspondant à cet emplacement XY. Ce processus devrait parcourir tous les pixels du premier raster sans créer de fichier de formes de points intermédiaires car cela va vraiment prendre beaucoup de temps car je dois gérer un raster avec près de 8 milliards de pixels. De plus, je dois le faire en utilisant Python dans ArcGIS 10.

@JamesS: Merci beaucoup pour votre suggestion. Oui, cela fonctionnerait pour un raster, mais je dois également collecter les valeurs des cellules pour plusieurs autres rasters. Le problème est qu'après avoir obtenu les coordonnées X et Y du premier pixel du premier raster, j'ai besoin d'obtenir la valeur de cellule du deuxième raster correspondant à cet emplacement X, Y du premier raster, puis du troisième raster et ainsi de suite. Donc, je pense que lorsque vous parcourez le premier raster, obtenir l'emplacement X et Y d'un pixel et obtenir les valeurs de cellule de l'autre raster correspondant à cet emplacement doit être fait simultanément, mais je ne suis pas sûr. Cela peut être fait en convertissant le premier raster en un fichier de formes de points et en effectuant la fonction d'extraction de plusieurs valeurs en points dans ArcGIS 10, mais je '

@hmfly: Merci, oui, cette méthode (RastertoNumpyarray) fonctionnera si je peux obtenir les coordonnées d'une valeur de ligne et de colonne connue du tableau.

@whuber: Je ne veux pas effectuer de calculs, tout ce que je dois faire est d'écrire les coordonnées XY et les valeurs des cellules dans un fichier texte et c'est tout

En suivant l'idée de @ Dango, j'ai créé et testé (sur de petits rasters avec la même étendue et la même taille de cellule) le code suivant:

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

Basé sur le code @hmfly, vous pouvez avoir accès aux valeurs souhaitées:

(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

Malheureusement, il y a un «mais» - le code convient aux tableaux NumPy qui peuvent être gérés par la mémoire système. Pour mon système (8 Go), la plus grande baie était d'environ 9 000 9 000.

Comme mon expérience ne me permet pas de fournir plus d'aide, vous pouvez envisager quelques suggestions sur le traitement avec de grands tableaux: /programming/1053928/python-numpy-very-large-matrices

arcpy.RasterToNumPyArraypermet de spécifier le sous-ensemble de raster converti en tableau NumPy ( page d'aide d'ArcGIS10 ) ce qui peut être utile lors du découpage d'un grand ensemble de données en sous-matrices.

Si vous voulez simplement faire passer les valeurs des pixels (ligne, colonne), vous pouvez écrire un script arcpy comme ceci:

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

Mais, si vous souhaitez obtenir les coordonnées du pixel, NumPyArray ne peut pas vous aider. Vous pouvez convertir le raster en point par l'outil RasterToPoint, puis vous pouvez obtenir les coordonnées par forme déposée.

La méthode la plus simple pour générer des coordonnées et des valeurs de cellule dans un fichier texte dans ArcGIS 10 est l' exemple de fonction , pas besoin de code et surtout pas besoin de boucler sur chaque cellule. Dans la calculatrice raster ArcGIS <= 9.3x, elle était aussi simple que celle outfile.csv = sample(someraster)qui produirait un fichier texte de toutes les valeurs et coordonnées (non nulles) des cellules (au format z, x, y). Dans ArcGIS 10, il semble que l'argument "in_location_data" soit désormais obligatoire, vous devez donc utiliser la syntaxe Sample(someraster, someraster, outcsvfile).

Edit: Vous pouvez également spécifier plusieurs rasters: Sample([someraster, anotherraster, etc], someraster, outcsvfile). Que cela fonctionne sur 8 milliards de cellules, je n'ai aucune idée ...

Edit: Remarque, je n'ai pas testé cela dans ArcGIS 10, mais j'ai utilisé l'exemple de fonction pendant des années dans <= 9.3 (et Workstation).

Edit: j'ai maintenant testé dans ArcGIS 10 et il ne sortira pas dans un fichier texte. L'outil modifie automatiquement l'extension de fichier en ".dbf". Cependant ... le code python suivant fonctionne car les instructions d'algèbre de carte SOMA et MOMA sont toujours prises en charge dans ArcGIS 10:

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

Une façon de procéder consiste à utiliser l' outil Raster_To_Point suivi de l' outil Add_XY_Coordinates . Vous vous retrouverez avec un fichier de formes où chaque ligne de la table attributaire représente un pixel de votre raster avec des colonnes pour X_Coord , Y_Coord et Cell_Value . Vous pouvez ensuite parcourir ce tableau à l'aide d'un curseur (ou l'exporter vers quelque chose comme Excel si vous préférez).

Si vous n'avez qu'un seul raster à traiter, cela ne vaut probablement pas la peine de faire un script - utilisez simplement les outils d'ArcToolbox. Si vous devez effectuer cette opération pour de nombreux rasters, vous pouvez essayer quelque chose comme ceci:

[ Remarque: je n'ai pas ArcGIS 10 et je ne connais pas ArcPy, il ne s'agit donc que d'un aperçu très approximatif. Il n'a pas été testé et devra presque certainement être peaufiné pour le faire fonctionner.]

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)

Vous pouvez ensuite parcourir les tables d'attributs de fichiers de formes à l'aide d'un curseur de recherche ou (peut-être plus simple) à l'aide de dbfpy . Cela vous permettra de lire les données de votre raster (maintenant stockées dans une table shapefile .dbf) dans des variables 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

Vous pouvez peut-être créer un fichier universel pour le raster, convertir le raster en un tableau numpy. puis si vous bouclez sur le tableau, vous obtiendrez les valeurs des cellules et si vous mettez à jour progressivement les x, y à partir du fichier mondial, vous aurez également les coordonnées de chaque valeur de cellule. j'espère que c'est utile.

Le code de Marcin a bien fonctionné, sauf un problème dans les fonctions rasCentrX et rasCentrY qui faisait apparaître les coordonnées de sortie à une résolution différente (comme l'a observé Grazia). Ma solution était de changer

coordX = rasXmin + (0.5*rasMeanCellWidth + rasWidth)

à

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

et

  coordY = rasYmin + (0.5*rasMeanCellHeight + rasHeight)

à

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

J'ai utilisé le code pour convertir une grille ESRI en un fichier CSV. Ceci a été réalisé en supprimant la référence à inRaster2, puis en utilisant un csv.writer pour sortir les coordonnées et les valeurs:

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])

Je n'ai pas non plus trouvé que la transposition était nécessaire dans

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

tellement converti en

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

Moche mais très efficace:

1. Créez une nouvelle entité ponctuelle avec 4 points en dehors des coins du raster en question. Assurez-vous dans le même système de coordonnées que le raster en question.
2. Ajouter des champs doubles «xcor» et «ycor»
3. Calculer la géométrie pour obtenir les coordonnées de ces champs
4. Analyste spatial-> Interpolation-> Tendance -> Régression linéaire
5. Paramètres d'environnement: alignement du raster et de la taille de cellule sur le raster en question
6. Exécuter séparément pour 'xcor' et 'ycor'
7. Il en résulte des évaluateurs avec des coordonnées comme valeurs de cellule, à utiliser comme entrée pour les scripts.

Une solution simple utilisant des packages open source python:

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 est pratique car vous pouvez ouvrir un fichier de formes, parcourir les fonctionnalités et (comme je l'ai) les ajouter à un dictobjet. En effet, le Fiona featurelui-même est comme un dictaussi, il est donc facile d'accéder aux propriétés. Si mes points avaient des attributs, ils apparaîtront dans ce dict avec les coordonnées, l'id, etc.

Rasterio est pratique car il est facile à lire dans le raster comme un tableau numpy, un type de données léger et rapide. Nous avons également accès à une dictdes propriétés du raster, y compris le affine, qui est toutes les données dont nous avons besoin pour convertir les coordonnées x, y du raster en ligne de tableau, coordonnées col. Voir l'excellente explication de @ perrygeo ici .

On se retrouve avec un pt_datade type dictqui a des données pour chaque point et l'extrait raster_value. Nous pourrions facilement réécrire le fichier de formes avec les données extraites si nous le voulions.