ArcGIS和Fragstats的脚本化调用 ------以ArcPy和命令行的方式
目录
前言
一、思路
二、源码
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')
后续的操作等有空了再发。
三、问题
不管是哪个问题您知道答案，亦或者是能给我一个方向，都恳请您能不吝赐教，我一定洗耳恭听。
总结
一开始是抱着偷懒的心态去搞，结果被从来没接触过的多进程卡了两周，好在最终结果收获颇丰：不仅仅节省了人力，出乎意料的竟然还及大幅度地提高了运行的速度，可以说是意外之喜，只是可惜我的暂停功能可能用不上了hh。这也令我唏嘘不已，所谓工欲善其事必先利其器，我们这两年就如同在拿一块生锈发钝的刀来砍一块肉，费了老鼻子劲砍完三分之二了，结果突发奇想磨了磨刀，一刀就把剩下的三分之一砍完了，结果是哭也不是，笑也不是。
总结
以上是生活随笔为你收集整理的ArcGIS和Fragstats的脚本化调用 ------以ArcPy和命令行的方式的全部内容，希望文章能够帮你解决所遇到的问题。
- 上一篇： 用极域课堂管理系统软件批量格式化D盘
- 下一篇： WordPress 最新RiPro9.0