四色填充是數學領域比較有名的定理，大概意思是說對於任意無飛地的多邊形區域，總能選取四種顏色對每個多邊形進行填充，保證相鄰的多邊形具有不同的顏色。在地圖製圖中，該定理被用於地圖著色，保證只採用四種顏色而使得每個省/市/縣與相鄰區域具有不同的顏色。

一、專案需求

最近因專案需要，研究了地圖四色填充演算法。
主要問題：在web製圖中，地圖被切成向量瓦片用於渲染和顯示，這就失去了地圖的拓撲關係，如何來計算不同省/市/縣的顏色以對其進行填充？
解決思路：地圖只有在切片前具有完整的拓撲關係，因此可以假設四種顏色代號為1、2、3、4，在地圖切片前提前計算好顏色代號填入屬性表，這樣就能在web地圖中直接根據屬性表賦顏色值。
具體方案：

地圖為shp格式，分省、市、縣三個級別，因此在arcgis中進行顏色計算，Arcgis10以後版本主要基於Python製作指令碼工具，因此利用Python很容易直接對shp資料進行操作。

二、相關的方法

1 四色填充演算法——回溯法

由於不需要過於複雜的演算法，僅能實現相關功能就行，因此採用回溯演算法，具體演算法參考一篇部落格和王清平的一篇論文。

2 ArcGIS生成鄰接表

回溯演算法首先需要一個鄰接矩陣，用來表示哪些節點是相鄰關係，這種鄰接關係可以用ArcGIS生成，可參考這篇文章：如何把相鄰圖斑的屬性新增到某個欄位中。
(1) 將目標shp檔案匯出一份副本備用（以省級為例，縣級和市級類似）

(2) 在Toolbox中選擇分析工具——疊置分析——空間連線（Spatial Join），輸入要連線的兩個圖層和輸出資訊如下圖所示：

(3) 在“連線要素的欄位對映中刪除多餘的欄位，並點選右邊的”+”新增一個connection欄位，用於儲存鄰接表，欄位型別為文字，長度根據需要選擇，合併規則為“連線(join)”，分隔符為空格。

(4) 右鍵單擊剛才新建的欄位，選擇“新增輸入欄位”，這裡新增用於計算鄰接資訊的標識欄位，即每個省的ID，本例中為“PAC”欄位，該欄位為省編碼。


(5) 點選確定進行空間連線操作，完成後屬性表中“connection”欄位中即為與之相鄰的省份編號。

3 基於Python編寫工具計算每個省份的顏色

(1) 引數輸入
在ArcGIS中使用acrpy模組訪問ArcGIS的相關操作，因此首先要引入acrpy模組。該工具需要四個引數：需要進行操作的要素、要素中每個省份的標識欄位（ID）、鄰接表字段、計算單位（可選引數，該欄位主要考慮到市級和縣級多邊形太多會導致演算法不收斂，需要逐省或逐市計算）。

import arcpy
import math
import string
from numpy import *
InputFeature = arcpy.GetParameterAsText(0)
UniqueField = arcpy.GetParameterAsText(1)
ConnectField = arcpy.GetParameterAsText(2)
level = arcpy.GetParameterAsText(3)

(2) 新建“color”欄位，用於標識顏色

arcpy.AddField_management(InputFeature,"color","SHORT")

(3) 讀取欄位的屬性值

rows = arcpy.UpdateCursor(InputFeature)
for row in rows:
    N=N+1
    U.append(str(row.getValue(UniqueField)))
    C.append(str(row.getValue(ConnectField)))
    S.append(row.getValue(level))

(4) 計算鄰接矩陣
此處注意事項：ArcGIS中生成的鄰接表包含節點本身，例如1號多邊形的鄰接多邊形中會有1號自己，這與實際是不符的，解決方法可以在屬性表中提前去掉其本身的ID，或者在指令碼中忽略這個值，這裡採用後者。

n=len(u)
    mat=zeros([n,n],int)
    for i in range(0,n):
        #arcpy.AddMessage(c[i])
        tem = c[i].split(" ")
        for j in tem:
            if (j in u):
                ind=u.index(j)
                if(ind != i):
                    mat[i][ind]=1

