How do I iterate through every cell in a continuous raster?

Question

See this link for more details.

The Problem:

I want to loop through a continuous raster (one that has no attribute table), cell by cell, and get the value of the cell. I want to take those values and run conditionals on them, emulating the map algebra steps detailed below without actually using the raster calculator.

Per request of comments below, I have added details providing background to the problem and justifying the need to implement a method as such in the section below called "The analysis needed:".

The analysis proposed below, while being relevant to my problem by providing background,does not need to be implemented in an answer. The scope of the question only pertains to iterating through a continuous raster to get/set the cell values.

The analysis needed:

If ANY of the following conditions are fulfilled, give output cell a value of 1. Only give output cell a value of 0 if none of the conditions are fulfilled.

Condition 1: If cell value is greater than top and bottom cells, give value of 1:

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

Where kernel file looks like this:

3 3 
0 1 0
0 0 0
0 1 0

Condition 2: If cell value is greater than left and right cells, give value of 1:

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

Where kernel file looks like this:

3 3 
0 0 0
1 0 1
0 0 0  

Condition 3: If cell value is greater than topleft and bottomright cells, give value of 1:

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

Where kernel file looks like this:

3 3 
1 0 0
0 0 0
0 0 1 

Condition 4: If cell value is greater than bottomleft and topright cells, give value of 1:

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

Where kernel file looks like this:

3 3 
0 0 1
0 0 0
1 0 0 

Condition 5: If any one of the adjacent cells has a value EQUAL to the center cell, give the output raster a value of 1 (using focal variety with two nearest neighborhood calculations)

Why not use map algebra?

It has been noted below that my problem could be solved using map algebra but as seen above this is a grand total of six raster calculations, plus one to combine all of the rasters created together. It seems to me that it is much more efficient to go cell-by-cell and do all of the comparisons at once in each cell instead of looping through each individually seven times and utilizing quite a bit of memory to create seven rasters.

How should the problem be attacked?

The link above advises to use IPixelBlock interface, however it is unclear from ESRI documentation whether you are actually accessing a single cell value itself through IPixelBlock, or if you are accessing multiple cell values from the size of the IPixelBlock you set. A good answer should suggest a method for accessing the cell values of a continuous raster and provide an explanation of the methodology behind the code, if not apparently obvious.

In summary:

What is the best method to loop through every cell in a CONTINUOUS raster (which has no attribute table) to access its cell values?

A good answer need not implement the analysis steps described above, it needs only to provide a methodology to access cell values of a raster.

Answer 1

I see this has already been solved by the Original Poster (OP), but I'll post a simple solution in python just in case anyone in the future is interested in different ways to solve this problem. I'm partial to open source software, so here's a solution using GDAL in 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)

Implement the function like this:

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

Then, iterate through your data with a nested loop:

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

Or maybe you want to flatten your 2-D array with a list comprehension:

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

Anyways, while iterating through the data on a cell-by-cell basis its possible to throw some conditionals into your loop to change/edit values. See this script I wrote for different ways to access the data: https://github.com/azgs/hazards-viewer/blob/master/python/zonal_stats.py.

Answer 2

Update! The numpy solution:

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

So, getting the finished array back to raster using arcpy is troublesome. arcpy.NumPyArrayToRaster is squirrelly and tends to redefine extents even if you feed it your LL coordinates.

I prefer save as text.

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

I am running Python as 64 bit for speed - as of right now this means I can't feed numpy.savetxt a header. So I have to open the output and add the ASCII header that Arc wants before converting ASCII to Raster

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

The numpy version runs my shift raster, multiplications, and addition much faster (1000 iterations in 2 minutes) than the arcpy version (1000 iterations in 15 min)

OLD VERSION I may delete this later I just wrote a similar script. I tried converting to points and using the search cursor. I got only 5000 iterations in 12 hours. So, I looked for another way.

