在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

You may also like...