Comment puis-je parcourir chaque cellule dans un raster continu?

Voir ce lien pour plus de détails.

Le problème:

Je veux parcourir un raster continu (un qui n'a pas de table attributaire), cellule par cellule, et obtenir la valeur de la cellule. Je veux prendre ces valeurs et exécuter des conditions dessus, en émulant les étapes d'algèbre de carte détaillées ci-dessous sans réellement utiliser la calculatrice raster.

Par demande de commentaires ci-dessous, j'ai ajouté des détails fournissant un contexte au problème et justifiant la nécessité de mettre en œuvre une méthode en tant que telle dans la section ci-dessous intitulée "L'analyse nécessaire:".

L'analyse proposée ci-dessous, tout en étant pertinente pour mon problème en fournissant un contexte, n'a pas besoin d'être implémentée dans une réponse. La portée de la question ne concerne que l'itération à travers un raster continu pour obtenir / définir les valeurs des cellules.

L'analyse avait besoin:

Si TOUTES les conditions suivantes sont remplies, donnez à la cellule de sortie une valeur de 1. Ne donnez à la cellule de sortie une valeur de 0 que si aucune des conditions n'est remplie.

Condition 1: si la valeur de la cellule est supérieure aux cellules supérieure et inférieure, donnez la valeur 1:

Con("raster" > FocalStatistics("raster", NbrIrregular("C:\filepath\kernel_file.txt"), "MAXIMUM"), 1, 0)

Où le fichier noyau ressemble à ceci:

3 3 
0 1 0
0 0 0
0 1 0

Condition 2: si la valeur de la cellule est supérieure aux cellules gauche et droite, donnez la valeur 1:

Con("raster" > FocalStatistics("raster", NbrIrregular("C:\filepath\kernel_file.txt"), "MAXIMUM"), 1, 0)

Où le fichier noyau ressemble à ceci:

3 3 
0 0 0
1 0 1
0 0 0  

Condition 3: si la valeur de la cellule est supérieure aux cellules de gauche et de droite, donnez la valeur 1:

Con("raster" > FocalStatistics("raster", NbrIrregular("C:\filepath\kernel_file.txt"), "MAXIMUM"), 1, 0)

Où le fichier noyau ressemble à ceci:

3 3 
1 0 0
0 0 0
0 0 1 

Condition 4: si la valeur de la cellule est supérieure aux cellules de gauche et de droite, donnez la valeur 1:

Con("raster" > FocalStatistics("raster", NbrIrregular("C:\filepath\kernel_file.txt"), "MAXIMUM"), 1, 0)

Où le fichier noyau ressemble à ceci:

3 3 
0 0 1
0 0 0
1 0 0 

Condition 5: si l' une des cellules adjacentes a une valeur ÉGALE à la cellule centrale, donnez au raster en sortie une valeur de 1 (en utilisant la variété focale avec deux calculs de voisinage les plus proches )

Pourquoi ne pas utiliser l'algèbre cartographique?

Il a été noté ci-dessous que mon problème pourrait être résolu en utilisant l'algèbre de la carte, mais comme vu ci-dessus, il s'agit d'un total de six calculs raster, plus un pour combiner tous les rasters créés ensemble. Il me semble qu'il est beaucoup plus efficace de procéder cellule par cellule et de faire toutes les comparaisons à la fois dans chaque cellule au lieu de parcourir sept fois individuellement et d'utiliser beaucoup de mémoire pour créer sept rasters.

Comment le problème doit-il être attaqué?

Le lien ci-dessus conseille d'utiliser l'interface IPixelBlock, mais il n'est pas clair dans la documentation ESRI si vous accédez réellement à une seule valeur de cellule via IPixelBlock, ou si vous accédez à plusieurs valeurs de cellule à partir de la taille de l'IPixelBlock que vous définissez. Une bonne réponse devrait suggérer une méthode pour accéder aux valeurs de cellule d'un raster continu et fournir une explication de la méthodologie derrière le code, si ce n'est apparemment évident.

En résumé:

Quelle est la meilleure méthode pour parcourir chaque cellule d'un raster CONTINU (qui n'a pas de table attributaire ) pour accéder à ses valeurs de cellule?

Une bonne réponse n'a pas besoin de mettre en œuvre les étapes d'analyse décrites ci-dessus, il suffit de fournir une méthodologie pour accéder aux valeurs de cellule d'un raster.

