時(shí)啟龍，邱 琳，喻 俊

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

森林資源年度更新是摸清森林資源家底、進(jìn)行森林資源管理的一項(xiàng)基礎(chǔ)性工作[1]，也是掌握林地及林木資源空間分布及動(dòng)態(tài)變化情況，保持林地及林木資源數(shù)據(jù)準(zhǔn)確性和時(shí)效性的有效手段，可為本地區(qū)生態(tài)文明試驗(yàn)區(qū)建設(shè)目標(biāo)評(píng)價(jià)考核、科學(xué)發(fā)展觀綜合考核、領(lǐng)導(dǎo)干部自然資源資產(chǎn)離任審計(jì)、自然資源資產(chǎn)負(fù)債表編制等工作提供及時(shí)準(zhǔn)確的基礎(chǔ)數(shù)據(jù)，也為當(dāng)?shù)亟ㄔO(shè)項(xiàng)目使用林地行政許可和林地保護(hù)行政執(zhí)法、生態(tài)公益林與天然林資源保護(hù)管理等提供執(zhí)法依據(jù)，并可為當(dāng)?shù)卣土謽I(yè)主管部門科學(xué)決策、規(guī)范管理提供重要支撐。

森林資源年度更新每年都需要進(jìn)行大量的地理空間數(shù)據(jù)處理工作，而且這種處理過(guò)程往往具有批量化和重復(fù)性的特點(diǎn)，在以往的數(shù)據(jù)處理中，一般是運(yùn)用Arcgis軟件對(duì)單一數(shù)據(jù)集或單一過(guò)程進(jìn)行處理，不僅效率低下，而且容易出錯(cuò)[2]。Python是一種不受局限、跨平臺(tái)的開源編程語(yǔ)言，它功能強(qiáng)大且簡(jiǎn)潔易讀，可伸縮性強(qiáng)，既適用于大型項(xiàng)目，又可應(yīng)用于一次性程序[3]，目前已延伸到Arcgis軟件中，成為一種進(jìn)行數(shù)據(jù)提取、數(shù)據(jù)轉(zhuǎn)換、數(shù)據(jù)處理以及數(shù)據(jù)分析和地圖制圖自動(dòng)化的腳本語(yǔ)言[4]，在Arcgis中運(yùn)用Python語(yǔ)言可實(shí)現(xiàn)對(duì)地理數(shù)據(jù)處理的自動(dòng)化和批量化，從而大大提高了工作效率，同時(shí)計(jì)算機(jī)視覺庫(kù)OpenCV、三維可視化庫(kù)VTK和醫(yī)學(xué)圖像處理庫(kù)ITK等眾多開源的科學(xué)計(jì)算軟件包都提供了Python調(diào)用接口[5]。因此，我們?cè)谏仲Y源年度更新多個(gè)環(huán)節(jié)中進(jìn)行了基于Acrgis的Python腳本批處理操作，也取得了滿意效果。

1 遙感影像裁剪

森林資源年度更新一個(gè)重要的基礎(chǔ)性工作就是以縣域?yàn)閱卧獙?duì)遙感影像進(jìn)行裁剪，在以往的工作過(guò)程中，影像裁剪工作均是手動(dòng)操作，首先需要參照影像圖幅框中各影像的空間位置，在全省的影像中找到與某縣域相交的所有影像，然后在Arcgis中進(jìn)行裁剪，一個(gè)單元裁剪完成后，再手動(dòng)進(jìn)行下一個(gè)單元的裁剪配置，不僅工作量大，效率低，而且在查找影像過(guò)程中容易混淆出錯(cuò)。而采用Python腳本語(yǔ)言，則只需要設(shè)置好全省縣界面、影像圖幅框和全省遙感影像圖的存放位置，就可以實(shí)現(xiàn)對(duì)全省一百多個(gè)單位和幾百景遙感影像進(jìn)行自動(dòng)裁剪。具體實(shí)現(xiàn)代碼如下：

#導(dǎo)入站點(diǎn)包和環(huán)境參數(shù)

import arcpy

from arcpy import env

import os

#定義初始條件

arcpy.CheckOutExtension("spatial")

arcpy.gp.overwriteOutput = 1

env.workspace = "c:/zygx"

xjm = "c:/zygx/xjm.shp"

tfk = "c:/zygx/tufukuang.shp"

xnum = "XIAN"

#加載全省縣界面

arcpy.MakeFeatureLayer_management(xjm, "xjm")

