前言

自从大一参加研究开始，到现在大三，终于有能力将两年前的问题（如下）解决掉，虽然可能解决方案不是完全的普适。你是否也身陷于数据处理过程中与软件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的模型构建器构建的可迭代模型，由于只支持一层迭代，所以需要人工进行参数配置，我一开始想要写代码解决的正是这一步。可当上述的代码写出来我发现相比于模型，ArcPy的方式效率高了好多倍，原本需要断断续续跑两周的工作，现在只需1.64h。那么这是为什么呢？二者本质上应该还都是调用的ArcGIS的C源码吧？怎么会有如此大的不同？观察发现，模型一开始跑的也很快，但跑到后期在任务管理器中可以看到其cpu、内存占用率都接近于0，跑取一个shp的时间也从开始的一秒不到跌落到1~2min。这到底是为什么呢？

不管是哪个问题您知道答案，亦或者是能给我一个方向，都恳请您能不吝赐教，我一定洗耳恭听。

总结

一开始是抱着偷懒的心态去搞，结果被从来没接触过的多进程卡了两周，好在最终结果收获颇丰：不仅仅节省了人力，出乎意料的竟然还及大幅度地提高了运行的速度，可以说是意外之喜，只是可惜我的暂停功能可能用不上了hh。这也令我唏嘘不已，所谓工欲善其事必先利其器，我们这两年就如同在拿一块生锈发钝的刀来砍一块肉，费了老鼻子劲砍完三分之二了，结果突发奇想磨了磨刀，一刀就把剩下的三分之一砍完了，结果是哭也不是，笑也不是。

ArcGIS和Fragstats的脚本化调用 ------以ArcPy和命令行的方式相关推荐

android 能调用gcc_如何在命令行下使用Android NDK交叉编译工具
我们知道,在Linux下可以使用gcc来把一份C代码编译成为Linux上的可执行程序, 如: $ gcc -o main.out main.c 而Android平台提供了NDK工具包来交叉编译可以运行 ...
android java调用参数,如何从命令行调用Android JNI函数并传递Java对象参数
一.前言当我们对某个使用原生库(native library)的恶意软件或者应用进行分析或渗透测试时,如果能够对库函数进行隔离和执行是再好不过的事情,这样做我们就可以使用其自身的代码来调试对抗恶意软 ...
Mac电脑调用自带的命令行窗口
前言提示:我们都熟悉Windows电脑通过win+r即可以快速启动命令行窗口,但是Mac电脑中如何实现调用命令行窗口呐? 一.快捷键:command + 空格二.敲入te,terinmal的缩写即 ...
nutch java_Nutch：用Java调用，而不是命令行？
如果您查看bin/nutch脚本,您将看到它调用与您的命令对应的Java类: # figure out which class to run if [ "$COMMAND" = & ...
wkhtmltopdf使用方法，页脚加页码 pdftk合并pdf命令行操作
wkhtmltopdf pdftk使用方法 htmltopdf 安装配置路径 htmltopdf 打开官网:https://wkhtmltopdf.org/downloads.html 安装配置路径 ...
VBA调用bat，doc 命令行窗口关闭之后，VBA代码再继续执行
VBA调用doc命令,doc窗口关闭之后,继续执行代码 Option Explicit'Docワィンドワ閉じるした後.後続けの処理実施 Public Declare Function WaitForS ...
Ubuntu 调用查看USB摄像头命令行介绍
1.查看当前插入的摄像头 ls -ltrh /dev/video* 输出 $ ls -ltrh /dev/video*crw-rw----+ 1 root video 81, 1 12月 29 10: ...
第17章　脚本化CSS
第17章脚本化CSS CSS脚本化是网页交互效果的技术基础,使用CSS和JavaScript可以设计网页动画.利用脚本化CSS可以动态地改变HTML属性,如字体颜色.字体大小等,还可以用它设置和改变 ...
python调用php命令行,python调用php函数 python怎样调用php文件中的函数详解
前言 python调用php代码实现思路:php文件可通过在terminal中使用php命令行进行调用,因此可使用python开启子进程执行命令行代码.函数所需的参数可通过命令行传递. 测试环境 1. ...
打开高效文本编辑之门_调用Linux的sed命令
Linux sed命令执行方式汇总案例声明与简介 sed:Stream Editor文本流编辑,sed是一个"非交互式的"面向字符流的编辑器.Sed的命令执行主要介绍如何引用se ...

ArcGIS和Fragstats的脚本化调用 ------以ArcPy和命令行的方式

前言