前一段时间做提取坡度的问题,当时首先想到的是使用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]