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

相似文件 换一批

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

如题，朋友做关于流行病学项目研究，使用到了全国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文档的格式如图

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

1. #批处理参考模板
   
2. import arcpy
   
3. def transform(years):
   
4.     filename_pattern = "F:/data/%d/NDVI_CN_%d_%02d.tif"
   
5.     new_filename_pattern = "F:/newdata/%d/NDVI_CN_%d_%02d.tif"
   
6.     for year in years:
   
7.        for month in range(1, 13):
   
8.            filename = filename_pattern % (year, year%100, month)
   
9.            new_filename = new_filename_pattern % (year, year%100, month)
   
10.            # 转换文件filename
    
11.            arcpy.Resample_management(filename, new_filename, "1000 1000", "NEAREST")
    
12. if __name__=='__main__':
    
13.     import sys
    
14.     start, end = map(int, sys.argv[1:])
    
15.     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在本数据中是最为重要的。

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

复制代码

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

1. #参考来源：http://www.faqssys.info/creating-shp-from-txt-files-using-arcpy/#comments
   
2. 
   
3. import os, sys, arcpy
   
4. 
   
5. InFile = sys.argv[1]
   
6. OutShp = sys.argv[2]
   
7. EPSG   = sys.argv[3]
   
8. 
   
9. try:
   
10.     SR = arcpy.SpatialReference(int(EPSG))
    
11.    
    
12.     # create a spatial reference from EPSG definition
    
13. 
    
14. except:
    
15.     arcpy.AddWarning(arcpy.GetMessages())
    
16.     arcpy.AddError("Unable to create spatial reference with code %s" % EPSG)
    
17.     sys.exit(0)
    
18. 
    
19. #break up OutShp to folder and name
    
20. OutFolder = os.path.dirname(OutShp)
    
21. OutName   = os.path.basename(OutShp)
    
22. 
    
23. # Create the output feature class/创建函数
    
24. arcpy.CreateFeatureclass_management(OutFolder,OutName,"POINT",has_z="ENABLED",spatial_reference=SR)
    
25. # add fields to store the values
    
26. arcpy.AddField_management(OutShp,"PntID","INTEGER")
    
27. arcpy.AddField_management(OutShp,"Xcoord","DOUBLE")
    
28. arcpy.AddField_management(OutShp,"Ycoord","DOUBLE")
    
29. arcpy.AddField_management(OutShp,"Zcoord","DOUBLE")
    
30. 
    
31. with open(InFile,'r') as srcFile:
    
32.     with arcpy.da.InsertCursor(OutShp,["SHAPE@","PntID","Xcoord","Ycoord","Zcoord"]) as InsCur:
    
33.         for fileLine in srcFile:
    
34.         
    
35.             # split the line up into a list
    
36.             # 1 337172.8 4278952.4 85.7334909091 becomes [1,337172.8,4278952.4,85.7334909091]
    
37.             # try a few common delimiters: tab, space, comma, pipe../确定数据间隔符类型：tab、空格、逗号
    
38. 
    
39.             lSplit = fileLine.split("t")
    
40.             if len(lSplit) == 1:
    
41.                 lSplit = fileLine.split(" ")
    
42.             if len(lSplit) == 1:
    
43.                 lSplit = fileLine.split(",")
    
44.             if len(lSplit) == 1:
    
45.                 lSplit = fileLine.split("|")
    
46. 
    
47.             if len(lSplit) > 1:
    
48.                 # more than just one word on the line
    
49.                 pointsOK = True
    
50.                 try:
    
51.                     ID     = int(lSplit[0])
    
52.                     Xcoord = float(lSplit[1])
    
53.                     Ycoord = float(lSplit[2])
    
54.                     Zcoord = float(lSplit[3])
    
55.                 except:
    
56.                     arcpy.AddWarning("Unable to translate points")
    
57.                     pointsOK = False
    
58. 
    
59.                 if pointsOK:
    
60.                     # create a point geometry from the 3 coordinates
    
61.                     newGeom  = arcpy.PointGeometry(arcpy.Point(Xcoord,Ycoord,Zcoord))
    
62.                     newGeom.spatial_reference = SR # set spatial reference
    
63.                     # you could project here using newGeom.projectAs(SpatialRef)
    
64.                     # if you needed them in a different spatial reference                        
    
65.                     InsCur.insertRow([newGeom,ID,Xcoord,Ycoord,Zcoord])# insert this point into the feature class

复制代码

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

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

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：Python编程 python ARCGIS Shape GIS 流行病学 朋友 模型 气象 项目

做一个excel 包含XYZ 然后倒入arcgis ---右键----display xy 就可以转成点 编什么程序

weiyuqiangweiyu 发表于 2016-9-12 19:10
做一个excel 包含XYZ 然后倒入arcgis ---右键----display xy 就可以转成点 编什么程序

上千个文件不能用那个方法，必须用编程啊！

Chengweing 发表于 2016-9-12 19:16
上千个文件不能用那个方法，必须用编程啊！

哦 我以为就一个excel表格呢  好厉害

不知楼主解决没有。没有的话可联系我有偿帮忙。