1. 程式人生 > >python在arcgis中案例開發(空間求交、連線以及excel資料匯出)

python在arcgis中案例開發(空間求交、連線以及excel資料匯出)

這是使用python寫的第二工具了,可以說這麼計算機語言也是初次接觸,還好使用過c、 c#、JAVA等計算機語言,所以在使用python的使用也不是完全找不到北。

這是國慶之前做的一個專案。首先我來說一下我所用到的資料:山最高點(shapefile點資料)、山底座資料(shapefile面數據,屬於一座山底部範圍面)、一個區域內的地形資料(有喀斯特沖積平原、丘陵等等多種型別資料),該資料為shapefile面數據,地理國情資料(各種地類,水田,旱地之類的shapefile面數據)。

我最終需要提交的成果是按照山底座、地形資料、地理國情資料利用arcgis提供的求交函式,將結果按照地形分類統計該地形類具有的各種地理國情要素分別佔有總面積分別為多少以excel資料表格匯出。按照arcgis在arcpy提供的api介面可以將三種資料進行求交,但是需要注意求交的順序,如果某些資料放在求交陣列的最開始有可能會在求交的過程中丟失一些欄位。我們在這裡根據需求,需要山地形資料名稱、地理國情要素編碼欄位用於地形分類、和所屬哪種地類(根據地類程式碼進行歸類)。按照山底座資料,山地形資料、國情資料一起求交之後可獲取相交資料,然後再利用arcpy的遊標,遍歷求交後的屬性shapefile資料,利用@shapearea標識求得對應屬性面的面積,再利用python的字典型別,分別儲存各種山地形資料。按照判斷字典中是否將已經歸類的地形資料型別載入到字典中的模式,即如果山地形資料地類沒有,則以山地形名稱作為鍵和該山地形內的某一個國情地類面積作為值新增對應的鍵值資料,如果已經載入有了之前的資料,那麼只要重新提取之前某種山地類的資料,然後將新值與該值累加,刪掉字典中對應的鍵值資料,又重新載入對應鍵值該山地類資料即可。整個過程編寫程式碼,經過測試,開發的外掛在求交方面執行比較花費時間,大約需要6到7分鐘,而統計山地形資料部分並不需要花太多的時間。

該基於arcgis外掛工具在執行中,由於本人是利用山地形是否有山頂點,然後將鍵值存放到python的字典中。所以再後來的統計過程中,將各種地理國情地類要素一一相加後,並沒有和山底座面積相等。經過排查,是兩方面的問題導致的。第一個原因是,地類程式碼沒有完全列舉,導致在統計的時候,程式將不在列舉範圍內的國情地類資料累計起來。第二個原因是,統計的山頂點資料中有些地形資料是沒有相應的山頂點的,比如在水域資料中的水庫地形就沒有山頂點,導致在國情資料的統計中少了某一無山頂點型別的山地形資料。在外掛工具開發過程中結合了資料,進行了國情要素程式碼的補充,以及地形資料的鍵值追加,最後才使得國情要素資料相加後,總面積與山底座資料面積相等。

在excel表格匯出方面,是藉助了python的xlwt類庫,對excel進行表格的建立以及填寫對應值。這裡有遇到一些問題就是,一定要注意編碼,我們這需要將編碼轉為UTF-8資料。為什麼這麼說,這會在字串的判斷中,如果按照其他計算機語言,C#、JAVA之類思想來判斷是否字元相同,是極易出錯的。再一個就是,使用阿拉伯數字插入到xlwt建立的表格似乎不報錯,但是如果是ascil碼,那就會報錯了。好了說了這麼多,我們來看一下,最後開發的工具介面,我相信這才是每一個從事地理資訊的GISer最為迫切與激動的。

以上的工具是利用之前講解的pycharm在arcpy開發中arcgis工具箱打包進行了打包,具體程式碼在arcgis中打包流程,大家也可以查閱相關的資料。下面來說明一下提供的原始碼。

首先是利用了空間連線函式SpatialJoin_analysis,判斷某山地形類中所具有的山頭(山頂點)數。然後在利用Intersect_analysis函式進行求交判斷。最終在選擇的【最終處理結果】資料夾中分別建立了shp、txt、xls資料夾,儲存了處理的中間檔案shapefile,統計的txt檔案、以及表格資料xls檔案。可參看如下圖所示。

其中txt的統計,本質就是講xls的內容寫到txt檔案而已,無他。下面來分別看一下txt和xls都是些什麼資料結果。