#定義查詢游標(biāo)，遍歷縣界面中的元素

xcursors = arcpy.da.SearchCursor("xjm", ["XIAN","XIAN_NAME"])

for xcursor in xcursors:

qry= xnum + "=" + "'" + xcursor[0] + "'"

#選中裁剪單元

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

#定義裁剪后影像的存放路徑

folder = "c:/zygx/cj/ " + xcursor[1]

os.makedirs(folder)

#定義查詢游標(biāo)，遍歷圖幅框中的元素

qryy = "XIAN" + "= '" + xcursor[0] + "'"

tcursors = arcpy.da.SearchCursor(tfk, ["XIAN","TNAME"],qryy)

for tcursor in tcursors :

inpt = "c:/zygx/yx/" + tcursor[1] + ".ers"

outpt = folder + "/" + tcursor[1] + "_clip.img"

#裁剪影像

arcpy.Clip_management(inpt,"",outpt,"xjm","","NONE")

print(tcursor[1] + "Has Done")

print(xcursor[1] + "All Done")

print("All Done")

2 遙感影像拼接

遙感影像拼接作為影像處理過(guò)程中的一個(gè)重要環(huán)節(jié)，比較費(fèi)事費(fèi)力，而借助Arcgis的Python腳本語(yǔ)言可以實(shí)現(xiàn)對(duì)全省一百多個(gè)單位的遙感影像圖進(jìn)行自動(dòng)化拼接，拼接過(guò)程簡(jiǎn)單易行，只需要設(shè)置好輸入和輸出數(shù)據(jù)的存放路徑即可。具體實(shí)現(xiàn)代碼如下：

import arcpy

from arcpy import env

import os

arcpy.CheckOutExtension("spatial")

arcpy.gp.overwriteOutput = 1

env.workspace = "c:/zygx/"

#遍歷影像存放文件夾列表

filelists = os.listdir("c:/zygx/cj")

for filelist in filelists:

#獲取文件夾中影像列表

env.workspace = "c:/zygx/cj/" + filelist

rasters = arcpy.ListRasters("*","img")

#影像拼接

yxNmae = filelist + ".img"

arcpy.MosaicToNewRaster_management(rasters, env.workspace, yxName, "", "8_BIT_UNSIGNED", "2", "3", "LAST", "FIRST")

3 森林資源小班數(shù)據(jù)庫(kù)分發(fā)

森林資源小班數(shù)據(jù)庫(kù)分發(fā)是森林資源年度更新的一項(xiàng)基礎(chǔ)性準(zhǔn)備工作，整個(gè)森林資源年度更新過(guò)程均在分發(fā)的森林資源小班數(shù)據(jù)庫(kù)上進(jìn)行。在以往的小班數(shù)據(jù)分發(fā)過(guò)程中，一般是先建立以單位代碼和名稱為文件名的文件地理數(shù)據(jù)庫(kù)或個(gè)人地理數(shù)據(jù)庫(kù)，然后在全省數(shù)據(jù)庫(kù)中按照單位代碼或單位名稱選中本單位數(shù)據(jù)并導(dǎo)入所建地理數(shù)據(jù)庫(kù)中，過(guò)程雖簡(jiǎn)單，但比較費(fèi)時(shí)，且效率低。而借助Python腳本語(yǔ)言可以實(shí)現(xiàn)對(duì)這一操作的自動(dòng)化運(yùn)行。具體代碼如下：

import arcpy

from arcpy import env

arcpy.CheckOutExtension("spatial")

arcpy.gp.overwriteOutput = 1

env.workspace = "c:/zygx/"

zym = "c:/zygx/jxzym2018.gdb/jxzym2018"

xjm = "c:/zygx/xjm.shp"

zyname = u'資源面2018'

xnum = "XIAN"

#加載森林資源小班數(shù)據(jù)庫(kù)

arcpy.MakeFeatureLayer_management(zym,"zym")

#將加載的資源小班數(shù)據(jù)庫(kù)隱藏不顯示

mxd = arcpy.mapping.MapDocument("CURRENT")

ryx = arcpy.mapping.ListDataFrames(mxd,"layers")[0]

for lyr in arcpy.mapping.ListLayers(mxd, "", ryx):

if lyr.name == "zym":

lyr.visible = "false"

#定義查詢游標(biāo)，遍歷數(shù)據(jù)庫(kù)中元素

cursors = arcpy.da.SearchCursor(xjm,["XIAN","XIAN_NAME"])

