@jo Hi jo, thank you so much for the help, I'll try to rewrite by using the dictionary.

I'm not sure if I got the "one quarter" approach, like how could I count the neighbor numbers?

After @jo 's 1D cellular automaton, I've explored the 2d CA considering 4 neighbors and starting from one black cell in the center.

My question: if there is a better way to implement the rules without using a lot of "if statements"?

The rule is: The base 2 digits of the rule number determines the CA evolution. The last bit specifies the state of a cell if all neighbors are OFF and it too is OFF. The next to last bit specifies the state of a cell if all neighbors are OFF but the cell itself is ON. Then each earlier pair of bits specifies what should happen if progressively (Totalistic) more neighbors are black. So, bits 2^0 and 2^1 apply if none of the four neighbors are ON, bits 2^2 and 2^3 apply if one neighbor is ON, bits 2^4 and 2^5 apply if two neighbors are ON, bits 2^6 and 2^7 apply if three neighbors are ON and bits 2^8 and 2^9 apply if all four neighbors are ON. (http://oeis.org/wiki/Index_to_2D_5-Neighbor_Cellular_Automata)

Thank you!

'''
2D 5-neighbor cellular automata with drawbot, by Juan Feng

See "A New Kind of Science" by Stephen Wolfram (p.170 - 179)
http://www.wolframscience.com/nks/p171--cellular-automata/
'''

# square size
cellSize = 10

# odd numbers only 
numCell = 55 
if numCell % 2 == 0:
    numCell += 1

canvas = numCell * cellSize
size (canvas, canvas)

# type the number: 0 - 1023
rule = 462 
ruleSet = bin(rule)[2:].zfill(10)
print (ruleSet)

# starting from one active black cell in the center
grid = [[0 for x in range(numCell)]for y in range(numCell)]
grid[int((numCell-1)/2)][int((numCell-1)/2)] = 1 

# count neighbor numbers (considering Von Neumann neighbors)
def nbs(grid, r, c):
    def get(r, c):
        if 0 <= r < len(grid) and 0 <= c < len(grid[r]):
            return grid[r][c]
        else:
            return 0
    neighbors_list = [get(r-1, c), get(r, c-1), get(r, c+1), get(r+1, c)]
    return sum(map(bool, neighbors_list))

# step: 0 - inf
step = 22    
for i in range(step):
    newGrid = []
    for r,u in enumerate(grid): 
        newGrid.append([]) 
        for c,v in enumerate(u):
            if grid[r][c] == 0 and nbs(grid,r,c) == 0:
                newGrid[r].append(int(ruleSet[9]))
            if grid[r][c] == 1 and nbs(grid,r,c) == 0:
                newGrid[r].append(int(ruleSet[8]))
            if grid[r][c] == 0 and nbs(grid,r,c) == 1:
                newGrid[r].append(int(ruleSet[7]))
            if grid[r][c] == 1 and nbs(grid,r,c) == 1:
                newGrid[r].append(int(ruleSet[6]))
            if grid[r][c] == 0 and nbs(grid,r,c) == 2:
                newGrid[r].append(int(ruleSet[5]))
            if grid[r][c] == 1 and nbs(grid,r,c) == 2:
                newGrid[r].append(int(ruleSet[4]))
            if grid[r][c] == 0 and nbs(grid,r,c) == 3:
                newGrid[r].append(int(ruleSet[3]))
            if grid[r][c] == 1 and nbs(grid,r,c) == 3:
                newGrid[r].append(int(ruleSet[2]))
            if grid[r][c] == 0 and nbs(grid,r,c) == 4:
                newGrid[r].append(int(ruleSet[1]))
            if grid[r][c] == 1 and nbs(grid,r,c) == 4:
                newGrid[r].append(int(ruleSet[0])) 
    grid = newGrid 

# draw cells     
yy = 0
while yy * cellSize <= canvas - cellSize:
    xx = 0
    while xx * cellSize <= canvas - cellSize:
        if grid[xx][yy] == 1:
            fill(0)
        else:
            fill(.7)
        rect(xx*cellSize, yy*cellSize, cellSize, cellSize)
        xx += 1
    yy += 1
    
# saveImage('~/Desktop/2d_CA_test.png')