一种基于Python语言的林业专题图批量制图方法实现

时启龙，张 权，邱 琳，喻 俊

(江西省林业调查规划研究院，南昌 330046)

在现代林业管理中，林业专题图制作是一项十分重要的工作，几乎涉及到林业管理的各个方面。传统的林业专题图制作首先需要设置好模板，然后再将图纸导出，每一次导出都需要重新设置数据源或者重新设置制图范围，并修改制图要素，无法实现批量出图，对大批量制图不利，不仅工作量巨大，而且需要占用大量时间和人力，时间和人力成本比较高。Python语言作为一种解释型脚本语言，为ArcGis定制了专门的批量制图功能，不仅可以利用已有的模板制图，还可以重新设定制图模板，从而大大减轻林业专题图制作的工作量，为林业管理工作带来很大便利[1]。本文将以公益林小班分布图为例，详细介绍基于Python语言的林业专题图批量制图方法的实现过程。

1 Python语言的优势

Python是一种不受局限、跨平台的开源编程语言，功能强大且简洁易读，可伸缩性强，既适用于大型项目，和其他编程语言黏合在一起，又可单独使用，应用于一次性程序[2]。众多开源的科学计算软件包都提供了Python的调用接口，例如著名的计算机视觉库OpenCV、三维可视化库VTK和医学图像处理库ITK等[3]。目前Python已延伸到ArcGis软件中，Arcpy是ArcGis的一个Python包，其中包含了对地图代数、地图操作的支持，不仅如此，它还支持地图的编辑处理和相关几何操作。Python已经实现了与 ArcGis软件的高度集成化，方便用户和开发者实现利用ArcGis自动化和批处理各种数据需求，大大提高了工作效率[4]。

2 数据准备

数据准备主要包括数据预处理和出图模板制作，数据预处理主要是新增几个中间过渡字段，以便参与后期的程序计算，而出图模板制作主要是确定图幅尺寸、标注样式，以及相关的图面整饰。

2.1 数据预处理

数据预处理，即在公益林小班矢量数据库中增加YMIN、YMAX、XMIN、XMAX和BJ五个字段，前四个字段用于通过ArcGis自带的字段计算器，计算公益林小班的四至坐标，即每个公益林小班的南至、北至、西至和东至坐标。第五个字段用于在全村公益林小班分布范围超过设定比例尺需要分幅时，依次循环计算并记录分幅页码，首赋值“A”以作初始标记。

2.2 制作出图模板

出图模板采用A3图幅和1∶1万标准比例尺，制图元素主要包括图框、图名、图幅号、图幅页码和标注说明等其他图面整饰元素。公益林小班标注，主要包括标注内容、标注颜色和字体大小等，模板效果如图1所示。

图1 出图模板

3 批量制图代码实现

首先打开代码编辑窗口，导入Arcpy、os等模块，设置初始环境和数据源并定义全局性变量。Arcpy是ArcGis的一个Python包，主要由数据访问模块、制图模块、网络分析模块和空间分析模块等组成,可以完成空间地理数据处理、分析和制图的功能。

import arcpy #导入站点包

from arcpy import env

import os

arcpy.CheckOutExtension("spatial")

arcpy.gp.overwriteOutput = 1

env.workspace = "c:/zygx/" #设置工作空间

gyljx = "c:/zygx/slgy.mdb/gyljx" #设置数据源

mxd = arcpy.mapping.MapDocument("CURRENT")#定义当前文档

公益林分布图以村为基本单位制图，并按乡镇和村分级存放，首先要访问数据库中的乡镇信息，按乡镇代码执行循环，然后访问村信息，并在乡镇循环中按村代码执行内循环，循环初始，分别建立对应的文件夹以存放MXD文档和专题图,主要代码如下。

xcursors=arcpy.da.SearchCursor("gyljx",["XIANG","CUN","XIANGCUN"])#定义查询游标，遍历公益林小班数据库中的元素

for xcursor in xcursors: #按乡镇代码建立循环

os.mkdir(r"c:/zygx/zhitu/" + xcursor[0])#建乡镇文件夹

xqry = "XIANG" + "=" + "'" + xcursor[0]+ "'" #设置查询条件

cursors=arcpy.da.SearchCursor("gyljx",["CUN","XIANGCUN","TFH"],xqry)