for cursor in cursors:

#建立文件地理數(shù)據(jù)庫(kù)

arcpy.CreateFileGDB_management("c:/zygx/zysj/",cursor[0] + cursor[1])

#選中森林資源小班數(shù)據(jù)庫(kù)

qry = xnum + "=" + "'" + cursor[0] + "'"

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

#導(dǎo)出所選數(shù)據(jù)

outpt = "c:/zygx/zysj/" + cursor[0] + cursor[1] + ".gdb"

tname = cursor[1] + zyname

arcpy.FeatureClassToFeatureClass_conversion("zym",outpt,tname)

print(cursor[1] + "Has Done")

print("All Done")

4 外業(yè)調(diào)查圖批量制圖

調(diào)查圖制作是森林資源調(diào)查中的一項(xiàng)重要工作，幾乎應(yīng)用于林業(yè)管理的各個(gè)方面，在傳統(tǒng)的森林資源調(diào)查過(guò)程中，外業(yè)調(diào)查圖紙的制作均是在設(shè)置好模板的基礎(chǔ)上，一張一張的導(dǎo)出，而無(wú)法實(shí)現(xiàn)批量出圖，這對(duì)于大批量的制圖工作是不利的，不僅工作量巨大，而且需要占用大量的時(shí)間，時(shí)間成本較高。Python腳本語(yǔ)言為Arcgis定制了專門的批量地圖制圖功能，不僅可以利用已有的模板制圖，而且可以重新設(shè)定制圖模板。由于森林資源年度更新外業(yè)調(diào)查圖圖式一致，本項(xiàng)目則借助已有的制圖模板，利用Python腳本的DataDrivenPages類并結(jié)合Arcgis中的Data Driven Pages工具，來(lái)實(shí)現(xiàn)對(duì)外業(yè)調(diào)查圖紙的批量化和自動(dòng)化制作。具體實(shí)現(xiàn)代碼如下：

import arcpy

from arcpy import env

arcpy.CheckOutExtension("spatial")

arcpy.gp.overwriteOutput = 1

env.workspace = "c:/zygx/"

#定義當(dāng)前文檔

mxd = arcpy.mapping.MapDocument("CURRENT")

#依據(jù)標(biāo)準(zhǔn)圖幅框切換制圖范圍

for pageNum in range(1,mxd.dataDrivenPages.pageCount):

mxd.dataDrivenPages.currentPageID=pageNum

mapName = mxd.dataDrivenPages.pageRow.getValue(mxd.dataDrivenPages.pageNameField.name)

raster = "c:/zygx/ditu" + mapName + ".tif"

#依據(jù)圖幅框加載底圖

arcpy.MakeRasterLayer_management(raster,mapName)

#修改圖名

elm = arcpy.mapping.ListLayoutptElements(mxd)[0]

for elm in arcpy.mapping.ListLayoutptElements(mxd):

if elm.name == "ZTH":

elm.text = mapName

#導(dǎo)出地圖

arcpy.mapping.ExportToJPEG(mxd, r"c:/zygx/zhitu/" + mapName + ".jpg",resolution = 300)

#保存MXD文件

scy = r"c:/zygx/zhitu/" + mapName + ".mxd"

mxd.saveACopy(scy)

#刪除所加載底圖

ryx = arcpy.mapping.ListDataFrames(mxd,"layers")[0]

for lyr in arcpy.mapping.ListLayers(mxd, "", ryx):

if lyr.name == mapName:

arcpy.mapping.RemoveLayer(ryx, lyr)

print(mapName + "Has Done")

print("All Done")

5 結(jié)語(yǔ)

Python語(yǔ)言是當(dāng)今最受程序設(shè)計(jì)語(yǔ)言之一，它是一種動(dòng)態(tài)的、面向?qū)ο蟮哪_本語(yǔ)言，是不受局限、跨平臺(tái)的開源泉編程語(yǔ)言，功能強(qiáng)大，簡(jiǎn)潔易懂，應(yīng)用廣泛，不僅與Web、Internet聯(lián)合開發(fā)和軟件開發(fā)等，而且在科學(xué)計(jì)算、統(tǒng)計(jì)分析和教育方面發(fā)揮了獨(dú)特的作用。本研究在Arcgis軟件系統(tǒng)中應(yīng)用Python語(yǔ)言解決了森林資源數(shù)據(jù)的年度更新，取得了滿意的效果。