一種基于Python語言的林業(yè)專題圖批量制圖方法實現(xiàn)

時啟龍，張 權(quán)，邱 琳，喻 俊

(江西省林業(yè)調(diào)查規(guī)劃研究院，南昌 330046)

在現(xiàn)代林業(yè)管理中，林業(yè)專題圖制作是一項十分重要的工作，幾乎涉及到林業(yè)管理的各個方面。傳統(tǒng)的林業(yè)專題圖制作首先需要設(shè)置好模板，然后再將圖紙導(dǎo)出，每一次導(dǎo)出都需要重新設(shè)置數(shù)據(jù)源或者重新設(shè)置制圖范圍，并修改制圖要素，無法實現(xiàn)批量出圖，對大批量制圖不利，不僅工作量巨大，而且需要占用大量時間和人力，時間和人力成本比較高。Python語言作為一種解釋型腳本語言，為ArcGis定制了專門的批量制圖功能，不僅可以利用已有的模板制圖，還可以重新設(shè)定制圖模板，從而大大減輕林業(yè)專題圖制作的工作量，為林業(yè)管理工作帶來很大便利[1]。本文將以公益林小班分布圖為例，詳細介紹基于Python語言的林業(yè)專題圖批量制圖方法的實現(xiàn)過程。

1 Python語言的優(yōu)勢

Python是一種不受局限、跨平臺的開源編程語言，功能強大且簡潔易讀，可伸縮性強，既適用于大型項目，和其他編程語言黏合在一起，又可單獨使用，應(yīng)用于一次性程序[2]。眾多開源的科學(xué)計算軟件包都提供了Python的調(diào)用接口，例如著名的計算機視覺庫OpenCV、三維可視化庫VTK和醫(yī)學(xué)圖像處理庫ITK等[3]。目前Python已延伸到ArcGis軟件中，Arcpy是ArcGis的一個Python包，其中包含了對地圖代數(shù)、地圖操作的支持，不僅如此，它還支持地圖的編輯處理和相關(guān)幾何操作。Python已經(jīng)實現(xiàn)了與 ArcGis軟件的高度集成化，方便用戶和開發(fā)者實現(xiàn)利用ArcGis自動化和批處理各種數(shù)據(jù)需求，大大提高了工作效率[4]。

2 數(shù)據(jù)準備

數(shù)據(jù)準備主要包括數(shù)據(jù)預(yù)處理和出圖模板制作，數(shù)據(jù)預(yù)處理主要是新增幾個中間過渡字段，以便參與后期的程序計算，而出圖模板制作主要是確定圖幅尺寸、標注樣式，以及相關(guān)的圖面整飾。

2.1 數(shù)據(jù)預(yù)處理

數(shù)據(jù)預(yù)處理，即在公益林小班矢量數(shù)據(jù)庫中增加YMIN、YMAX、XMIN、XMAX和BJ五個字段，前四個字段用于通過ArcGis自帶的字段計算器，計算公益林小班的四至坐標，即每個公益林小班的南至、北至、西至和東至坐標。第五個字段用于在全村公益林小班分布范圍超過設(shè)定比例尺需要分幅時，依次循環(huán)計算并記錄分幅頁碼，首賦值“A”以作初始標記。

2.2 制作出圖模板

出圖模板采用A3圖幅和1∶1萬標準比例尺，制圖元素主要包括圖框、圖名、圖幅號、圖幅頁碼和標注說明等其他圖面整飾元素。公益林小班標注，主要包括標注內(nèi)容、標注顏色和字體大小等，模板效果如圖1所示。

圖1 出圖模板

3 批量制圖代碼實現(xiàn)

首先打開代碼編輯窗口，導(dǎo)入Arcpy、os等模塊，設(shè)置初始環(huán)境和數(shù)據(jù)源并定義全局性變量。Arcpy是ArcGis的一個Python包，主要由數(shù)據(jù)訪問模塊、制圖模塊、網(wǎng)絡(luò)分析模塊和空間分析模塊等組成,可以完成空間地理數(shù)據(jù)處理、分析和制圖的功能。

import arcpy #導(dǎo)入站點包

from arcpy import env

import os

arcpy.CheckOutExtension("spatial")

arcpy.gp.overwriteOutput = 1

env.workspace = "c:/zygx/" #設(shè)置工作空間

gyljx = "c:/zygx/slgy.mdb/gyljx" #設(shè)置數(shù)據(jù)源

mxd = arcpy.mapping.MapDocument("CURRENT")#定義當(dāng)前文檔

公益林分布圖以村為基本單位制圖，并按鄉(xiāng)鎮(zhèn)和村分級存放，首先要訪問數(shù)據(jù)庫中的鄉(xiāng)鎮(zhèn)信息，按鄉(xiāng)鎮(zhèn)代碼執(zhí)行循環(huán)，然后訪問村信息，并在鄉(xiāng)鎮(zhèn)循環(huán)中按村代碼執(zhí)行內(nèi)循環(huán)，循環(huán)初始，分別建立對應(yīng)的文件夾以存放MXD文檔和專題圖,主要代碼如下。