下面來看最為核心的部分,就是程式碼了。首先大家可以看一下,python中字典、類、函式,以及檔案的建立與刪除相關知識。以及在軟體設計的時候,最好是將變化的部分給抽離出來,我這裡將我們需要的地理國情要素地理程式碼作為變化部分給抽離出來,這樣是方便對地理程式碼的配置,以及錯誤的檢查,當然這裡談不上模式,只是程式設計中的一些經驗積累而已。

然後下面的這一部分程式碼就是關於,arcgis中python程式碼如何打包的相應配置,這對於初學arcpy的夥伴而言,可以借鑑一下,當然更多引數方面的配置,大家還是可以結合之前的部落格。好了,不管怎麼說,大家可以細細去體會程式碼其中所要表達的意思。至此整個使用python方面在本次arcgis開發就說明完了,可能有些小夥伴還不知所云。可以私底下交流。

功能核心程式碼:

# coding:utf-8
import arcpy
import xlrd
import xlwt
import sys
import os
reload(sys)
sys.setdefaultencoding('utf8')

def debug(messages, txt):
    if messages is None:
        print txt
    else:
        messages.addMessage(txt)

class LandInfo:
    # 山頭個數
    STCNT=0
    # 種植土地
    ZZTD=0
    # 林草覆蓋
    LCFG=0
    # 房屋建築區
    FWJZQ=0
    # 鐵路與道路
    TLYDL=0
    # 構築物
    GZW=0
    # 人工挖掘地
    RGWJD=0
    # 荒漠與裸露地
    FMYLLD=0
    #水域
    SY=0;

ZZTDListCode=['0100','0110','0120','0130','0131','0132','0133','0140','0150','0160','0170','0180','0190','0191','0192','0193','0210','0211','0212','0213','0220','0230','0240','0250','0260','0290','0291','0292','0293']#種植土地程式碼
LCFGListCode=['0300','0310','0311','0312','0313','0320','0321','0322','0323','0330','0340','0350','0360','0370','0380','0390','0391','0392','0393','03A0','03A1','03A2','03A3','03A4','03A9','0410','0411','0412','0413','0420','0421','0422','0424','0429'] # 林草覆蓋程式碼
FWJZQListCode=['0500','0510','0511','0512','0520','0521','0522','0530','0540','0541','0542','0543','0544','0550']#房屋建築區程式碼
TLYDLListCode=['0600','0610','0601']#鐵路與道路程式碼
GZWListCode=['0700','0710','0711','0712','0713','0714','0715','0716','0717','0718','0719','0720','0721','0740','0750','0760','0761','0762','0763','0769','0770','0780','0790']#構築物程式碼
RGWJDListCode=['0800','0810','0811','0812','0813','0814','0815','0819','0820','0821','0822','0829','0830','0831','0832','0833','0839','0890']#人工挖掘地程式碼
FMYLLDListCode=['0900','0910','0920','0930','0940','0950']#荒漠與裸露地程式碼
SYListCode=['1000','1001','1012','1050']#水域程式碼

#ZZTDListCode=['0210','0211','0212','0213','0220','0230','0240','0250','0260','0290','0291','0292','0293','0100','0110','0120']#種植土地
#LCFGListCode=['0310','0311','0312','0313','0320','0321','0330','0340','0360','0370','0380','0400','0410','0411','0412','0413','0420','0421','0422','0424','0429']# 林草覆蓋程式碼
#FWJZQListCode=['0510','0511','0512','0520','0521','0522','0530','0540','0541','0542','0543','0544','0550','']

#計算面積
def caculateArea(exportData,flagSXLB,flagCode,flagArea,ZZTDListCode,LCFGListCode,FWJZQListCode,TLYDLListCode,
                 GZWListCode,RGWJDListCode,FMYLLDListCode,SYListCode):
    #有些地形並沒有山頭,重追加key
    if flagSXLB not in exportData.keys():
        myLandInfo = LandInfo()
        exportData[flagSXLB] = myLandInfo
        whichClassCode(myLandInfo, flagCode, flagArea, ZZTDListCode, LCFGListCode, FWJZQListCode, TLYDLListCode,
                       GZWListCode, RGWJDListCode, FMYLLDListCode, SYListCode)
        exportData[flagSXLB] = myLandInfo
    else:
      for key, value in exportData.items():
         if key==flagSXLB:
            myLandInfo=value
            del exportData[key]
            whichClassCode(myLandInfo,flagCode,flagArea,ZZTDListCode,LCFGListCode,FWJZQListCode,TLYDLListCode,
                 GZWListCode,RGWJDListCode,FMYLLDListCode,SYListCode)
            exportData[key]=myLandInfo
            break