Je vois que cela a déjà été résolu par l'affiche originale (OP), mais je publierai une solution simple en python au cas où quelqu'un à l'avenir serait intéressé par différentes manières de résoudre ce problème. Je suis amateur de logiciels open source, voici donc une solution utilisant GDAL en python:

import gdal

#Set GeoTiff driver
driver = gdal.GetDriverByName("GTiff")
driver.Register()

#Open raster and read number of rows, columns, bands
dataset = gdal.Open(filepath)
cols = dataset.RasterXSize
rows = dataset.RasterYSize
allBands = dataset.RasterCount
band = dataset.GetRasterBand(1)

#Get array of raster cell values.  The two zeros tell the 
#iterator which cell to start on and the 'cols' and 'rows' 
#tell the iterator to iterate through all columns and all rows.
def get_raster_cells(band,cols,rows):
    return band.ReadAsArray(0,0,cols,rows)

Implémentez la fonction comme ceci:

#Bind array to a variable
rasterData = get_raster_cells(band,cols,rows)

#The array will look something like this if you print it
print rasterData
> [[ 1, 2, 3 ],
   [ 4, 5, 6 ],
   [ 7, 8, 9 ]]

Ensuite, parcourez vos données avec une boucle imbriquée:

for row in rasterData:
    for val in row:
        print val
> 1
  2
  3
  4...

Ou peut-être souhaitez-vous aplatir votre réseau 2D avec une compréhension de liste:

flat = [val for row in rasterData for val in row]

Quoi qu'il en soit, tout en parcourant les données cellule par cellule, il est possible de jeter des conditions dans votre boucle pour modifier / modifier les valeurs. Voir ce script que j'ai écrit pour différentes façons d'accéder aux données: https://github.com/azgs/hazards-viewer/blob/master/python/zonal_stats.py .

Mise à jour! La solution numpy:

import arcpy
import numpy as np

in_ras = path + "/rastername"

raster_Array = arcpy.RasterToNumPyArray(in_ras)
row_num = raster_Array.shape[0]
col_num = raster_Array.shape[1]
cell_count = row_num * row_num

row = 0
col = 0
temp_it = 0

while temp_it < cell_count:
    # Insert conditional statements
    if raster_Array[row, col] > 0:
        # Do something
        val = raster_Array[row, col]
        print val
    row+=1
    if col > col_num - 1:
        row = 0
        col+=1

Il est donc difficile de remettre le tableau fini sur le raster à l'aide d'arcpy. arcpy.NumPyArrayToRaster est écureuil et a tendance à redéfinir les étendues même si vous lui donnez vos coordonnées LL.

Je préfère enregistrer en tant que texte.

np.savetxt(path + "output.txt", output, fmt='%.10f', delimiter = " ")

J'utilise Python en 64 bits pour la vitesse - en ce moment, cela signifie que je ne peux pas nourrir numpy.savetxt un en-tête. Je dois donc ouvrir la sortie et ajouter l'en-tête ASCII que Arc veut avant de convertir ASCII en Raster

File_header = "NCOLS xxx" + '\n'+ "NROWS xxx" + '\n' + "XLLCORNER xxx"+'\n'+"YLLCORNER xxx"+'\n'+"CELLSIZE xxx"+'\n'+"NODATA_VALUE xxx"+'\n'

La version numpy exécute mon raster shift, mes multiplications et mes ajouts beaucoup plus rapidement (1000 itérations en 2 minutes) que la version arcpy (1000 itérations en 15 min)

ANCIENNE VERSION Je peux supprimer cela plus tard, je viens d'écrire un script similaire. J'ai essayé de convertir en points et d'utiliser le curseur de recherche. J'ai obtenu seulement 5000 itérations en 12 heures. J'ai donc cherché une autre voie.

Ma façon de procéder consiste à parcourir les coordonnées du centre de chaque cellule. Je commence dans le coin supérieur gauche et me déplace de droite à gauche. À la fin du rang, je descends d'un rang et recommence à gauche. J'ai un raster de 240 m avec 2603 colonnes et 2438 lignes donc un total de 6111844 cellules au total. J'utilise une variable itérateur et une boucle while. Voir ci-dessous

Quelques notes: 1 - vous devez connaître les coordonnées de l'étendue

2 - exécuter avec les coordonnées du point pour le centre de la cellule - déplacer de la moitié de la taille de la cellule à partir des valeurs d'étendue

