求助arcGIS中Python编程批处理问题——坐标XYZM转shape图层

壮志凌云MAN · 发表于 2016-9-12 14:42

如题，朋友做关于流行病学项目研究，使用到了全国830多个全国气象监测站点的日监测数据，有7~8年的数据量，原始数据为数值型，包括XYZ,及一个气象指标。先需要使用监测数据做全国范围的插值提取，思路：使用原始数据转shape点，使用插值方法进行插值，然后使用全国主要城市坐标提取各城市的气象数据。
由于数据文件多（365*10），所以只能求助使用批处理方法进行，前期使用模型构建器进行了实验。具体如下：
--------------------------------模型构建器方法--------------------------------------------------
使用的toolbox有【sample】create feature from text file\【插值分析】反距离权重法等，关于插值方法的选择不做深入讨论。在数据处理流程中的插值分析和插值提取过程是简单顺利的，使用模型构建起可进行批处理【批处理使用迭代，具体参考arcgis中帮助文档即可】，该方法的前面原始数据XYZM转点shape图层不完善【使用9.3版的sample-create feature from text file】工具，该方法输入数据文件为text，且对数据格式有严格要求，而且转换的数据只能包括XY坐标，高程Z和属性M均不能插入，这一点和单个坐标文件转换完全不同。使用模型构建器方法，前部分处理不能完成，所以想到使用Python来处理。
输入text文档的格式如图
txt文件格式.JPG

-------------------------------python方法--------------------------------------------------------
使用Python处理要解决2个问题：批处理，以及ZM字段的插入。
------------------------------------------------------------------------------------------------------------------------
【1】批处理：没有使用Python做个批处理{之前一直使用R}，所以还得求助。
在小木虫查到一篇关于批处理的帖子：http://muchong.com/html/201505/8981718.html
查到的资料如下【类似的批处理，文件名是有日期构成】

#批处理参考模板
import arcpy
def transform(years):
filename_pattern = "F:/data/%d/NDVI_CN_%d_%02d.tif"
new_filename_pattern = "F:/newdata/%d/NDVI_CN_%d_%02d.tif"
for year in years:
for month in range(1, 13):
filename = filename_pattern % (year, year%100, month)
new_filename = new_filename_pattern % (year, year%100, month)
# 转换文件filename
arcpy.Resample_management(filename, new_filename, "1000 1000", "NEAREST")
if __name__=='__main__':
import sys
start, end = map(int, sys.argv[1:])
transform(range(start, end+1))

复制代码

我要处理的文件是

------------------------------------------------------------------------------------------------------------------------------------
【2】ZM字段的插入
从找到的资料来看，要在坐标转shape点图层中实现属性字段的插入，有2类方法：修改sample中create feature from text file脚本 -/-使用【数据管理】/【要素类】/【创建要素类】，然后将坐标插入已建的要素类图层中。
-------------------------------------------------------------------------------------
【2-1】create feature from text file方法（部分脚本代码如下）中对ZM字段也有处理，但是在实际处理中却未把ZM插入，从结果的属性表中为发现，而ZM在本数据中是最为重要的。

def createPoint(point, geometry):
try:
point.id = geometry[0]
point.x = geometry[1]
point.y = geometry[2]
# When empty values are written out from pyWriteGeomToTextFile, they come as 1.#QNAN
# Additionally, the user need not supply these values, so if they aren't in the list don't add them
if len(geometry) > 3:
if geometry[3].lower().find("nan") == -1: point.z = geometry[3]
if len(geometry) > 4:
if geometry[4].lower().find("nan") == -1: point.m = geometry[4]
return point
except:
raise Exception, msgErrorCreatingPoint

复制代码

---------------------------------------------------------------
【2-2】使用【数据管理】/【要素类】/【创建要素类】，然后将坐标插入已建的要素类图层中。
从网上找的一段代码，可惜没看到输入输出函数啊。。。求牛人指教！代码贴下面了。

#参考来源：http://www.faqssys.info/creating-shp-from-txt-files-using-arcpy/#comments
import os, sys, arcpy
InFile = sys.argv[1]
OutShp = sys.argv[2]
EPSG = sys.argv[3]
try:
SR = arcpy.SpatialReference(int(EPSG))
# create a spatial reference from EPSG definition
except:
arcpy.AddWarning(arcpy.GetMessages())
arcpy.AddError("Unable to create spatial reference with code %s" % EPSG)
sys.exit(0)
#break up OutShp to folder and name
OutFolder = os.path.dirname(OutShp)
OutName = os.path.basename(OutShp)
# Create the output feature class/创建函数
arcpy.CreateFeatureclass_management(OutFolder,OutName,"POINT",has_z="ENABLED",spatial_reference=SR)
# add fields to store the values
arcpy.AddField_management(OutShp,"PntID","INTEGER")
arcpy.AddField_management(OutShp,"Xcoord","DOUBLE")
arcpy.AddField_management(OutShp,"Ycoord","DOUBLE")
arcpy.AddField_management(OutShp,"Zcoord","DOUBLE")
with open(InFile,'r') as srcFile:
with arcpy.da.InsertCursor(OutShp,["SHAPE@","PntID","Xcoord","Ycoord","Zcoord"]) as InsCur:
for fileLine in srcFile:
# split the line up into a list
# 1 337172.8 4278952.4 85.7334909091 becomes [1,337172.8,4278952.4,85.7334909091]
# try a few common delimiters: tab, space, comma, pipe../确定数据间隔符类型：tab、空格、逗号
lSplit = fileLine.split("t")
if len(lSplit) == 1:
lSplit = fileLine.split(" ")
if len(lSplit) == 1:
lSplit = fileLine.split(",")
if len(lSplit) == 1:
lSplit = fileLine.split("|")
if len(lSplit) > 1:
# more than just one word on the line
pointsOK = True
try:
ID = int(lSplit[0])
Xcoord = float(lSplit[1])
Ycoord = float(lSplit[2])
Zcoord = float(lSplit[3])
except:
arcpy.AddWarning("Unable to translate points")
pointsOK = False
if pointsOK:
# create a point geometry from the 3 coordinates
newGeom = arcpy.PointGeometry(arcpy.Point(Xcoord,Ycoord,Zcoord))
newGeom.spatial_reference = SR # set spatial reference
# you could project here using newGeom.projectAs(SpatialRef)
# if you needed them in a different spatial reference
InsCur.insertRow([newGeom,ID,Xcoord,Ycoord,Zcoord])# insert this point into the feature class

复制代码

--------------------------------------------------------------------------------------------------------------------------
总结下，如上问题，主要是解决坐标转点shape图层的批处理问题，涉及2点。希望各位好友能给出解答，不胜感激！

后话：关键是我现在研三，深入学习Python也需要些时日，而朋友的这个任务也是比较紧。之前对arcgis的批处理不清楚，没想到前面环节卡住了。。。现在要是等我自学段时间又不允许，手头杂事也多。求助各位啦！