#判斷屬於哪個地類
def whichClassCode(landInfo,flagCode,flagArea,ZZTDListCode,LCFGListCode,FWJZQListCode,TLYDLListCode,
                 GZWListCode,RGWJDListCode,FMYLLDListCode,SYListCode):
    if isSpecificLand(ZZTDListCode,flagCode):
        landInfo.ZZTD=landInfo.ZZTD+flagArea


    elif isSpecificLand(LCFGListCode,flagCode):
        landInfo.LCFG=landInfo.LCFG+flagArea

    elif isSpecificLand(FWJZQListCode,flagCode):
        landInfo.FWJZQ=landInfo.FWJZQ+flagArea

    elif isSpecificLand(TLYDLListCode,flagCode):
        landInfo.TLYDL=landInfo.TLYDL+flagArea

    elif isSpecificLand(GZWListCode,flagCode):
        landInfo.GZW=landInfo.GZW+flagArea

    elif isSpecificLand(RGWJDListCode,flagCode):
        landInfo.RGWJD=landInfo.RGWJD+flagArea

    elif isSpecificLand(FMYLLDListCode,flagCode):
        landInfo.FMYLLD=landInfo.FMYLLD+flagArea

    elif isSpecificLand(SYListCode,flagCode):
        landInfo.SY=landInfo.SY+flagArea
    else:
        0

#刪除某一目錄下的所有檔案
def delDirAllFiles(dir):
    if os.path.exists(dir):
        files = list_all_files(dir)
        for tempfile in files:
            os.remove(tempfile)

def list_all_files(rootdir):
    _files = []
    list = os.listdir(rootdir)  # 列出資料夾下所有的目錄與檔案
    for i in range(0, len(list)):
        path = os.path.join(rootdir, list[i])
        if os.path.isdir(path):
            _files.extend(list_all_files(path))
        if os.path.isfile(path):
            _files.append(path)
    return _files

def text_create(fileName,dir):
    delDirAllFiles(dir)
    full_path=dir+fileName+".txt"
    file = open(full_path, 'w')
    file.close()
    return full_path

#判斷是否為特定地類
def isSpecificLand(listDLCode,flagCode):
       isSpecific=False
       for tmpCode in listDLCode:
         if tmpCode==flagCode:
             isSpecific=True
             break
       return isSpecific


def insertExcelData(strSXLB,row,worksheet,landInfo):
    worksheet.write(row, 0, strSXLB.encode("utf-8"))  # 山形類別
    areaTotal=landInfo.ZZTD+landInfo.LCFG+landInfo.FWJZQ+landInfo.TLYDL+landInfo.GZW+landInfo.RGWJD+landInfo.FMYLLD+landInfo.SY
    worksheet.write(row,1,landInfo.STCNT)#山頭數
    worksheet.write(row,2, areaTotal)    #山頭面積
    worksheet.write(row,3, areaTotal)#地類總面積
    worksheet.write(row,4, landInfo.ZZTD)#種植土地
    worksheet.write(row,5, landInfo.LCFG)#林草覆蓋
    worksheet.write(row,6, landInfo.FWJZQ)#房屋建築區
    worksheet.write(row,7, landInfo.TLYDL)#鐵路與道路
    worksheet.write(row,8, landInfo.GZW)#構築物
    worksheet.write(row,9, landInfo.RGWJD)#人工堆掘地
    worksheet.write(row,10, landInfo.FMYLLD)#荒漠與裸露地表
    worksheet.write(row,11, landInfo.SY)#水域



def caculteFinalArea(landInfo,totalLandinfo):
    totalLandinfo.STCNT=totalLandinfo.STCNT+landInfo.STCNT
    totalLandinfo.ZZTD=totalLandinfo.ZZTD+landInfo.ZZTD
    totalLandinfo.LCFG=totalLandinfo.LCFG+landInfo.LCFG
    totalLandinfo.FWJZQ=totalLandinfo.FWJZQ+landInfo.FWJZQ
    totalLandinfo.TLYDL=totalLandinfo.TLYDL+landInfo.TLYDL
    totalLandinfo.GZW=totalLandinfo.GZW+landInfo.GZW
    totalLandinfo.RGWJD=totalLandinfo.RGWJD+landInfo.RGWJD
    totalLandinfo.FMYLLD=totalLandinfo.FMYLLD+landInfo.FMYLLD
    totalLandinfo.SY=totalLandinfo.SY+landInfo.SY

def mkdir(path):
    #path = path.strip()
    #path = path.rstrip("\\")

    if not os.path.isdir(path):
        print "xxx"
    isExists = os.path.exists(path)
    if not isExists:
        os.makedirs(path)
        print path + ' 建立成功'
        return True
    else:
        print path + ' 目錄已存在'
        return False

