ArcGIS和Fragstats的脚本化调用 ------以ArcPy和命令行的方式
目錄
前言
一、思路
二、源碼
1.工作文件夾結構
2.掩膜提取
3.fragstas中的計算
三、問題
總結
前言
? ? ? ?自從大一參加研究開始，到現(xiàn)在大三，終于有能力將兩年前的問題（如下）解決掉，雖然可能解決方案不是完全的普適。你是否也身陷于數(shù)據處理過程中與軟件UI交互的囹圄？如何利用編程解放自己的雙手和雙眼？我的問題中涉及到兩個應用程序ArcGIS和fragstats，前者我使用其提供的ArcPy庫解決，后者我在其官網上的R包使用教程中發(fā)現(xiàn)了端倪，找到了其命令行方式調用的格式及參數(shù)意義，最終使用python將整個流程串了起來。
?之前的問題：
最近在跟學長做項目，很多事情都是機械化的操作，但涉及到多個應用程序，tedious?and?repetitive。于是不免yy，要是能把這個操作流程給編定該多好，我就輸輸參數(shù)，躺著就是了。So，請問現(xiàn)在有沒有可以實現(xiàn)如此功能的操作？如果沒有，有沒有有可行性的方向？
大體流程是先在arcgis中將數(shù)據進行掩膜提取，然后將得到的tif用py跑取路徑，導入fragstats進行指標計算，然后再導入arcgis連接表并輸出。
?
一、思路
? ? ? ? 這個問題的實質是怎樣用一門編程語言來使用一個程序，請注意這里說的是“使用”，即完全替代或者說模擬人與UI的交互。雖說替代這個詞用來總覺得有些退步的味道，畢竟計算機發(fā)展的歷史上，應該是先人們費了不少心血才用GUI替代了對一般用戶不友好的命令行交互操作。但當涉及到軟件設計者未能提供的功能時，我們還是要回到更親近計算機的語言來指揮軟件要干什么。
????????最初我的想法是：應用程序有其內在的結構與邏輯，對外表現(xiàn)為一個整體，協(xié)調應用程序間的關系，應該由操作系統(tǒng)來實現(xiàn)，于是我想到了shell。如果我能夠使用命令行來操作各個程序，以命令行的語言來寫一個腳本，順序化調用各個應用程序，并讓其在循環(huán)中自動生成相應參數(shù)，那么就能實現(xiàn)“躺著”了。可我并不知道是不是每一個程序都可以用命令行加參數(shù)的形式來調用，這一塊的原理及實現(xiàn)我并不清楚（如有了解的，望不吝賜教，即“是不是所有程序都可以用命令行加參數(shù)的形式來調用，以代替它在GUI上的交互”，如果一個程序可以采用這種方式調用，那么這塊在編寫出程序的過程中是怎么實現(xiàn)的呢？）。最終問題的解決過程中，也隱含地涉及到了這一方式。
二、源碼
1.工作文件夾結構
? ? ? ? 為了方便想理解源碼的同學，把文件夾結構貼在下面：
?
NEW
? ? ? ? ?2002
? ? ? ? ?2003
? ? ? ? ?...
? ? ? ? ?2015? ? ? ? ? ? ? # 存放各年份掩膜提取的結果
? ? ? ? ?CSV? ? ? ? ? ? ? # 存放fragstats計算出來的最終結果
? ? ? ? ?Feature? ? ? ? ?# 連接各年份表后分別導出至此文件夾
Raster? ? ? ? ? ?# 存放最終導出的柵格結果
? ? ? ? ?tes? ? ? ? ? ? ? ? ?# 測試性能確定掩膜提取時的進程數(shù)
2.掩膜提取
? ? ? ? 參考了不少人的博客及文檔如下：
使用ArcPy的環(huán)境配置：
[1]?Python2.7.14配置ArcGIS10.2自帶的arcpy環(huán)境
使用的具體函數(shù)調用格式：
[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): ?# 文件夾個數(shù)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): ?# 文件夾個數(shù)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"# 進程函數(shù) 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序號, ? ? ? ?同一數(shù)字序號的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的名稱#同一數(shù)字序號的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) ?# 創(chuàng)建進程池，進程數(shù)目為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)#性能測試已確定最優(yōu)進程數(shù)，或者直接設置成物理內核數(shù)也可以 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() ?# 邏輯內核數(shù)print psutil.cpu_count(False) ?# 物理內核數(shù)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')
?
? ? ? ? 后續(xù)的操作等有空了再發(fā)。
三、問題
? ? ? ? 不管是哪個問題您知道答案，亦或者是能給我一個方向，都懇請您能不吝賜教，我一定洗耳恭聽。
總結
? ? ? ? 一開始是抱著偷懶的心態(tài)去搞，結果被從來沒接觸過的多進程卡了兩周，好在最終結果收獲頗豐：不僅僅節(jié)省了人力，出乎意料的竟然還及大幅度地提高了運行的速度，可以說是意外之喜，只是可惜我的暫停功能可能用不上了hh。這也令我唏噓不已，所謂工欲善其事必先利其器，我們這兩年就如同在拿一塊生銹發(fā)鈍的刀來砍一塊肉，費了老鼻子勁砍完三分之二了，結果突發(fā)奇想磨了磨刀，一刀就把剩下的三分之一砍完了，結果是哭也不是，笑也不是。
總結
以上是生活随笔為你收集整理的ArcGIS和Fragstats的脚本化调用 ------以ArcPy和命令行的方式的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 用极域课堂管理系统软件批量格式化D盘
- 下一篇： WordPress 最新RiPro9.0