去年花了3天时间，学习了Python 的语法，除当时写了几个小游戏之外，一直没有用于实践解决实际问题。解决实际问题是学习编程的最好方法。学习Python最初的目的是用python写arcgis脚本，完成批处理。刚好，一个师妹最近要处理大批量的global30的数据，于是我试着用Python写了个批处理。做一个记录，记下写代码过程中遇到到问题和解决办法，方便以后查阅。
程序目标
将800多个分副的Global30数据中的人工建设用地提取出来，转成Shaplefile文件，再合成一个完整的shapefile文件。
用到的arcgis 模块有 Feature Clip, Reclassify, PolygontoRaster, Feature Merge.
程序实现
要素遍历
批处理程序第一步肯定是要获取需要处理的所有文件，也就是文件遍历。Arcpy站点包中，为Arcgis的Feature\Raster\Table等要素专门设置了函数，分别为arcpy.ListFeatureClasses()、arcpy.ListRasters()、arcpy.ListTable()等。返回的是文件名的list集合。
功能实现
ArcGis 中对每一个ArcToolBox 中的工具基本都有对应的Python语句，可以从帮助中查询。
FeatureClip
arcpy.Clip_management(precRaster, "#", outfile, shapeList[i], "0", "ClippingGeometry")
Reclassfiy
Reclassify(outfile, "Value", RemapRange([[0, 79, "NODATA"], [80, 80, 1], [81, 255, "NODATA"]]), "NODATA")
Reclassify 中提供两种重分类的赋值方法，意思范围RemapRange,另外是一对一赋值RemapValue.
Reclassify 如果没有栅格属性表，可能回报错。因此在Reclassify 前要创建一个栅格属性表。可用arcpy.BuildRasterAttributeTable_management(outfile, "Overwrite")创建。
PolygontoRaster
arcpy.RasterToPolygon_conversion(outReclass, outPolygons, "NO_SIMPLIFY", "Value")
Feature Merge
arcpy.Merge_management(featureList, "OutFiles.shp")
辅助函数
文件名解析
文件名解析也是批处理中重要的一个内容，可以方便地控制输出文件。可用OS 模块或者，用split分割输入文件名。再用 + 连接新的字符串，形成输出文件名。
调试
调试是写程序的重要工作，“行百里者半九十”，调试就是那最后的占一半的10percent。
Python 中调试的语句是：
try:
语句
except:
错误处理
可以提示错误在哪儿，也可以进行下一步操作。
完整代码
# -*- coding: utf-8 -*-
##Clip
import arcpy
from arcpy import env
from arcpy.sa import *
# arcpy.env.workspace = r"D:\JM\HelpOthers\global20\data"
arcpy.env.workspace = r"E:\goble"
rasterList = arcpy.ListRasters("*")
shapeList = arcpy.ListFeatureClasses("*")
# print rasterList, shapeList
n = len(rasterList)
m = len(shapeList)
print "共有" + str(n) + "副栅格"
error = open("E:\\goble\\error.txt", "w")
eundo = open("E:\\goble\\eundo.txt", "w")
enoarti = open("E:\\goble\\enoarti.txt", "w")
# outfile = "test_clip.tif"
if m == n:
for i in range(0, n, 1):
# 栅格裁剪 如果裁剪失败，提示错误但不终止重新
precRaster = rasterList[i]
# print precRaster, shapeList[i]
outfile = arcpy.env.workspace + "\\temp\\" + precRaster.split(".")[0] + "_clip.tif"
try:
print "clip Raster"
arcpy.Clip_management(precRaster, "#", outfile, shapeList[i], "0", "ClippingGeometry")
except:
print precRaster, shapeList[i],"Clip Failed"
print arcpy.GetMessages()
error.write(precRaster+ " " +shapeList[i]+" Clip Failed")
error.write(arcpy.GetMessages())
outfile = precRaster
pass
# 重分类
arcpy.CheckOutExtension("Spatial")
try:
print "Reclassify"
outReclass = Reclassify(outfile, "Value", RemapRange([[0, 79, "NODATA"], [80, 80, 1], [81, 255, "NODATA"]]), "NODATA")
outfile = arcpy.env.workspace + "\\temp\\" + precRaster.split(".")[0] + "_Reclass.tif"
outReclass.save(outfile)
except:
print precRaster, shapeList[i],"Reclassify Failed"
print arcpy.GetMessages()
error.write(precRaster+" "+ shapeList[i]+" Reclassify Failed")
error.write(arcpy.GetMessages())
str1 = precRaster+" "+ shapeList[i] + " undo"
eundo.write(str1)
continue
print "判断个数"
# 获取分类后的栅格独立值个数，判断是否有80这个值，有的话输出shape,没有则不输出
arcpy.BuildRasterAttributeTable_management(outfile, "Overwrite")
# 获取唯一值函数
uniq = arcpy.GetRasterProperties_management(outfile, "UNIQUEVALUECOUNT")
# 栅格转矢量
print "转栅格"
if str(uniq) != str(0):
outPolygons = arcpy.env.workspace + "\\ShapeFile\\" + precRaster.split(".")[0] + ".shp"
arcpy.RasterToPolygon_conversion(outReclass, outPolygons, "NO_SIMPLIFY", "Value")
# else:
# print precRaster + "NO 80"
#合并所有shapefiles
else:
print precRaster + "no 80"
str2 = precRaster + " "
enoarti.write(str2)
arcpy.env.workspace = arcpy.env.workspace + "\\ShapeFile"
featureList = arcpy.ListFeatureClasses("*")
arcpy.Merge_management(featureList, "Global30all.shp")
print "完成啦!"
eundo.close()
error.close()
enoarti.close()
else:
print "个数不一致"