3 - Mon script utilise la valeur de cellule pour extraire un raster spécifique à une valeur, puis décaler ce raster pour le centrer sur la cellule d'origine. Cela ajoute à un raster zéro pour étendre l'étendue avant de l'ajouter à un raster final. C'est juste un exemple. Vous pouvez mettre vos instructions conditionnelles ici (deuxième instruction if dans la boucle while).

4 - Ce script suppose que toutes les valeurs raster peuvent être converties en entiers. Cela signifie que vous devez d'abord vous débarrasser de l'absence de données. Con IsNull.

6 - Je ne suis toujours pas satisfait de cela et je travaille pour le supprimer complètement d'Arcpy. Je préfère lancer des tableaux numpy et faire le calcul là-bas, puis le ramener à Arc.

ULx = 959415 ## coordinates for the Upper Left of the entire raster 
ULy = 2044545
x = ULx ## I redefine these if I want to run over a smaller area
y = ULy
temp_it = 0

while temp_it < 6111844: # Total cell count in the data extent
        if x <= 1583895 and y >= 1459474: # Coordinates for the lower right corner of the raster
           # Get the Cell Value
           val_result = arcpy.GetCellValue_management(inraster, str(x)+" " +str(y), "1")
           val = int(val_result.getOutput(0))
        if val > 0: ## Here you could insert your conditional statements
            val_pdf = Raster(path + "pdf_"str(val))
            shift_x  =  ULx - x # This will be a negative value
            shift_y = ULy - y # This will be a positive value
            arcpy.Shift_management(val_pdf, path+ "val_pdf_shift", str(-shift_x), str(-shift_y))
            val_pdf_shift = Raster(path + "val_pdf_shift")
            val_pdf_sh_exp = CellStatistics([zeros, val_pdf_shift], "SUM", "DATA")
            distr_days = Plus(val_pdf_sh_exp, distr_days)
        if temp_it % 20000 == 0: # Just a print statement to tell me how it's going
                print "Iteration number " + str(temp_it) +" completed at " + str(time_it)
        x += 240 # shift x over one column
        if x > 1538295: # if your at the right hand side of a row
            y = y-240 # Shift y down a row
            x = 959415 # Shift x back to the first left hand column
        temp_it+=1

distr_days.save(path + "Final_distr_days")

Essayez d'utiliser IGridTable, ICursor, IRow. Cet extrait de code sert à mettre à jour les valeurs des cellules raster, mais il montre les bases de l'itération:

Comment puis-je ajouter un nouveau champ dans une table attributaire raster et le parcourir?

Public Sub CalculateArea(raster As IRaster, areaField As String)
    Dim bandCol As IRasterBandCollection
    Dim band As IRasterBand

    Set bandCol = raster
    Set band = bandCol.Item(0)

    Dim hasTable As Boolean
    band.hasTable hasTable
    If (hasTable = False) Then
        Exit Sub
    End If    

    If (AddVatField(raster, areaField, esriFieldTypeDouble, 38) = True) Then
        ' calculate cell size
        Dim rstProps As IRasterProps
        Set rstProps = raster

        Dim pnt As IPnt
        Set pnt = rstProps.MeanCellSize

        Dim cellSize As Double
        cellSize = (pnt.X + pnt.Y) / 2#

        ' get fields index
        Dim attTable As ITable
        Set attTable = band.AttributeTable

        Dim idxArea As Long, idxCount As Long
        idxArea = attTable.FindField(areaField)
        idxCount = attTable.FindField("COUNT")

        ' using update cursor
        Dim gridTableOp As IGridTableOp
        Set gridTableOp = New gridTableOp

        Dim cellCount As Long, cellArea As Double

        Dim updateCursor As ICursor, updateRow As IRow
        Set updateCursor = gridTableOp.Update(band.RasterDataset, Nothing, False)
        Set updateRow = updateCursor.NextRow()
        Do Until updateRow Is Nothing
            cellCount = CLng(updateRow.Value(idxCount))
            cellArea = cellCount * (cellSize * cellSize)

            updateRow.Value(idxArea) = cellArea
            updateCursor.updateRow updateRow

            Set updateRow = updateCursor.NextRow()
        Loop

    End If
End Sub

Une fois que vous parcourez la table, vous pouvez obtenir la valeur de ligne de champ spécifique en utilisant row.get_Value(yourfieldIndex). Si vous Google

arcobjects row.get_Value

vous devriez pouvoir obtenir de nombreux exemples le montrant.

J'espère que cela pourra aider.

Que diriez-vous d'une idée radicale, cela nécessiterait que vous programmiez en python ou ArcObjects.

