在Arcgis中基于Python对地图分级别进行四色填充

四色填充是数学领域比较有名的定理，大概意思是说对于任意无飞地的多边形区域，总能选取四种颜色对每个多边形进行填充，保证相邻的多边形具有不同的颜色。在地图制图中，该定理被用于地图着色，保证只采用四种颜色而使得每个省/市/县与相邻区域具有不同的颜色。

一、项目需求

最近因项目需要，研究了地图四色填充算法。
主要问题：在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

转载自：https://blog.csdn.net/wan_yanyan528/article/details/49175673