topic Re: Finding adjacent neighbors for Cellular Automata in Python Questions

Finding adjacent neighbors for Cellular Automata

LindaVasil — Mon, 26 Sep 2011 19:39:01 GMT

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

Re: Finding adjacent neighbors for Cellular Automata

StacyRendall1 — Sat, 11 Dec 2021 14:39:35 GMT

I think you need to make your feature class into a feature layer before you can select by attributes (making a feature layer is like adding it to the map - it makes it active). i.e.

arcpy.MakeFeatureLayer_management(shp, "shp_lyr")
#Select the sequential order of Object ID
arcpy.SelectLayerByAttribute_management("shp_lyr", "New_Selection", sql)

After that, if you need it, refer to "shp_lyr"...

Re: Finding adjacent neighbors for Cellular Automata

LindaVasil — Tue, 27 Sep 2011 12:35:31 GMT

Hi Stacy,

Thank you very much! That got me through that part of the code. I hope I can ask you something else. The select layer by location is getting an error now. Do you know what I would need to do to resolve this error:

Traceback (most recent call last):
File "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 37, in <module>
Adj_Neighbors = arcpy.SelectLayerbyLocation_management("shp_lyr", "Boundary_Touches", "shp_lyr ", " ", "New Selection")
AttributeError: 'module' object has no attribute 'SelectLayerbyLocation_management'

Does this mean that the new selection has no attributes associated with it? I'm not sure how to solve this?

Thanks for any help!

Linda

Re: Finding adjacent neighbors for Cellular Automata

AlessandroCinnirella — Tue, 27 Sep 2011 12:55:44 GMT

arcpy.SelectLayerByLocation_management
not
arcpy.SelectLayerbyLocation_management

😉

ciao,
AC

Re: Finding adjacent neighbors for Cellular Automata

LindaVasil — Tue, 27 Sep 2011 15:51:06 GMT

Thank you! Good eye...you guys are awesome!

Linda

Re: Finding adjacent neighbors for Cellular Automata

LindaVasil — Tue, 27 Sep 2011 16:08:34 GMT

So I know I'm like a child that's bugging you guys (anyone please!). So hopefully you guys will continue to pick out my silly errors. Ok, so I changed the select by location to right syntax. I changed the boundary touches to intersect because I read boundary touches only works for line and polygon shapefiles and I want to do this on a point shapefile. So the new command I have is:

#Select only the Point ID that are surrounding the cursor cell
Adj_Neighbors = arcpy.SelectLayerByLocation_management("shp_lyr", "INTERSECT", "shp_lyr", "0 meter", "New Selection")

I put 0 meter in for search distance because when I just left it as " ", the error comes back as "cannot set input parameter search_distance". Now with the 0 meter as the search distance, I get the following error:

Adj_Neighbors = arcpy.SelectLayerByLocation_management("shp_lyr", "INTERSECT", "shp_lyr", "0 meter", "New Selection")
File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\management.py", line 4381, in SelectLayerByLocation
raise e
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000800: The value is not a member of NEW_SELECTION | ADD_TO_SELECTION | REMOVE_FROM_SELECTION | SUBSET_SELECTION | SWITCH_SELECTION.
Failed to execute (SelectLayerByLocation).

I have to say "huh?" to that error. Again I appreciate anyone who would hold my hand like I'm a child and help me decipher what I'm doing wrong.

Thanks again for any help,
Linda

Re: Finding adjacent neighbors for Cellular Automata

LindaVasil — Tue, 27 Sep 2011 16:19:10 GMT

Never mind about the last post. I just deleted the last two parameters from the command and it ran!

Re: Finding adjacent neighbors for Cellular Automata

AlessandroCinnirella — Tue, 27 Sep 2011 19:28:36 GMT

FYI,

Adj_Neighbors = arcpy.SelectLayerByLocation_management("shp_lyr", "INTERSECT", "shp_lyr", "0 meter", "NEW_SELECTION")

not

Adj_Neighbors = arcpy.SelectLayerByLocation_management("shp_lyr", "INTERSECT", "shp_lyr", "0 meter", "New Selection")

ehehee,

ciao,
AC

Re: Finding adjacent neighbors for Cellular Automata

ChrisSnyder — Tue, 27 Sep 2011 20:20:22 GMT

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

Re: Finding adjacent neighbors for Cellular Automata

ElijahKafer — Mon, 28 Nov 2011 17:33:15 GMT

Hi all,

I am fresh to a lot of this subject matter, but a project I am working on in my spare time (I'm actually just a designer) is getting me into things I have no idea where to start in, but I think it is in the same vein as the land use CA model you are discussing here.

I have a basic understanding of programming (Python, Java and Processing) and I programmed a CA in Processing to test its functionality but I need to find out where to learn to program in an actual geospatial program so my analysis can be properly visualized. Would anyone point me in the right direction to start?

Many thanks,

Elijah Kafer

Re: Finding adjacent neighbors for Cellular Automata

LindaVasil — Mon, 05 Mar 2012 19:10:59 GMT

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

Re: Finding adjacent neighbors for Cellular Automata

ChrisSnyder — Mon, 05 Mar 2012 22:14:08 GMT

'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.

Re: Finding adjacent neighbors for Cellular Automata

curtvprice — Mon, 05 Mar 2012 22:51:05 GMT

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

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

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.

Re: Finding adjacent neighbors for Cellular Automata

ChrisSnyder — Mon, 05 Mar 2012 23:04:23 GMT

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

Re: Finding adjacent neighbors for Cellular Automata

RaphaelR — Sat, 11 Dec 2021 14:39:38 GMT

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

Re: Finding adjacent neighbors for Cellular Automata

JamalNUMAN — Sat, 11 Dec 2021 14:39:41 GMT

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