for cursor in cursors: #按村代码建立循环

os.mkdir(r"c:/zygx/zhitu/" + xcursor[0]+ "/" + cursor[1])#在乡镇文件夹下面建村文件夹

公益林小班分布图的出图比例设定为1∶1万标准比例尺，在出图时，首先要判断在当前图幅内全村小班分布范围比例是否大于或等于1∶1万，若大于或等于该比例尺，则意味着全村小班均在当前图幅范围内，无需分幅，可直接出图，主要代码如下。

tc = arcpy.mapping.ListDataFrames(mxd,"layers")[0]#定义当前图层组

for lyr in arcpy.mapping.ListLayers(mxd, "", tc): #遍历图层组中的图层

if lyr.name == "gyljx":

lyr.definitionQuery = qry

arcpy.SelectLayerByAttribute_management("gyljx","NEW_SELECTION",qry)

dataframe = arcpy.mapping.ListDataFrames(mxd, "layers")[0]#定义当前数据框

dataframe.zoomToSelectedFeatures()#缩放至当前所选小班范围

dataframe.scale = dataframe.scale +(1000-dataframe.scale % 1000)#计算当前比例尺

if dataframe.scale <= 10000:

dataframe.scale = 10000

mtname = cursor[0]+ cursor[1]+ wname

tname = arcpy.mapping.ListLayoutElements(mxd)[0]

for tname in arcpy.mapping.ListLayoutElements(mxd):

if tname.name == "TMC":

tname.text = cursor[0]+ cursor[1]+ wname #修改图名

if tname.name == "TYM":

tname.text = u"共1幅，第1幅" #修改图幅页码

if tname.name == "TFHM":

tname.text = cursor[4]#修改图幅号

arcpy.mapping.ExportToJPEG(mxd, r"c:/zygx/zhitu/" + xcursor[0]+ "/" + cursor[1]+ "/" + mtname + ".jpg",resolution = 300)#导出专题图

scy = r"c:/zygx/zhitu/" + xcursor[0]+ "/" + cursor[1]+ "/" + mtname + ".mxd"

mxd.saveACopy(scy)#保存mxd文档

若在当前图幅内全村公益林小班分布范围比例小于1∶1万比例尺，则意味着，在当前图幅范围内无法涵盖所有小班，需要对全村公益林小班分布范围进行分幅计算。首先要计算全村最北面小班北至点的Y值，然后以Y值为图幅纵向范围的最大值，并以图幅大小为约束范围，循环计算位于该范围内的公益林小班。以下代码为第一次计算图幅范围内的小班位置。

if dataframe.scale > 10000:

stdness = "true"

while stdness == "true":

zxbqry = "XIANGCUN" + "=" + "'" + cursor[2]+ "'" + "AND " + "BJ" + "=" + "'" + "A" + "'"

acursors = arcpy.da.SearchCursor("gyljx",["XIANGCUN","YMAX","XMIN","XMAX","BJ"],zxbqry)

for acursor in acursors:

ymaxsz.append(acursor[1])

ymaxzb = str(max(ymaxsz))#计算最北侧小班北至点的Y值

maxqry = "XIANGCUN" + "=" + "'" + cursor[2]+ "'" + "AND " + "YMAX" + "=" + ymaxzb + "AND " + "BJ" + "=" + "'" + "A" + "'"

bcursors = arcpy.da.SearchCursor("gyljx",["XIANGCUN","YMAX","XMIN","XMAX","BJ"],maxqry)

for bcursor in bcursors:

yd = bcursor[1]-2400 #计算图幅框的纵向范围

xl =(bcursor[2]+ bcursor[3])/2-1900 #计算当前图幅框内最西侧小班的西至点X值

xr =(bcursor[2]+ bcursor[3])/2 + 1900 #计算当前图幅框内最东侧小班的东至点X值

fxbqry = "YMIN" + ">=" + str(yd)+ "AND " + "YMAX" + "<=" + ymaxzb + "AND "+ "XMIN" + ">=" + str(xl)+ "AND " + "XMAX" + "<=" + str(xr)+ "AND " + "BJ" + "=" + "'" + "A" + "'"

arcpy.SelectLayerByAttribute_management("gyljx","NEW_SELECTION",fxbqry)