1. Convertissez votre grille en points de caractéristiques.
2. Créez des champs XY et remplissez.
3. Chargez les points dans un dictionnaire où la clé est une chaîne de X, Y et l' élément est la valeur de la cellule.
4. Parcourez votre dictionnaire et déterminez pour chaque point les 8 XY de cellule environnants.
5. Récupérez-les dans votre dictionnaire et testez avec vos règles, dès que vous trouvez une valeur vraie, vous pouvez ignorer les autres tests.
6. Écrivez les résultats dans un autre dictionnaire, puis reconvertissez-les en grille en créant d'abord un point FeatureClass, puis convertissez les points en grille.

Une solution:

J'ai résolu cela plus tôt dans la journée. Le code est une adaptation de cette méthode . Le concept derrière cela n'a pas été terriblement difficile une fois que j'ai compris ce que les objets utilisés pour interfacer avec le raster font réellement. La méthode ci-dessous prend deux jeux de données d'entrée (inRasterDS et outRasterDS). Ils sont tous deux le même ensemble de données, je viens de faire une copie de inRasterDS et de le passer dans la méthode en tant que outRasterDS. De cette façon, ils ont tous deux la même étendue, la même référence spatiale, etc. La méthode lit les valeurs dans inRasterDS, cellule par cellule, et effectue des comparaisons avec les voisins les plus proches. Il utilise les résultats de ces comparaisons comme valeurs stockées dans outRasterDS.

Le processus:

J'ai utilisé IRasterCursor -> IPixelBlock -> SafeArray pour obtenir les valeurs en pixels et IRasterEdit pour en écrire de nouvelles dans le raster. Lorsque vous créez IPixelBlock, vous indiquez à la machine la taille et l'emplacement de la zone dans laquelle vous souhaitez lire / écrire. Si vous souhaitez uniquement modifier la moitié inférieure d'un raster, vous le définissez comme vos paramètres IPixelBlock. Si vous souhaitez boucler sur l'ensemble du raster, vous devez définir IPixelBlock égal à la taille de l'ensemble du raster. Je fais cela dans la méthode ci-dessous en passant la taille à IRasterCursor (pSize) puis en obtenant le PixelBlock du curseur raster.

