Jak iterować po każdej komórce w ciągłym rastrze?

Zobacz ten link, aby uzyskać więcej informacji.

Problem:

Chcę przejść przez ciągły raster (taki, który nie ma tabeli atrybutów), komórka po komórce i uzyskać wartość komórki. Chcę wziąć te wartości i uruchomić na nich warunki warunkowe, emulując kroki algebry mapy wyszczególnione poniżej, bez korzystania z kalkulatora rastrowego.

Na prośbę o komentarze poniżej dodałem szczegółowe informacje przedstawiające tło problemu i uzasadniające potrzebę wdrożenia metody jako takiej w poniższej sekcji zatytułowanej „Wymagana analiza:”.

Przedstawiona poniżej analiza, choć jest odpowiednia dla mojego problemu, ponieważ stanowi tło, nie musi być realizowana w odpowiedzi. Zakres pytania dotyczy tylko iteracji przez ciągły raster, aby uzyskać / ustawić wartości komórek.

Potrzebna analiza:

Jeśli DOWOLNY z poniższych warunków jest spełniony, nadaj komórce wyjściowej wartość 1. Daj komórce wyjściowej wartość 0, tylko jeśli żaden z warunków nie zostanie spełniony.

Warunek 1: Jeśli wartość komórki jest większa niż górna i dolna komórka, podaj wartość 1:

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

Gdzie plik jądra wygląda następująco:

3 3 
0 1 0
0 0 0
0 1 0

Warunek 2: Jeśli wartość komórki jest większa niż lewa i prawa komórka, podaj wartość 1:

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

Gdzie plik jądra wygląda następująco:

3 3 
0 0 0
1 0 1
0 0 0  

Warunek 3: Jeśli wartość komórki jest większa niż górne i dolne komórki, podaj wartość 1:

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

Gdzie plik jądra wygląda następująco:

3 3 
1 0 0
0 0 0
0 0 1 

Warunek 4: Jeśli wartość komórki jest większa niż dolnej i lewej górnej komórki, podaj wartość 1:

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

Gdzie plik jądra wygląda następująco:

3 3 
0 0 1
0 0 0
1 0 0 

Warunek 5: Jeśli którakolwiek z sąsiednich komórek ma wartość RÓWNOMOCNĄ w stosunku do komórki środkowej, nadaj wyjściowemu rastrowi wartość 1 ( używając odmiany ogniskowej z dwoma obliczeniami najbliższego sąsiedztwa )

Dlaczego nie skorzystać z algebry map?

Zauważono poniżej, że mój problem można rozwiązać za pomocą algebry mapy, ale jak widać powyżej, jest to w sumie sześć obliczeń rastrowych plus jedno, aby połączyć wszystkie rastry utworzone razem. Wydaje mi się, że o wiele bardziej wydajne jest przechodzenie między komórkami i wykonywanie wszystkich porównań jednocześnie w każdej komórce zamiast zapętlania każdego z nich osobno siedem razy i wykorzystywania całkiem sporej ilości pamięci do utworzenia siedmiu rastrów.

Jak należy zaatakować problem?

Powyższy link zaleca użycie interfejsu IPixelBlock, jednak z dokumentacji ESRI nie jest jasne, czy faktycznie uzyskuje się dostęp do samej wartości pojedynczej komórki za pośrednictwem IPixelBlock, czy też uzyskuje się dostęp do wielu wartości komórek z ustawionego rozmiaru IPixelBlock. Dobra odpowiedź powinna sugerować metodę dostępu do wartości komórek ciągłego rastra i wyjaśnić metodologię kodu, jeśli nie jest to oczywiste.

W podsumowaniu:

Jaka jest najlepsza metoda na przejście przez każdą komórkę w CIĄGŁYM rastrze (który nie ma tabeli atrybutów ), aby uzyskać dostęp do jej wartości komórek?

Dobra odpowiedź nie musi implementować opisanych powyżej etapów analizy, wystarczy jedynie zapewnić metodologię dostępu do wartości komórek rastra.