arcpy.CalculateField_management("gyljx","BJ","'T'","PYTHON")#对查询数据赋值

sxbqry = "XIANGCUN" + "=" + "'" + cursor[2]+ "'" + "AND " + "BJ" + "=" + "'" + "T" + "'"

ccursors = arcpy.da.SearchCursor("gyljx",["XMIN","XMAX","BJ"],sxbqry)

for ccursor in ccursors:

xminsz.append(ccursor[0])

xmaxsz.append(ccursor[1])

xminzb = min(xminsz)

xmaxzb = max(xmaxsz)

dvalue = xmaxzb-xminzb

ldvalue =(3800-dvalue)/2

arcpy.SelectLayerByAttribute_management("gyljx","NEW_SELECTION",sxbqry)

arcpy.CalculateField_management("gyljx","BJ","'A'","PYTHON")

lstdness = "true"

zzlvalue = xl-ldvalue

以第一次所计算公益林小班的中心位置为基准，循环计算图幅范围内最西侧的小班，每次循环后，均更新当前图幅范围的中心位置，以作为下次循环中心点位置的基准，直到当前图幅范围包含了最西侧的小班为止，主要代码如下。

while lstdness == "true":

txbqry = "YMIN" + ">=" + str(yd)+ "AND " + "YMAX" + "<=" + ymaxzb + "AND " + "XMIN" + ">=" + str(zzlvalue)+ "AND " + "XMAX" + "<=" + str(xmaxzb)+ "AND " + "BJ" + "=" + "'" + "A" + "'"

arcpy.SelectLayerByAttribute_management("gyljx","NEW_SELECTION",txbqry)

arcpy.CalculateField_management("gyljx","BJ","'T'","PYTHON")

dcursors = arcpy.da.SearchCursor("gyljx",["XMIN","XMAX","BJ"],sxbqry)

for dcursor in dcursors:

xminsz.append(dcursor[0])

xmaxsz.append(dcursor[1])

xminzb = min(xminsz)

xmaxzb = max(xmaxsz)

sdvalue = xmaxzb-xminzb

if sdvalue == dvalue:

lstdness = "fasle"

else:

lstdness = "true"

zzlvalue = zzlvalue-(3800-sdvalue)/2

dvalue = sdvalue

arcpy.SelectLayerByAttribute_management("gyljx","NEW_SELECTION",sxbqry)

arcpy.CalculateField_management("gyljx","BJ","'A'","PYTHON")

rstdness = "true"

zzrvalue = xmaxzb +(3800-sdvalue)/2

以上述所计算最西侧小班西至点的X值为横轴最小值，开始循环计算图幅范围内最东侧小班的位置，每次循环后，均更新当前图幅范围的中心位置，以作为下次循环中心点位置的基准，直到当前图幅范围包含了最东侧的小班为止，代码同上。

通过上述循环计算过程，即可最大限度计算出当前图幅范围内所有的公益林小班，然后在“BJ”字段内，赋予页码值，以便后续出图时分幅出图，并显示分幅页码，主要代码如下。

arcpy.SelectLayerByAttribute_management("gyljx","NEW_SELECTION",sxbqry)

arcpy.CalculateField_management("gyljx","BJ",k,"PYTHON")#k为页码值，初始1

fcursors = arcpy.da.SearchCursor("gyljx",["BJ"],qry)

for fcursor in fcursors:

yzsz.append(fcursor[0])

if 'A' in yzsz: #条件语句，判断分幅是否完成

stdness = "true"

else:

stdness ="false"

按照以上所述操作，循环计算全村的公益林小班，直到完成分幅，最后按照分幅后的公益林小班范围依次出图，并进入下一村的数据计算，直至完成全县各个村的公益林专题图制作。

4 讨论

传统林业专题图制作模式需要手动设置并逐张导出，费时费力成本高。基于ArcGis的Python语言可以为林业专题图制作提供自动化和批量化操作，极大提高林业制图能力和效率，且不易出错，这在我们的工作实践中得以验证。此外，Python语言还可为其他空间地理数据处理提供自动化处理操作，给当前各行业海量数据处理分析提供技术支持，是人工智能和大数据分析可资利用的编程语言。作为一种开源的解释型脚本语言，今后应不断加大开发深度和应用广度，为提高数据处理能力和处理效率提供有力保障。