Python初遇见

去年花了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 "个数不一致"