Widzę, że problem został już rozwiązany przez Original Poster (OP), ale opublikuję proste rozwiązanie w pythonie na wypadek, gdyby ktoś w przyszłości był zainteresowany różnymi sposobami rozwiązania tego problemu. Jestem częściowo zwolennikiem oprogramowania typu open source, więc oto rozwiązanie wykorzystujące GDAL w pythonie:

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)

Zaimplementuj funkcję w następujący sposób:

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

Następnie iteruj dane za pomocą zagnieżdżonej pętli:

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

A może chcesz spłaszczyć tablicę 2-D ze zrozumieniem listy:

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

W każdym razie, podczas iteracji danych po komórce, możliwe jest wrzucenie do pętli niektórych warunków warunkowych w celu zmiany / edycji wartości. Zobacz ten skrypt który napisałem dla różnych sposobów, aby uzyskać dostęp do danych: https://github.com/azgs/hazards-viewer/blob/master/python/zonal_stats.py .

Aktualizacja! Rozwiązanie 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

Tak więc odzyskanie gotowej tablicy z powrotem do rastra przy użyciu arcpy jest kłopotliwe. arcpy.NumPyArrayToRaster jest wiewiórczy i ma tendencję do redefiniowania zakresu, nawet jeśli podajesz mu swoje współrzędne LL.

Wolę zapisać jako tekst.

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

Używam Pythona jako 64-bitowej prędkości - od tej chwili oznacza to, że nie mogę nakarmić numpy.savetxt nagłówka. Więc muszę otworzyć wyjście i dodać nagłówek ASCII, który Arc chce przed konwersją ASCII na Raster

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

Wersja numpy uruchamia mój raster shift, mnożenie i dodawanie znacznie szybciej (1000 iteracji w 2 minuty) niż wersja arcpy (1000 iteracji w 15 minut)

STARE WERSJE Mogę to później usunąć właśnie napisałem podobny skrypt. Próbowałem przekonwertować na punkty i użyć kursora wyszukiwania. Mam tylko 5000 iteracji w ciągu 12 godzin. Szukałem więc innej drogi.

Moim sposobem na to jest iteracja przez współrzędne centrum komórki każdej komórki. Zaczynam w lewym górnym rogu i przechodzę od prawej do lewej. Na końcu rzędu przechodzę w dół o rząd i zaczynam od lewej. Mam raster 240 mz 2603 kolumnami i 2438 rzędami, więc łącznie 6111844 komórek. Używam zmiennej iteratora i pętli while. Patrz poniżej

Kilka uwag: 1 - musisz znać współrzędne zasięgu

2 - uruchom ze współrzędnymi punktowymi dla centrum komórki - przesuń o 1/2 rozmiaru komórki od wartości zasięgu

3 - Mój skrypt używa wartości komórki do pobrania rastra zależnego od wartości, a następnie przesunięcie tego rastra na środek oryginalnej komórki. Dodaje się to do zera rastrowego, aby zwiększyć zakres przed dodaniem do końcowego rastra. To tylko przykład. Tutaj możesz umieścić swoje instrukcje warunkowe (drugie polecenie if w pętli while).

4 - Ten skrypt zakłada, że ​​wszystkie wartości rastrowe mogą być rzutowane jako liczby całkowite. Oznacza to, że najpierw musisz pozbyć się żadnych danych. Con IsNull.

6 - Nadal nie jestem z tego zadowolony i pracuję nad tym, aby całkowicie usunąć to z arkady. Wolałbym rzucić się jak dziwaczne tablice i zrobić tam matematykę, a potem sprowadzić ją z powrotem do 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")

Spróbuj użyć IGridTable, ICursor, IRow. Ten fragment kodu służy do aktualizowania wartości komórek rastrowych, jednak pokazuje podstawy iteracji:

Jak mogę dodać nowe pole do tabeli atrybutów rastrowych i przejść przez nią?

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

Po przejściu przez tabelę można uzyskać określoną wartość wiersza pola za pomocą row.get_Value(yourfieldIndex). Jeśli ty Google