# stmData="C:/Users/qrb_PC/Desktop/STProject/statistics/福泉山頭資料/reslut.gdb/fwm"#山頭範圍面
# stdData="C:/Users/qrb_PC/Desktop/STProject/statistics/福泉山頭資料/reslut.gdb/std"#山頭點資料
# sxData="C:/Users/qrb_PC/Desktop/STProject/statistics/合併BGLA.gdb/BGLA"#山形資料
# gqData="C:/Users/qrb_PC/Desktop/STProject/statistics/福泉市國情資料/522702FQN.gdb/LCA"#國情資料
#
# rootPath="C:/Users/qrb_PC/Desktop/STProject/statistics/result2/"





def do(stdData,stmData,sxData,gqData,rootPath):
    messages = None
    debug(messages,"get root folder is "+rootPath)
    shpRoot = rootPath + "/shp/"
    xlsRootPath = rootPath + "/xls/"
    txtRootPath = rootPath + "/txt/"

    debug(messages,"get shpRoot folder is "+shpRoot)
    debug(messages,"get xlsRootPath folder is "+xlsRootPath)
    debug(messages,"get txtRootPath folder is "+txtRootPath)

    mkdir(shpRoot)
    mkdir(xlsRootPath)
    mkdir(txtRootPath)

    insertGQShp = shpRoot + "insertGQshp.shp"
    insertSTShp = shpRoot + "insertSTshp.shp"

    delDirAllFiles(shpRoot)
    insertGQShp = shpRoot + "insertGQshp.shp"
    insertSTShp = shpRoot + "insertSTshp.shp"

    insertSTSXShp = shpRoot + "insertSTSXshp.shp"
    #刪掉txt目錄下所有的檔案
    delDirAllFiles(txtRootPath)
    txtPath = text_create("result", txtRootPath)
    #輸出資料
    exportData = {}

    #連線求交
    arcpy.SpatialJoin_analysis(sxData,stdData,insertSTShp,
                               join_operation="JOIN_ONE_TO_ONE",
                               join_type="KEEP_COMMON",
                               match_option="INTERSECT",
                               search_radius="#",
                               distance_field_name="#")

    cursor = arcpy.da.SearchCursor(insertSTShp, ['NAME','Join_Count'])
    for row in cursor:
        flagSXLB = row[0]
        flagSXLB=flagSXLB.encode("utf-8")

        flagCNT = row[1]
        if flagSXLB not in exportData.keys():
            landInfo=LandInfo()
            landInfo.STCNT=flagCNT
            exportData[flagSXLB]=landInfo
        else:
            tmpLandInfo=exportData[flagSXLB]
            del exportData[flagSXLB]
            tmpLandInfo.STCNT=tmpLandInfo.STCNT+flagCNT
            exportData[flagSXLB]=tmpLandInfo

    f=file(txtPath,"a+")# 追加的方式讀寫
    f.write("------------------------山頭點個數統計開始------------\n")
    for key,value in exportData.items():
        #print "山形型別:"+key
        f.write("山形型別::" + str(key) + "\n")
        tmpSTCNT=str(value.STCNT)
        #print "山頭數量:"+tmpSTCNT
        f.write("山頭數量:" + tmpSTCNT + "\n")

    f.write("------------------------山頭點個數統計結束------------\n")
    f.close()



    arcpy.Intersect_analysis([stmData, gqData,sxData], insertGQShp, "ALL", "#", "INPUT")#山頭面資料和國情資料
    cursor=arcpy.da.SearchCursor(insertGQShp, ["NAME", 'CC', "[email protected]", "[email protected]"])

    for row in cursor:
        flagSXLB = row[0]
        flagSXLB = flagSXLB.encode("utf-8")
        flagCode = row[1]
        flagArea = row[2]

        caculateArea(exportData,flagSXLB,flagCode,flagArea,ZZTDListCode,LCFGListCode,FWJZQListCode,TLYDLListCode,
                     GZWListCode,RGWJDListCode,FMYLLDListCode,SYListCode)

    f=file(txtPath,"a+")# 追加的方式讀寫
    f.write("------------------------山頭面積統計開始------------\n")
    for key,value in exportData.items():
        kLandInfo=value
        f.write(str(key)+"種植土地面積為:" + str(kLandInfo.ZZTD) + "\n")
        f.write(str(key)+"林草覆蓋面積為:" + str(kLandInfo.LCFG) + "\n")
        f.write(str(key)+"房屋建築區面積為:" + str(kLandInfo.FWJZQ) + "\n")
        f.write(str(key)+"鐵路與道路面積為:" + str(kLandInfo.TLYDL) + "\n")
        f.write(str(key)+"構築物面積為:" + str(kLandInfo.GZW) + "\n")
        f.write(str(key)+"人工挖掘地面積為:" + str(kLandInfo.RGWJD) + "\n")
        f.write(str(key)+"荒漠與裸露地面積為:" + str(kLandInfo.FMYLLD) + "\n")
        f.write(str(key)+"水域面積為:" + str(kLandInfo.SY) + "\n")
        f.write("------------------------------------\n")
    f.write("------------------------山頭面積統計結束------------\n")
    f.close()

    workbook = xlwt.Workbook(encoding='utf-8')
    worksheet = workbook.add_sheet('My Sheet')
    #表頭部分
    worksheet.write(0,1,"山頭數".encode("utf-8"))
    worksheet.write(0,2,"山頭面積".encode("utf-8"))
    worksheet.write(0,3,"地類總面積".encode("utf-8"))
    worksheet.write(0,4,"種植土地".encode("utf-8"))
    worksheet.write(0,5,"林草覆蓋".encode("utf-8"))
    worksheet.write(0,6,"房屋建築區".encode("utf-8"))
    worksheet.write(0,7,"鐵路與道路".encode("utf-8"))
    worksheet.write(0,8,"構築物".encode("utf-8"))
    worksheet.write(0,9,"人工挖掘地".encode("utf-8"))
    worksheet.write(0,10,"荒漠與裸露地表".encode("utf-8"))
    worksheet.write(0,11,"水域".encode("utf-8"))

    index=1

    #最後的面積統計
    totalLandinfo=LandInfo()
    for key,value in exportData.items():
        mLandInfo=value
        strSXLB=str(key)
        strSXLB=strSXLB.encode("utf-8")
        insertExcelData(strSXLB,index,worksheet,mLandInfo)
        caculteFinalArea(mLandInfo,totalLandinfo)
        index=index+1

    #插入最後一行統計資料
    insertExcelData("合計",index,worksheet,totalLandinfo)

    #刪除xls根目錄下所有的檔案
    delDirAllFiles(xlsRootPath)
    workbook.save(xlsRootPath+'result.xls')

