目錄

前言

一、思路

二、源碼

1.工作文件夾結構

2.掩膜提取

3.fragstas中的計算

三、問題

總結

---

前言

? ? ? ?自從大一參加研究開始，到現在大三，終于有能力將兩年前的問題（如下）解決掉，雖然可能解決方案不是完全的普適。你是否也身陷于數據處理過程中與軟件UI交互的囹圄？如何利用編程解放自己的雙手和雙眼？我的問題中涉及到兩個應用程序ArcGIS和fragstats，前者我使用其提供的ArcPy庫解決，后者我在其官網上的R包使用教程中發現了端倪，找到了其命令行方式調用的格式及參數意義，最終使用python將整個流程串了起來。

?之前的問題：

最近在跟學長做項目，很多事情都是機械化的操作，但涉及到多個應用程序，tedious?and?repetitive。于是不免yy，要是能把這個操作流程給編定該多好，我就輸輸參數，躺著就是了。So，請問現在有沒有可以實現如此功能的操作？如果沒有，有沒有有可行性的方向？

大體流程是先在arcgis中將數據進行掩膜提取，然后將得到的tif用py跑取路徑，導入fragstats進行指標計算，然后再導入arcgis連接表并輸出。

---

?

一、思路

? ? ? ? 這個問題的實質是怎樣用一門編程語言來使用一個程序，請注意這里說的是“使用”，即完全替代或者說模擬人與UI的交互。雖說替代這個詞用來總覺得有些退步的味道，畢竟計算機發展的歷史上，應該是先人們費了不少心血才用GUI替代了對一般用戶不友好的命令行交互操作。但當涉及到軟件設計者未能提供的功能時，我們還是要回到更親近計算機的語言來指揮軟件要干什么。

????????最初我的想法是：應用程序有其內在的結構與邏輯，對外表現為一個整體，協調應用程序間的關系，應該由操作系統來實現，于是我想到了shell。如果我能夠使用命令行來操作各個程序，以命令行的語言來寫一個腳本，順序化調用各個應用程序，并讓其在循環中自動生成相應參數，那么就能實現“躺著”了?？晌也⒉恢朗遣皇敲恳粋€程序都可以用命令行加參數的形式來調用，這一塊的原理及實現我并不清楚（如有了解的，望不吝賜教，即“是不是所有程序都可以用命令行加參數的形式來調用，以代替它在GUI上的交互”，如果一個程序可以采用這種方式調用，那么這塊在編寫出程序的過程中是怎么實現的呢？）。最終問題的解決過程中，也隱含地涉及到了這一方式。

二、源碼

1.工作文件夾結構

? ? ? ? 為了方便想理解源碼的同學，把文件夾結構貼在下面：

?

NEW

? ? ? ? ?2002

? ? ? ? ?2003

? ? ? ? ?...

? ? ? ? ?2015? ? ? ? ? ? ? # 存放各年份掩膜提取的結果

? ? ? ? ?CSV? ? ? ? ? ? ? # 存放fragstats計算出來的最終結果

? ? ? ? ?Feature? ? ? ? ?# 連接各年份表后分別導出至此文件夾

Raster? ? ? ? ? ?# 存放最終導出的柵格結果

? ? ? ? ?tes? ? ? ? ? ? ? ? ?# 測試性能確定掩膜提取時的進程數

---

2.掩膜提取

? ? ? ? 參考了不少人的博客及文檔如下：

使用ArcPy的環境配置：

[1]?Python2.7.14配置ArcGIS10.2自帶的arcpy環境

使用的具體函數調用格式：

[2]?按掩膜提取 (Spatial Analyst)? ? ?(官方文檔)

[3]?兩個簡單的arcpy例子

[4]?ArcGIS Python編程案例(0)-前言

多進程方面的參考：

[5]?Python 進階必備：進程模塊 multiprocessing

[6]?Python的多線程和多進程——從一個爬蟲任務談起

另外還要特別感謝許向武老師和涂建光老師的解惑

? ? ? ? 代碼如下：