L'autre clé est que vous devez utiliser SafeArray pour interfacer avec les valeurs de cette méthode. Vous obtenez IPixelBlock d'IRasterCursor, puis SafeArray d'IPixelBlock. Ensuite, vous lisez et écrivez dans SafeArray. Lorsque vous avez terminé la lecture / écriture sur SafeArray, réécrivez l'intégralité de votre SafeArray sur IPixelBlock, puis écrivez votre IPixelBlock sur IRasterCursor, puis utilisez enfin IRasterCursor pour définir l'emplacement de démarrage de l'écriture et IRasterEdit pour effectuer l'écriture elle-même. Cette dernière étape est l'endroit où vous modifiez réellement les valeurs de l'ensemble de données.

    public static void CreateBoundaryRaster(IRasterDataset2 inRasterDS, IRasterDataset2 outRasterDS)
    {
        try
        {
            //Create a raster. 
            IRaster2 inRaster = inRasterDS.CreateFullRaster() as IRaster2; //Create dataset from input raster
            IRaster2 outRaster = outRasterDS.CreateFullRaster() as IRaster2; //Create dataset from output raster
            IRasterProps pInRasterProps = (IRasterProps)inRaster;
            //Create a raster cursor with a pixel block size matching the extent of the input raster
            IPnt pSize = new DblPnt();
            pSize.SetCoords(pInRasterProps.Width, pInRasterProps.Height); //Give the size of the raster as a IPnt to pass to IRasterCursor
            IRasterCursor inrasterCursor = inRaster.CreateCursorEx(pSize); //Create IRasterCursor to parse input raster 
            IRasterCursor outRasterCursor = outRaster.CreateCursorEx(pSize); //Create IRasterCursor to parse output raster
            //Declare IRasterEdit, used to write the new values to raster
            IRasterEdit rasterEdit = outRaster as IRasterEdit;
            IRasterBandCollection inbands = inRasterDS as IRasterBandCollection;//set input raster as IRasterBandCollection
            IRasterBandCollection outbands = outRasterDS as IRasterBandCollection;//set output raster as IRasterBandCollection
            IPixelBlock3 inpixelblock3 = null; //declare input raster IPixelBlock
            IPixelBlock3 outpixelblock3 = null; //declare output raster IPixelBlock
            long blockwidth = 0; //store # of columns of raster
            long blockheight = 0; //store # of rows of raster

            //create system array for input/output raster. System array is used to interface with values directly. It is a grid that overlays your IPixelBlock which in turn overlays your raster.
            System.Array inpixels; 
            System.Array outpixels; 
            IPnt tlc = null; //set the top left corner

            // define the 3x3 neighborhood objects
            object center;
            object topleft;
            object topmiddle;
            object topright;
            object middleleft;
            object middleright;
            object bottomleft;
            object bottommiddle;
            object bottomright;

            long bandCount = outbands.Count; //use for multiple bands (only one in this case)

            do
            {

                inpixelblock3 = inrasterCursor.PixelBlock as IPixelBlock3; //get the pixel block from raster cursor
                outpixelblock3 = outRasterCursor.PixelBlock as IPixelBlock3;
                blockwidth = inpixelblock3.Width; //set the # of columns in raster
                blockheight = inpixelblock3.Height; //set the # of rows in raster
                outpixelblock3.Mask(255); //set any NoData values

                for (int k = 0; k < bandCount; k++) //for every band in raster (will always be 1 in this case)
                {
                    //Get the pixel array.
                    inpixels = (System.Array)inpixelblock3.get_PixelData(k); //store the raster values in a System Array to read
                    outpixels = (System.Array)outpixelblock3.get_PixelData(k); //store the raster values in a System Array to write
                    for (long i = 1; i < blockwidth - 1; i++) //for every column (except outside columns)
                    {
                        for (long j = 1; j < blockheight - 1; j++) //for every row (except outside rows)
                        {
                            //Get the pixel values of center cell and  neighboring cells

                            center = inpixels.GetValue(i, j);

                            topleft = inpixels.GetValue(i - 1, j + 1);
                            topmiddle = inpixels.GetValue(i, j + 1);
                            topright = inpixels.GetValue(i + 1, j + 1);
                            middleleft = inpixels.GetValue(i - 1, j);
                            middleright = inpixels.GetValue(i + 1, j);
                            bottomleft = inpixels.GetValue(i - 1, j - 1);
                            bottommiddle = inpixels.GetValue(i, j - 1);
                            bottomright = inpixels.GetValue(i - 1, j - 1);


                            //compare center cell value with middle left cell and middle right cell in a 3x3 grid. If true, give output raster value of 1
                            if ((Convert.ToDouble(center) >= Convert.ToDouble(middleleft)) && (Convert.ToDouble(center) >= Convert.ToDouble(middleright)))
                            {
                                outpixels.SetValue(1, i, j);
                            }


                            //compare center cell value with top middle and bottom middle cell in a 3x3 grid. If true, give output raster value of 1
                            else if ((Convert.ToDouble(center) >= Convert.ToDouble(topmiddle)) && (Convert.ToDouble(center) >= Convert.ToDouble(bottommiddle)))
                            {
                                outpixels.SetValue(1, i, j);
                            }

                            //if neither conditions are true, give raster value of 0
                            else
                            {

                                outpixels.SetValue(0, i, j);
                            }
                        }
                    }
                    //Write the pixel array to the pixel block.
                    outpixelblock3.set_PixelData(k, outpixels);
                }
                //Finally, write the pixel block back to the raster.
                tlc = outRasterCursor.TopLeft;
                rasterEdit.Write(tlc, (IPixelBlock)outpixelblock3);
            }
            while (inrasterCursor.Next() == true && outRasterCursor.Next() == true);
            System.Runtime.InteropServices.Marshal.ReleaseComObject(rasterEdit);


        }
        catch (Exception ex)
        {
            MessageBox.Show(ex.Message);
        }

    }

Les données raster AFAIK peuvent être lues de trois manières:

• par cellule (inefficace);
• par image (assez efficace);
• par blocs (le moyen le plus efficace).

Sans réinventer la roue, je vous propose de lire ces diapositives éclairantes de Chris Garrard.

Ainsi, la méthode la plus efficace consiste à lire les données par bloc, mais cela entraînerait une perte de données dans la correspondance des pixels situés sur les limites du bloc lors de l'application du filtre. Une autre solution sûre devrait donc consister à lire l'image entière en une seule fois et à utiliser l'approche numpy.

Du côté informatique, je devrais plutôt utiliser gdalfilter.py et implicitement l'approche VRT KernelFilteredSource afin d'appliquer les filtres nécessaires et, surtout, d'éviter les calculs lourds.