(5) 演算法部分，通過回溯演算法計算顏色值，用colorIndex[]陣列儲存

    #計算顏色
    maxColor=4
    colorIndex = ones(n,int)
    I=1
    colorI=1
    #arcpy.AddMessage(maxColor)
    while (I<n and I>=0):
        arcpy.AddMessage(str(I))
        while (colorI<=maxColor and I<n):
            for k in range(0,I):
                if(mat[k][I] and colorIndex[k]==colorI):
                    k=k-1
                    break
            if((k+1)==I):
                colorIndex[I]=colorI
                colorI=1
                I=I+1
            else:
                colorI=colorI+1
        if(colorI>maxColor):
            #arcpy.AddMessage(str(I))
            I=I-1
            colorI=colorIndex[I]+1

(6) 給新建的“color”欄位賦值，將計算的顏色值填入color欄位

rows = arcpy.UpdateCursor(InputFeature)
    for row in rows:
        if (S[j]==each_sheng):
            row.color = colorIndex[i]
            rows.updateRow(row)
            i=i+1
        j=j+1

4 在ArcGIS中新增指令碼工具

(1) 在ArcMap的目錄中右鍵單擊Toolbox.tbx，選擇新增指令碼

(2) 根據嚮導新增指令碼工具，輸入的引數有四個，順序必須跟指令碼中的引數順序保持一致。也可以在新增指令碼後在屬性中設定引數。

5 執行指令碼工具

輸入引數如下：

點選確定後，指令碼開始執行，如果未出現錯誤會出現執行成功的提示，這時在屬性表中的color欄位就會看到計算好的顏色值。

利用color欄位對顏色進行符號化，如下圖所示，縣級和市級執行方法一樣，只是需要設定第四個引數，市級第四個引數應設定為”省”，縣級第四個引數應設定為”市”。
省級圖：

市級圖：

縣級圖：

附錄

Python指令碼的完整程式碼：

import arcpy
import math
import string
from numpy import *
InputFeature = arcpy.GetParameterAsText(0)
UniqueField = arcpy.GetParameterAsText(1)
ConnectField = arcpy.GetParameterAsText(2)
level = arcpy.GetParameterAsText(3)
#先檢查color欄位是否存在，不存在則建立該欄位
arcpy.AddField_management(InputFeature,"color","SHORT")

U = []
C = []
S = []
N = 0    #圖層中多邊形的個數
rows = arcpy.UpdateCursor(InputFeature)
#讀取資料
if (level):
    for row in rows:
        N=N+1
        U.append(str(row.getValue(UniqueField)))
        C.append(str(row.getValue(ConnectField)))
        S.append(row.getValue(level))
else:
    for row in rows:
        N=N+1
        U.append(str(row.getValue(UniqueField)))
        C.append(str(row.getValue(ConnectField)))
        S.append(u'全國')

sheng =list(set(S))
for each_sheng in sheng:
    u=[]
    c=[]
    arcpy.AddMessage(u'正在計算：'+each_sheng)
    for i in range(0,N):
        if (S[i]==each_sheng):
            u.append(U[i])
            c.append(C[i])
    #生成鄰接矩陣
    n=len(u)
    mat=zeros([n,n],int)
    for i in range(0,n):
        #arcpy.AddMessage(c[i])
        tem = c[i].split(" ")
        for j in tem:
            if (j in u):
                ind=u.index(j)
                if(ind != i):
                    mat[i][ind]=1
                    #arcpy.AddMessage("mat["+str(i)+"]["+str(ind)+"]")
    #計算顏色
    maxColor=4
    colorIndex = ones(n,int)
    I=1
    colorI=1
    #arcpy.AddMessage(maxColor)
    while (I<n and I>=0):
        arcpy.AddMessage(str(I))
        while (colorI<=maxColor and I<n):
            for k in range(0,I):
                if(mat[k][I] and colorIndex[k]==colorI):
                    k=k-1
                    break
            if((k+1)==I):
                colorIndex[I]=colorI
                colorI=1
                I=I+1
            else:
                colorI=colorI+1
        if(colorI>maxColor):
            #arcpy.AddMessage(str(I))
            I=I-1
            colorI=colorIndex[I]+1


    i=0
    j=0
    rows = arcpy.UpdateCursor(InputFeature)
    for row in rows:
        if (S[j]==each_sheng):
            row.color = colorIndex[i]
            rows.updateRow(row)
            i=i+1
        j=j+1