Finding adjacent neighbors for Cellular Automata

[LindaVasil]

Hi all,

I'm trying to write a simple cellular automata code on a reclassified land use point shapefile.  It is reclassified by likelihood of land use conversion based on the number of like adjacent neighbors.  I would like to the code to run cell by cell and do the following:
Find adjacent neighbors
Get Grid_Code value
Find out if there is 4 of more cells with same land use grid code
If there is and this value is greater than the center cell, then the center cell converts.

I've written the following code, but it gets caught in select by attribute command.  The error that comes back is the following:
"C:\Python26\ArcGIS10.0\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 325, in RunScript
    exec codeObject in __main__.__dict__
  File "C:\Users\Flash\Documents\School\CA\CellularAutomata.py", line 32, in <module>
    arcpy.SelectLayerByAttribute_management(shp, "New_Selection", sql)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\management.py", line 4259, in SelectLayerByAttribute
    raise e
ExecuteError: Failed to execute. Parameters are not valid.
The value cannot be a feature class
ERROR 000840: The value is not a Raster Layer.
ERROR 000840: The value is not a Mosaic Layer.
Failed to execute (SelectLayerByAttribute).

Here is the code itself: 

#Import arc geoprocessing module
import arcpy, sys

#Set variable to the point shapefile
shp = sys.argv[1]

#Add LRG_Value to the point shapefile
arcpy.AddField_management(shp, "LRG_Value", "Long")

#Place an update cursor in the point shapefile table
#move to the first record in the table
cur = arcpy.UpdateCursor(shp)
row = cur.next()

#Loop through each record in the table
while row:
    #Extract Grid Code and Object ID
    GCP = row.getValue("Grid_Code")       #GCP = Grid Code of Processing Cell
    PIP = row.getValue("PointID")            #PIP = Point ID of Processing Cell
   
    #Put ObjectID in sequential order
    sql = "PointID" + " = " + str(PIP)
   
    #Select the sequential order of Object ID
    arcpy.SelectLayerByAttribute_management(shp, "New_Selection", sql)  #Where the mess up is!
   
    #Select only the Object ID that are surrounding the cursor cell
    Adj_Neighbors = arcpy.SelectLayerbyLocation_management(shp, "Boundary_Touches", shp, " ", "New Selection")
   
    for codes in Adj_Neighbors:
        GCA = row.getValue("Grid_Code")                  #GCA = Grid Codes from adjacent cells           
   
        #Put Grid codes in a list
        GCA_List = list(GCA)
   
        #Sort Grid Codes from smallest to largest
        srt_GCA_List = sorted(GCA_List)
   
        #select the largest number in the list (last number in sequence)
        LR_Value = srt_GCA_List[-1]
   
        #count the number of times the largest occurs
        tie_LR = GCA_List.count(LR_Value)
   
        if tie_LR >= 4 and LR_Value > GCA:  #is this syntax correct...seems a little messy to me.
            row.setValue("LRG_Value", LR_Value)
        #else:
        #    LR_Value = srt_GCA_List[-2]
        #keep going until the end of the list -- I will write this later
       
        #Update the row values
        cur.updateRow(row)

    #Delect temporary Variables
    del GCP, OIP, sql, result, GCA, GCA_List, srt_GCA_List, LR_Value, tie_LR

    #move to the next row
    row = cur.next()

    #delete row and cursor object, arcpy, module
    del cur, row, arcpy


I would appreciate any help on that line, or any other lines that are obviously wrong.  I'm pretty new to the code writing.

Thanks!
Linda

[LindaVasil]

Thanks very much...that is exactly what I had to do.  My basin had 250,000 points in it.  I had no other choice but to code up as a matrix in python, not only for the time to run, but I would also run out of memory.  Now with python, it runs 85 simulations in 2:20 seconds.  I would suggest to anyone doing something similar is to export your FID and Grid Codes to a text file and have python read it in as a dictionary.  And do all the rest of through python as well.

I wish I would have come back and read your post sooner...it would have saved me some time and frustration.

Just a comment, but depending on how many points you have (I am guessing these points were derived from a raster), this sort of vector based analysis could take an extremely long time. For example, if your original grid was only 1000 x 1000 cells, your code would need to loop 1,000,000 times. Assuming it takes even just 1 second to do your SelectbyLocation and SelectByAttributes (most likely more than 1 second!), your process will take about 12 days to run.

