目录
前言
自从大一参加研究开始，到现在大三，终于有能力将两年前的问题（如下）解决掉，虽然可能解决方案不是完全的普适。你是否也身陷于数据处理过程中与软件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例子
多进程方面的参考：
[5] Python 进阶必备：进程模块 multiprocessing
另外还要特别感谢许向武老师和涂建光老师的解惑
代码如下：
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): # 文件夹个数72
if i == 0:
continue
elif 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): # 文件夹个数124
if i == 0:
continue
elif 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 -1
else:
return -1
# 供子进程调用保存，生成或修改日志文件
def check_out(path, year, filename):
if os.access(path + "/run.log", os.F_OK): # 存在之前的日志文件
lines = []
check = 0
with 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 string
clipf.write(string)
check = 1
else:
clipf.write(line_t)
if check == 0:
string = str(year) + "," + filename + '\n'
print string
clipf.write(string)
else: # 不存在
clipf = open(path + "/run.log", 'a+')
string = str(year) + "," + filename + '\n'
print string
clipf.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 # 设置工作空间，以便调用listfeature
shplist_1 = arcpy.ListFeatureClasses() # 读取工作空间内的要素类文件名（unicode）
#上次跑到的shp序号, 同一数字序号的shp可能存在两份，其中一份末尾会有一个'_' e.g."1.shp", "1_.shp"
lastfilename = readlog(path, year)
if lastfilename == -1:
print "no history of " + path
elif 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):
continue
else:
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):
continue
else:
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) # 创建进程池，进程数目为i
if 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 = 1
print "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，不是普通的int
print 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 os
root = "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 " + name
more = ",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): # 文件夹个数124
if i == 0:
continue
elif 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和命令行的方式分享到此就结束了，感谢您的阅读，如果确实帮到您，您可以动动手指转发给其他人。
版权声明：本文内容由互联网用户自发贡献，该文观点仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 举报，一经查实，本站将立刻删除。
如需转载请保留出处：https://bianchenghao.cn/24659.html
