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

04-24 16:04