import os import multiprocessing as mp import arcpy import psutil import time# 準備工作，生成輸出路徑文件夾 def gera_folders(paths):for path in paths:os.mkdir("C:\\research\\NEW\\" + path[10:])# 測試用 # 生成包含測試路徑的list def gera_test_path():paths = []for i in range(73): ?# 文件夾個數72if i == 0:continueelif i < 10:path = "E:/Cgrid_/tes/Folder 0" + str(i)paths.append(path)else:path = "E:/Cgrid_/tes/Folder " + str(i)paths.append(path)return paths# 子進程用，以供迭代 # 生成形如“E:/Cgrid_/Folder 01” ~ “E:/Cgrid_/Folder 124”的文件夾路徑 def gera_folder_path():paths = []for i in range(125): ?# 文件夾個數124if i == 0:continueelif i < 10:path = "E:/Cgrid_/Folder 0" + str(i)paths.append(path)else:path = "E:/Cgrid_/Folder " + str(i)paths.append(path)return paths# 暫停功能 #讀取之前跑到的shp名，若本年沒跑過則返回-1 def readlog(path, year):if os.access(path + "/run.log", os.F_OK):with open(path + "/run.log", 'r') as f:for line_t in f:print int(line_t[0:4])if int(line_t[0:4]) == year:print "Last time run shp " + line_t[5:-1]return line_t[5:-1]else:return -1else:return -1# 供子進程調用保存，生成或修改日志文件 def check_out(path, year, filename):if os.access(path + "/run.log", os.F_OK): ?# 存在之前的日志文件lines = []check = 0with open(path + "/run.log", 'r') as f:for line_t in f:lines.append(line_t)clipf = open(path + "/run.log", 'w+')for line_t in lines:if int(line_t[0:4]) == year:string = str(year) + "," + filename + '\n'print stringclipf.write(string)check = 1else:clipf.write(line_t)if check == 0:string = str(year) + "," + filename + '\n'print stringclipf.write(string)else: ?# 不存在clipf = open(path + "/run.log", 'a+')string = str(year) + "," + filename + '\n'print stringclipf.write(string)print path + " has been saved!\n"# 進程函數 def single(path, year, flag): ?# path 這個進程要跑的文件夾(124中的一個）， year本次工作的被裁柵格年份（2002~2015）# file_lists = []# file_lists = get_my_shp_paths(path)[0]arcpy.CheckOutExtension("spatial") ?# 權限檢查arcpy.gp.overwriteOutput = 1 ?# 裁出來的tif重復可覆蓋raster = "C:\\research\\ESACCI-LC-L4-LCCS-Map-300m-P1Y-" + str(year) + "-v2.0.7.tif" ?# 被裁的柵格影像所在目錄arcpy.env.workspace = path ?# 設置工作空間，以便調用listfeatureshplist_1 = arcpy.ListFeatureClasses() ?# 讀取工作空間內的要素類文件名（unicode）#上次跑到的shp序號, ? ? ? ?同一數字序號的shp可能存在兩份，其中一份末尾會有一個'_' ? e.g."1.shp", "1_.shp"lastfilename = readlog(path, year)if lastfilename == -1:print "no history of " + pathelif lastfilename[-1] == '_':lastfilename = lastfilename[0:-1]for Inputfeature in shplist_1:(filename, extension) = os.path.splitext(Inputfeature.encode('utf8')) ?# 獲取shp的名稱#同一數字序號的shp可能存在兩份，其中一份末尾會有一個'_' ? e.g."1.shp", "1_.shp"if filename[-1] != '_':if int(filename) < int(lastfilename):continueelse:out = "C:\\research\\NEW\\" + str(year) + "\\" + path[10:] + "\\" + filename ?# 將shp的名稱作為tif輸出時的名稱 ? ?str(year)******************* path[14:]if flag.value == 0:# print("OK")arcpy.gp.ExtractByMask_sa(raster, Inputfeature, out + ".tif")# outExtractByMask = ExtractByMask(raster, Inputfeature) ? ? ? ? ? ? ? ? ?#另一種調用方式 需要from arcpy.sa import *# outExtractByMask.save(out + ".tif")elif flag.value == 1:# print("not ok")check_out(path, year, filename)else:filename = filename[0:-1]if int(filename) < int(lastfilename):continueelse:out = "C:\\research\\NEW\\" + str(year) + "\\" + path[10:] + "\\" + filename + '_' ?# 將shp的名稱作為tif輸出時的名稱 ? ?str(year)******************* path[14:]if flag.value == 0:# print("OK")arcpy.gp.ExtractByMask_sa(raster, Inputfeature, out + ".tif")elif flag.value == 1:# print("not ok")check_out(path, year, filename)if flag.value == 0:print path + " has finished!\n"#掩膜提取主進程state = 1為正常運行，state = 0 為性能測試時使用 def Run(i,year,state = 1):mpp = mp.Pool(processes=i) ?# 創建進程池，進程數目為iif state == 1:file_list_all = gera_folder_path()elif state == 0:file_list_all = gera_test_path()else:print "State error!"for filelist in file_list_all:mpp.apply_async(single, args=(filelist, year, flag)) ?# 非阻塞提交新任務mpp.close() ?# 關閉進程池，意味著不再接受新的任務while (True):num = input("Input 1 to save and quit\n")try:if num == 1:flag.value = 1print "trying to save progress\n"else:print "illegal input\n"except Exception as ex:print ("表達式為空，請檢查/loadingtimeout")mpp.join()print ("whole year finished with %d", i)#性能測試已確定最優進程數，或者直接設置成物理內核數也可以 def testE(year):for i in range(mp.cpu_count()):if i > 3:begin_time = time.time()Run (i, year, 0)end_time = time.time()print (end_time - begin_time)if __name__ == '__main__':years = [2002，2003，2004，2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2014, 2015]for year in years:begin_time = time.time()# 供主進程傳遞保存指令給子進程，0為正常運行，1為開始保存flag = mp.Manager().Value('i', 0) ?# flag類型是ctypes.c_long，不是普通的intprint mp.cpu_count() ?# 邏輯內核數print psutil.cpu_count(False) ?# 物理內核數print "will be running "+str(year)+ "\n"Run(psutil.cpu_count(False),year)end_time = time.time()print(end_time - begin_time)

