前一段时间做提取坡度的问题，当时首先想到的是使用ArcEngine来做，因为记得有ITopoOperator接口可以构建缓冲带，用IExtractionRaster可以掩膜栅格数据，利用IPixelBlock3接口可以读取栅格信息，计算像元的平均值。当时花了一段时间实现了，有时间把AE的这段代码分享出来，但是效率不敢恭维。输入4个多边形数据，半个小时愣是只跑了3个结果出来，当时我果断放弃了这条路。此时wangye学长跟我说，python处理栅格数据效率很高，基于C为底层，别说长江中下游区域，就是全国也能很快实现。由此，花了两天学了一下，便动手做了。最后45个要素组成的要素类和坡度栅格数据运算仅仅花了11分钟，这差距！！
数据：湖泊数据，Slope。
要求：湖泊指定距离的缓冲带范围内提取平均坡度
输出：LakeName_BufferDistance_SlopeMean.xls
Arcgis10开始推出了Arcpy模块，利用python调用相应的方法写出脚本语言，可以快速的进行批处理，相对于Ae的环境配置麻烦而言，是一个不错的选择，学习起来很简单，这个和IDL非常类似。
代码：
#!usr/bin/env python
#-*- coding:utf-8 -*-
‘create a buffer’
__author__ = ‘zhigang’
import arcpy
from arcpy import env
from FilenameWithoutExtension import GetFilenameWithoutExtension
import time
#overwrite output
#
env.overwriteOutput = True
inputLayer = ‘D:\\zgcao\\Zhigang\\WaterSheld\\AllData\\database.gdb\\Lakes_50’
distance = [dis for dis in range(5,55,5)]
output = ‘D:\\output\\’
start_time=time.localtime(time.time())
nowTime = time.strftime(‘%Y-%m-%d %H:%M:%S’,(start_time))
print ‘Start:’+str(nowTime)
print ‘Processing….’
for item in distance:
buffer_distance = str(item)+’ Kilometers’
outputName = output+str(item)+’_KM.shp’
#create buffer
#
arcpy.Buffer_analysis(inputLayer, outputName,buffer_distance)
#difference
#
#Set Full filename
#
differenceName = GetFilenameWithoutExtension(outputName)+’_Diff.shp’
#Symdifference
#
arcpy.SymDiff_analysis(outputName,inputLayer,differenceName)
# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension(“Spatial”)
#ZonalStatics
#
slopeData = ‘D:/zgcao/Zhigang/WaterSheld/AllData/database.gdb/slope’
fieldName = ‘Name_CH’
Sta_Result = GetFilenameWithoutExtension(differenceName)+’_sta’
arcpy.sa.ZonalStatisticsAsTable(differenceName,fieldName,slopeData,Sta_Result,’NODATA’,’MEAN’)
print ‘END’
end=time.localtime(time.time())
print str(time.strftime(‘%Y-%m-%d %H:%M:%S’,(end)))
使用到的自定义函数：
def GetFilenameWithoutExtension(filenameFull):
filename = filenameFull.split(‘.’)
return filename[0]