My way of doing this is to iterate through the cell center coordinates of each cell. I start in the upper left corner and move right to left. At the end of the row I move down a row and start again at the left. I have a 240 m raster with 2603 columns and 2438 rows so a total of 6111844 total cells. I use an iterator variable and a while loop. See below

A few notes: 1 - you need to know the coordinates of the extent

2 - run with point coordinates for cell center - move in 1/2 the cell size from the extent values

3 - My script is using the cell value to pull a value specific raster, then shift this raster to center on the original cell. This adds to a zero raster to expand the extent before adding into a final raster. This is just an example. You can put your conditional statements in here (second if statement within the while loop).

4 - This script assumes all of the raster values can be cast as integers. This means you need to get rid of the no data first. Con IsNull.

6 - I'm still not happy with this and I am working to take this out of arcpy entirely. I would rather cast as numpy arrays and do the math over there then bring it back to 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")

Answer 3

How about this as a radical idea, it would require you to program in python or ArcObjects.

1. Convert your grid to a point featureclasss.
2. Create XY fields and populate.
3. Load the points into a dictionary where key is a string of X,Y and item is cell value..
4. Step through your dictionary and for each point work out the 8 surrounding cell XYs.
5. Retrieve these from your dictionary and test with your rules, as soon as you find a value that is true you can skip the rest of the tests.
6. Write the results to another dictionary and then convert back to a grid by first creating a point FeatureClass and then convert points to a grid.

Answer 4

Try using IGridTable, ICursor, IRow. This code snippet is for updating raster cell values, however it shows the basics of iterating:

How can I add a new field in a raster attribute table and loop through it?

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

Once you are traversing the table you can get the specific field row value by using row.get_Value(yourfieldIndex). If you Google

arcobjects row.get_Value

you should be able to get plenty of examples showing this.

Hope that helps.

Answer 5

A solution:

I solved this earlier today. The code is an adaptation of this method. The concept behind this wasn't terribly difficult once I figured out what the objects used to interface with the raster actually do. The below method takes two input datasets (inRasterDS and outRasterDS). They are both the same dataset, I just made a copy of inRasterDS and passed it into the method as outRasterDS. This way they both have the same extent, spatial reference, etc. The method reads the values from inRasterDS, cell by cell, and does nearest neighbor comparisons on them. It uses the results of those comparisons as stored values in outRasterDS.

The Process:

I used IRasterCursor -> IPixelBlock -> SafeArray to get at the pixel values and IRasterEdit to write new ones to the raster. When you create IPixelBlock, you are telling the machine the size and location of the area that you want to read/write to. If you only want to edit the bottom half of a raster, you set that as your IPixelBlock parameters. If you want to loop over the entire raster, you have to set IPixelBlock equal to the size of the entire raster. I do this in the method below by passing the size to IRasterCursor (pSize) then getting the PixelBlock from the raster cursor.

The other key is that you have to use SafeArray to interface with the values in this method. You get IPixelBlock from IRasterCursor, then SafeArray from IPixelBlock. Then you read and write to SafeArray. When you finish read/write to SafeArray, write your entire SafeArray back to IPixelBlock, then write your IPixelBlock to IRasterCursor, then finally use IRasterCursor to set the location to start write and IRasterEdit to do the write itself. This final step is where you actually edit the values of the dataset.

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

    }

Answer 6

AFAIK raster data can be read in three ways:

• by cell (inefficient);
• by image (quite efficient);
• by blocks (the most efficient way).

Without reinventing the wheel, I suggest to read these enlightening slides of Chris Garrard.

So the most efficient method is read data by block, however this would cause a data loss in correspondence of pixels located over the block boundaries while applying the filter. So a safe alternative way should consist into reading the entire image at once and using the numpy approach.

On the computational side, instead, I should use gdalfilter.py and implicitly the VRT KernelFilteredSource approach in order to apply the needed filters and, above all, avoid heavy calculations.