Finding adjacent neighbors for Cellular Automata

1447
15
09-26-2011 12:39 PM
LindaVasil
New Contributor III
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
Tags (2)
0 Kudos
15 Replies
StacyRendall1
Occasional Contributor III
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"...
0 Kudos
LindaVasil
New Contributor III
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
0 Kudos
AlessandroCinnirella
New Contributor III
arcpy.SelectLayerByLocation_management
not
arcpy.SelectLayerbyLocation_management

😉

ciao,
AC
0 Kudos
LindaVasil
New Contributor III
Thank you!  Good eye...you guys are awesome!

Linda
0 Kudos
LindaVasil
New Contributor III
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
0 Kudos
LindaVasil
New Contributor III
Never mind about the last post.  I just deleted the last two parameters from the command and it ran!
0 Kudos
AlessandroCinnirella
New Contributor III
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
0 Kudos
ChrisSnyder
Regular Contributor III
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...
0 Kudos
ElijahKafer
New Contributor
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
0 Kudos