Python：构建缓冲带提取区域平均坡度

时间：2023-12-20 14:32:20

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

相关文章