This analysis would be much better handled using a Python dictionary object (a numpy array is another option, but one I don't have any experience with). Dictionaries are an in memory data structure meaning that once you load your data from disk, the processing is, well, light speed. The basic idea is that you load your land use data into a dictionary by using the column and row order as the keys:

dict = {}
dict[1,1] = 4 #the cell existing at column = 1, row = 1 has a landuse code of 4
dict[2,1] = 4 #the cell existing at column = 2, row = 1 has a landuse code of 4
dict[1,2] = 1 #the cell existing at column = 1, row = 2 has a landuse code of 1
dict[2,2] = 2 #the cell existing at column = 2, row = 2 has a landuse code of 2

Anyway, look into using Python dictionaries or numpy arrays to do this sort of thing... It will be immensely faster and cooler. This sort or "neighbor notation" capability was actually built in functionality in good ole' Workstation ArcInfo Grid. However, ESRI chose not to port neighbor notation forward to ArcGIS, and so we are stuck devising our own data structures and algorithms to deal with ESRI apathy.:mad:

Hmmm...

[ChrisSnyder]

'Python Dictionary' is my answer to almost everything these days - It's an extremely powerful data structure for doing iterative-type analyses... as long as all your data fits.

[curtvprice]

Thanks very much...that is exactly what I had to do. 

1. Get Grid_Code value


2. Find out if there is 4 of more cells with same land use grid code


3. If there is and this value is greater than the center cell, then the center cell converts

Linda, can you share your code here? I'm curious how you loaded the info into the dict and how you worked with it once it was there. Neat that you got pretty nice performance with a generic python dictionary.

One thing to think about in this kind of analysis is whether the cell changes are "destructive", ie, when you move to the next cell to the right, and do the analysis again looking at the land use codes: does the analysis for the new location use the changed value, or the original?

Good old Workstation grid neighborhood notation does not do that, it always looks back at the original raster.

[ChrisSnyder]

I'll second that request. If you are willing to share, I'd also like to see that code too.

[RaphaelR]

if i understood all the posts correctly, i recently tried something similar. i went with a raster/numpy-array approach.
it works, but not nearly as fast as linda´s approach (hoping for her code/approach to be posted here as well)...i just slightly modified my stuff and timed a random 250x250 cell raster @ ~20s.

import numpy as np
import arcpy

def chkNeighbors(arr, x, y):
    try:
            
        left = arr[(x-1, y)] if x-1 >= 0 else 0
        right = arr[(x+1, y)] if x+1 <= (np.shape(arr)[0]-1) else 0
        top = arr[(x, y-1)] if y-1 >= 0 else 0
        bottom = arr[(x, y+1)] if y+1 <= (np.shape(arr)[1]-1) else 0
                
        hVal = 3
        
        if all([left > hVal, right > hVal, top > hVal, bottom > hVal]) \
            and left == right == top == bottom:       
                    
            return top
        else:
            return arr[x,y]
            
    except IndexError:
        print "error@", x,y

def doanalysis():    
    arr = arcpy.RasterToNumPyArray(r"C:\Default.gdb\Reclass_test500")
   
    print np.shape(arr)
        
    itemindex = zip(*np.where(arr))   
    i=0 
    
    for i in range(0,len(arr.flat)):
        
        arr[itemindex] = chkNeighbors(arr, *itemindex)
        
    
    ras = arcpy.NumPyArrayToRaster(arr)
    ras.save(r"C:\Default.gdb\testraster")

[JamalNUMAN]

if i understood all the posts correctly, i recently tried something similar. i went with a raster/numpy-array approach.
it works, but not nearly as fast as linda´s approach (hoping for her code/approach to be posted here as well)...i just slightly modified my stuff and timed a random 250x250 cell raster @ ~20s.

import numpy as np
import arcpy

def chkNeighbors(arr, x, y):
    try:
            
        left = arr[(x-1, y)] if x-1 >= 0 else 0
        right = arr[(x+1, y)] if x+1 <= (np.shape(arr)[0]-1) else 0
        top = arr[(x, y-1)] if y-1 >= 0 else 0
        bottom = arr[(x, y+1)] if y+1 <= (np.shape(arr)[1]-1) else 0
                
        hVal = 3
        
        if all([left > hVal, right > hVal, top > hVal, bottom > hVal]) \
            and left == right == top == bottom:       
                    
            return top
        else:
            return arr[x,y]
            
    except IndexError:
        print "error@", x,y

def doanalysis():    
    arr = arcpy.RasterToNumPyArray(r"C:\Default.gdb\Reclass_test500")
   
    print np.shape(arr)
        
    itemindex = zip(*np.where(arr))   
    i=0 
    
    for i in range(0,len(arr.flat)):
        
        arr[itemindex] = chkNeighbors(arr, *itemindex)
        
    
    ras = arcpy.NumPyArrayToRaster(arr)
    ras.save(r"C:\Default.gdb\testraster")

Hi all,

I would be very highly appreciated if you help me with writing up a code for the following statements:

1. IF the cell is water, THEN no growth is allowed at this cell.
2. IF the cell is road, THEN no growth is allowed at this cell.
3. IF cell is residential OR commercial, THEN keep this cell the same without any change.
4. IF the cell is either forest OR pasture OR row crops AND there are 4 commercial cells in the neighborhood, THEN change the cell to commercial.
5. IF the center cell is either forest OR pasture OR row crops AND there are 4 residential in the neighborhood, THEN change the cell to residential.

Many thanks for the help,

best

Jamal