打包程式碼

# coding:gbk
import arcpy

from StatisticsAreaClass import do


class Toolbox(object):
    def __init__(self):
        """Define the toolbox (the name of the toolbox is the name of the
        .pyt file)."""
        self.label = "Toolbox"
        self.alias = ""

        # List of tool classes associated with this toolbox
        self.tools = [Tool]


class Tool(object):
    def __init__(self):
        """Define the tool (tool name is the name of the class)."""
        self.label = "xxxx大資料中心"
        self.description = "xxxx大資料中心"
        self.canRunInBackground = False

    def getParameterInfo(self):
        """Define parameter definitions"""
        stdPath = arcpy.Parameter(
            displayName="山頂點資料",
            name="stdPath",
            datatype="GPFeatureLayer",
            parameterType="Required",
            direction="Input"
        )
        stfwxPath = arcpy.Parameter(
            displayName="山頭範圍線",
            name="stfwxPath",
            datatype="GPFeatureLayer",
            parameterType="Required",
            direction="Input"
        )
        sxPath = arcpy.Parameter(
            displayName="山形圖層",
            name="sxPath",
            datatype="GPFeatureLayer",
            parameterType="Required",
            direction="Input"
        )
        gqPath = arcpy.Parameter(
            displayName="國情要素圖層",
            name="gqPath",
            datatype="GPFeatureLayer",
            parameterType="Required",
            direction="Input"
        )

        resultFolder = arcpy.Parameter(
            displayName="最終處理結果",
            name="resultFolder",
            datatype="Folder",
            parameterType="Required",
            direction="Input"
        )
        #params = None
        params = [stdPath, stfwxPath, sxPath, gqPath, resultFolder]
        return params

    def isLicensed(self):
        """Set whether tool is licensed to execute."""
        return True

    def updateParameters(self, parameters):
        """Modify the values and properties of parameters before internal
        validation is performed.  This method is called whenever a parameter
        has been changed."""
        return

    def updateMessages(self, parameters):
        """Modify the messages created by internal validation for each tool
        parameter.  This method is called after internal validation."""
        return

    def execute(self, parameters, messages):
        """The source code of the tool."""
        stdPath = parameters[0].valueAsText
        stfwxPath = parameters[1].valueAsText
        sxPath = parameters[2].valueAsText
        gqPath = parameters[3].valueAsText
        resultFolder = parameters[4].valueAsText
        do(stdPath, stfwxPath, sxPath, gqPath, resultFolder)
        return

                                                                            更多內容,請關注公眾號