?

3.fragstas中的計算

? ? ? ? 這一塊的思路來自于對fragstats官網上的R包教程的解析：FRAGSTATS: Spatial Pattern Analysis Program for Categorical Maps Downloads，在其中用R語言調用了如下命令行：

#execute fragstats
system('frg -m fragmodelR.fca -b geotiffbatch.fbt -o c:\\work\\fragstats\\tutorial\\tutorial_5\\fragout')

# -*- coding: utf-8 -*- import osroot = "C:/research/NEW/2002"def gera_pathfbt(year):names = []for i in range(9):name = str(0) + str(i + 1)names.append(name)for i in range(9, 124):name = str(i + 1)names.append(name)for name in names:clip = "C:/research/NEW/" + str(year) + "/Folder " + namemore = ",x,999,x,x,1,x,IDF_GeoTIFF"pathlist = os.listdir(clip)clipf = open(clip + "/_" + "1" + ".fbt", 'w+')for line in pathlist:(filename, extension) = os.path.splitext(line)if (extension == ".tif"):clipf.write(clip + '/' + line + more + '\n')def Frg(year):root = "C:/research/NEW/"+str(year)paths = []for i in range(125): # 文件夾個數124if i == 0:continueelif i < 10:path = root + "/Folder 0" + str(i)paths.append(path)else:path = root + "/Folder " + str(i)paths.append(path)for path in paths:os.chdir(path)out = path +"/fragout"fbt = "_1.fbt"fca = "C:/research/unnamed8.fca"os.system('frg -m '+ fca +' -b '+ fbt + ' -o ' + "\"" + out + "\"")if __name__ == '__main__':years = [2008]for year in years:print "will be running "+str(year)+ "\n"gera_pathfbt(year)Frg(year)

?

? ? ? ? 后續的操作等有空了再發。

三、問題

是不是所有程序都可以用命令行加參數的形式來調用，以代替它在GUI上的交互？如果一個程序可以采用這種方式調用，那么這塊在編寫出程序的過程中是怎么實現的呢？

之前跑掩膜提取的時候，我們使用的是ArcGIS的模型構建器構建的可迭代模型，由于只支持一層迭代，所以需要人工進行參數配置，我一開始想要寫代碼解決的正是這一步?？僧斏鲜龅拇a寫出來我發現相比于模型，ArcPy的方式效率高了好多倍，原本需要斷斷續續跑兩周的工作，現在只需1.64h。那么這是為什么呢？二者本質上應該還都是調用的ArcGIS的C源碼吧？怎么會有如此大的不同？觀察發現，模型一開始跑的也很快，但跑到后期在任務管理器中可以看到其cpu、內存占用率都接近于0，跑取一個shp的時間也從開始的一秒不到跌落到1~2min。這到底是為什么呢？

? ? ? ? 不管是哪個問題您知道答案，亦或者是能給我一個方向，都懇請您能不吝賜教，我一定洗耳恭聽。

---

總結

? ? ? ? 一開始是抱著偷懶的心態去搞，結果被從來沒接觸過的多進程卡了兩周，好在最終結果收獲頗豐：不僅僅節省了人力，出乎意料的竟然還及大幅度地提高了運行的速度，可以說是意外之喜，只是可惜我的暫停功能可能用不上了hh。這也令我唏噓不已，所謂工欲善其事必先利其器，我們這兩年就如同在拿一塊生銹發鈍的刀來砍一塊肉，費了老鼻子勁砍完三分之二了，結果突發奇想磨了磨刀，一刀就把剩下的三分之一砍完了，結果是哭也不是，笑也不是。

總結

以上是生活随笔為你收集整理的ArcGIS和Fragstats的脚本化调用 ------以ArcPy和命令行的方式的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。