arcobjects row.get_Value

powinieneś być w stanie uzyskać wiele przykładów tego pokazujących.

Mam nadzieję, że to pomaga.

Co powiesz na to jako radykalny pomysł, wymagałoby to programowania w Pythonie lub ArcObjects.

1. Przekształć siatkę w punktowe klasy obiektów.
2. Utwórz pola XY i wypełnij.
3. Załaduj punkty do słownika, w którym klucz to ciąg X, Y, a element to wartość komórki.
4. Przejrzyj słownik i dla każdego punktu oblicz 8 otaczających komórek XY.
5. Pobierz je ze słownika i przetestuj zgodnie z regułami, gdy tylko znajdziesz prawdziwą wartość, możesz pominąć pozostałe testy.
6. Zapisz wyniki do innego słownika, a następnie przekonwertuj z powrotem na siatkę, najpierw tworząc punkt FeatureClass, a następnie konwertuj punkty na siatkę.

Rozwiązanie:

Rozwiązałem to dzisiaj dzisiaj. Kod jest adaptacją tej metody . Koncepcja tego nie była strasznie trudna, kiedy zorientowałem się, co faktycznie robią obiekty używane do interakcji z rastrem. Poniższa metoda pobiera dwa zestawy danych wejściowych (inRasterDS i outRasterDS). Oba są tym samym zestawem danych, właśnie utworzyłem kopię inRasterDS i przekazałem ją do metody jako outRasterDS. W ten sposób oba mają ten sam zasięg, odniesienie przestrzenne itp. Metoda odczytuje wartości z inRasterDS, komórka po komórce i dokonuje na nich porównań najbliższego sąsiada. Wykorzystuje wyniki tych porównań jako zapisane wartości w outRasterDS.

Proces:

Użyłem IRasterCursor -> IPixelBlock -> SafeArray, aby uzyskać wartości pikseli, a IRasterEdit, aby zapisać nowe wartości w rastrze. Podczas tworzenia IPixelBlock mówisz maszynie o wielkości i lokalizacji obszaru, w którym chcesz czytać / zapisywać. Jeśli chcesz edytować tylko dolną połowę rastra, ustaw to jako parametry IPixelBlock. Jeśli chcesz zapętlić cały raster, musisz ustawić IPixelBlock równy rozmiarowi całego rastra. Robię to w poniższej metodzie, przekazując rozmiar do IRasterCursor (pSize), a następnie pobierając PixelBlock z kursora rastrowego.

Drugi klucz polega na tym, że musisz użyć SafeArray do połączenia z wartościami w tej metodzie. Otrzymasz IPixelBlock od IRasterCursor, a następnie SafeArray od IPixelBlock. Następnie czytasz i piszesz w SafeArray. Po zakończeniu odczytu / zapisu w SafeArray, napisz całą SafeArray z powrotem do IPixelBlock, a następnie napisz swój IPixelBlock do IRasterCursor, a następnie użyj IRasterCursor, aby ustawić lokalizację, w której chcesz rozpocząć zapis, a IRasterEdit, aby sam napisał. Ten ostatni krok to miejsce, w którym faktycznie edytujesz wartości zestawu danych.

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

    }

Dane rastrowe AFAIK można odczytać na trzy sposoby:

• według komórki (nieefektywne);
• według obrazu (dość wydajny);
• według bloków (najbardziej wydajny sposób).

Bez wymyślania koła proponuję przeczytać te oświecające slajdy Chrisa Garrarda.

Zatem najbardziej wydajną metodą jest odczyt danych po bloku, jednak spowodowałoby to utratę danych w korespondencji pikseli znajdujących się ponad granicami bloku podczas stosowania filtra. Bezpieczny alternatywny sposób powinien polegać na jednoczesnym odczytaniu całego obrazu i zastosowaniu metody numpy.

Zamiast tego po stronie obliczeniowej powinienem użyć gdalfilter.py i domyślnie metody VRT KernelFilteredSource, aby zastosować potrzebne filtry, a przede wszystkim uniknąć ciężkich obliczeń.