xcursors=arcpy.da.SearchCursor("gyljx",["XIANG","CUN","XIANGCUN"])#定義查詢游標，遍歷公益林小班數(shù)據(jù)庫中的元素

for xcursor in xcursors: #按鄉(xiāng)鎮(zhèn)代碼建立循環(huán)

os.mkdir(r"c:/zygx/zhitu/" + xcursor[0])#建鄉(xiāng)鎮(zhèn)文件夾

xqry = "XIANG" + "=" + "'" + xcursor[0]+ "'" #設(shè)置查詢條件

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

for cursor in cursors: #按村代碼建立循環(huán)

os.mkdir(r"c:/zygx/zhitu/" + xcursor[0]+ "/" + cursor[1])#在鄉(xiāng)鎮(zhèn)文件夾下面建村文件夾

公益林小班分布圖的出圖比例設(shè)定為1∶1萬標準比例尺，在出圖時，首先要判斷在當(dāng)前圖幅內(nèi)全村小班分布范圍比例是否大于或等于1∶1萬，若大于或等于該比例尺，則意味著全村小班均在當(dāng)前圖幅范圍內(nèi)，無需分幅，可直接出圖，主要代碼如下。

tc = arcpy.mapping.ListDataFrames(mxd,"layers")[0]#定義當(dāng)前圖層組

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]#定義當(dāng)前數(shù)據(jù)框

dataframe.zoomToSelectedFeatures()#縮放至當(dāng)前所選小班范圍

dataframe.scale = dataframe.scale +(1000-dataframe.scale % 1000)#計算當(dāng)前比例尺

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)#導(dǎo)出專題圖

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

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

若在當(dāng)前圖幅內(nèi)全村公益林小班分布范圍比例小于1∶1萬比例尺，則意味著，在當(dāng)前圖幅范圍內(nèi)無法涵蓋所有小班，需要對全村公益林小班分布范圍進行分幅計算。首先要計算全村最北面小班北至點的Y值，然后以Y值為圖幅縱向范圍的最大值，并以圖幅大小為約束范圍，循環(huán)計算位于該范圍內(nèi)的公益林小班。以下代碼為第一次計算圖幅范圍內(nèi)的小班位置。

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))#計算最北側(cè)小班北至點的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 #計算當(dāng)前圖幅框內(nèi)最西側(cè)小班的西至點X值

xr =(bcursor[2]+ bcursor[3])/2 + 1900 #計算當(dāng)前圖幅框內(nèi)最東側(cè)小班的東至點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")#對查詢數(shù)據(jù)賦值

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

以第一次所計算公益林小班的中心位置為基準，循環(huán)計算圖幅范圍內(nèi)最西側(cè)的小班，每次循環(huán)后，均更新當(dāng)前圖幅范圍的中心位置，以作為下次循環(huán)中心點位置的基準，直到當(dāng)前圖幅范圍包含了最西側(cè)的小班為止，主要代碼如下。

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

以上述所計算最西側(cè)小班西至點的X值為橫軸最小值，開始循環(huán)計算圖幅范圍內(nèi)最東側(cè)小班的位置，每次循環(huán)后，均更新當(dāng)前圖幅范圍的中心位置，以作為下次循環(huán)中心點位置的基準，直到當(dāng)前圖幅范圍包含了最東側(cè)的小班為止，代碼同上。

通過上述循環(huán)計算過程，即可最大限度計算出當(dāng)前圖幅范圍內(nèi)所有的公益林小班，然后在“BJ”字段內(nèi)，賦予頁碼值，以便后續(xù)出圖時分幅出圖，并顯示分幅頁碼，主要代碼如下。

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"

按照以上所述操作，循環(huán)計算全村的公益林小班，直到完成分幅，最后按照分幅后的公益林小班范圍依次出圖，并進入下一村的數(shù)據(jù)計算，直至完成全縣各個村的公益林專題圖制作。

4 討論

傳統(tǒng)林業(yè)專題圖制作模式需要手動設(shè)置并逐張導(dǎo)出，費時費力成本高。基于ArcGis的Python語言可以為林業(yè)專題圖制作提供自動化和批量化操作，極大提高林業(yè)制圖能力和效率，且不易出錯，這在我們的工作實踐中得以驗證。此外，Python語言還可為其他空間地理數(shù)據(jù)處理提供自動化處理操作，給當(dāng)前各行業(yè)海量數(shù)據(jù)處理分析提供技術(shù)支持，是人工智能和大數(shù)據(jù)分析可資利用的編程語言。作為一種開源的解釋型腳本語言，今后應(yīng)不斷加大開發(fā)深度和應(yīng)用廣度，為提高數(shù)據(jù)處理能力和處理效率提供有力保障。