利用Python ArcPy计算地形起伏度最佳统计单元Word下载.docx
- 文档编号:6516538
- 上传时间:2023-05-06
- 格式:DOCX
- 页数:11
- 大小:22.46KB
利用Python ArcPy计算地形起伏度最佳统计单元Word下载.docx
《利用Python ArcPy计算地形起伏度最佳统计单元Word下载.docx》由会员分享,可在线阅读,更多相关《利用Python ArcPy计算地形起伏度最佳统计单元Word下载.docx(11页珍藏版)》请在冰点文库上搜索。
250万地貌图的地势起伏度是指16km2中的最高点与最低点的高差(Demek等,1989)[4]。
涂汉明等在编制1:
400万中国及其邻区地貌图时,定义地形起伏度(ReliefAmplitude)是指在某一确定面积内所有栅格中最大高程与最小高程之差,它是定量描述地貌形态、划分地貌类型的重要指标,在宏观区域内反映地面的起伏特征,并确定全国地形起伏最佳统计单元为21km2[5]。
21世纪以后,绝大多数研究都是采用一定范围内最大高程与最小高程之差来计算地形起伏度。
这个范围即统计单元,存在尺度问题,统计单元过小会切割完整的地貌单元,统计单元过大,会导致多个地貌单元出现而引起起伏度过大。
因此,起伏度统计单元的确定对于地貌分类及制图是非常重要的参数。
随着GIS技术的发展及大量高精度DEM数据的产生,许多学者进行了广泛的研究(见表1)。
根据研究可以发现,数据源以国家1:
25万地形图生产的DEM最为广泛,其次是SRTM(ShuttleRadarTopographyMission)DEM和ASTERGDEM数据。
SRTMDEM数据于2000年2月获取,全球数据产品包括1弧秒(30m)和3弧秒(90m)两种精度,数据覆盖了北纬60°
~南纬56°
之间的所有区域,数据绝对高程精度为16m,绝对平面精度为20m,2003年开始公布,陆续发布了多个数据版本,目前已经发布的最新版本为SRTM4.1。
ASTERGDEM12009年分别由NASA和METI分别发布,数据范围,北纬83°
~南纬83°
,数据空间分辨率为1弧秒(30m),GDEM2版本于2011年发布,绝对高程精度为20m,绝对平面精度为度30m[6]。
目前,确定起伏度单元大小的方法有4种:
人工作图法,最大高差法,模糊数学法以及均值变点法[7]。
分析范围从地区到全国乃至欧洲,分析单元大小也不尽相同。
但有规律可循,首先,数据比例尺越小,分析单元越大;
比例尺越大,分析单元越小。
按照比例尺,1:
5万尺度(25-30m),单元集中于2-3km2,1:
25万尺度(90-100m),单元集中在4-6km2;
从1:
100万地貌制图角度,4-6km2统计单元太小,1:
100万制图选择的数据源为1:
50万,DEM分辨率为200m,则相应的最小统计单元应该在20km2左右。
由表1可知,由于数据的精度,数据范围等的不同,统计单元会有所变化。
因此,在进行具体地区地貌制图时,首先要对使用数据进行分析从而确定最佳统计单元,而不能采用以往的经验。
本文根据选取均值变点法来确定最佳统计单元,并利用Python语言及ArcPy、numpy、Scipy、matplotlib完成。
1开发及运行环境配置
首先要安装ArcGIS10.1,其自带Python2.7,然后再下载其他几个安装包。
ArcGIS的Python环境是32位,可从网站上下载win32版本的几个软件包,并进行安装。
具体如下:
matplotlib-1.5.0.win32-py2.7.exe
numpy-1.9.2-win32-superpack-python2.7.exe
scikit-learn-0.18.1.win32-py2.7.exe
scipy-0.15.0-win32-superpack-python2.7.exe
2计算原理及方法
选择均值变点法确定移动窗口大小。
其过程如下:
1)首先统计n×
n(n=2,3,4,…,39,40)窗口内像元的最大值max、最小值min;
然后计算出max与min的差值,该值即为n×
n窗口的地形起伏度值;
2)从步骤1得到的每个起伏度栅格grid读取平均起伏度值,并计算单位起伏度=平均起伏度/移动窗口面积(本文采用圆形移动窗口);
3)将步骤2的结果构建一个序列;
式中:
Ti代表分析窗口下的单位起伏度;
Ti代表分析窗口下的平均起伏度;
Si代表分析窗口下的网格单元面积;
i为网格单元的大小。
4)对序列T取对数In(T),得序列X,X为{Xi,i=1,2,3,…,39},利用Python绘制单元面积与平均起伏度对应关系拟合曲线。
5)在上述处理的基础上,利用式
(2)计算样本序列X的算术平均值X,利用式(3)和式(4)分别计算样本序列X的统计量S与Si的值,Si值计算过程如下:
首先,令i=2,3,…,N,对每个i将样本Xi分为两段:
x1,x2,…,xi-1和xi,xi+1,…,xN。
然后计算每段样本的算术平均值及统计量,分别计算每个序列的序列值xt与算数平均值的离差的平方和,并相加得到Si。
6)计算S-Si,得到总样本与两样本的离差平方和的差值。
该序列值表示样本S与Si样本的偏差。
绘制S_Si与移动窗口半径R关系图,对应的最大值所对应的移动窗口单元大小为最佳统计单元。
3程序代码
1)采用圆形移动窗口,半径依次为2,3,…,39,41,求各个窗口下地形起伏度GRID(grd),平均起伏度(w_mean)。
采用循环程序进行计算,其具体代码如下:
foriinrange(2,41):
grd=dir1+str(i)
neighborhood=NbrCircle(i,"
CELL"
)
outFocalStatistics=FocalStatistics(inRaster,neighborhood,"
RANGE"
"
"
outFocalStatistics.save(grd)
#平均起伏度w_mean
w_mean=arcpy.GetRasterProperties_management(grd,"
MEAN"
2)求序列Ti的对数lnTi,Ti=ti/si,其中ti为平均起伏度(w_mean),si=移动窗口面积(areakm)。
代码如下:
lnTi=math.log(float(w_mean)/areakm)
3)前两步得到的数据输出为文件(f_mean_R.txt),根据文件进行对数拟合,并计算拟合度R2,绘制图形。
#用对数形式来拟合
f_mean_R=file('
D:
\\demclass\\class\\testd\\f_mean_R.txt'
x=y=[]
foriinf_mean_R.readlines():
j1=float(i.split()[1])
x.append(j1)
j2=float(i.split()[2])
y.append(j2)
f_mean_R.close()
x=np.array(x)
y=np.array(y)
deffunc(x,a,b):
returna*log(x)+b
popt,pcov=curve_fit(func,x,y)
a=popt[0]#popt里面是拟合系数
b=popt[1]
yvals=func(x,a,b)
defR2_fun(y,y_forecast):
#拟合优度R^2
y_mean=np.mean(y)
return1-(np.sum((y_forecast-y)**2))/(np.sum((y-y_mean)**2))
y_forecast=func(x,a,b)
R2=R2_fun(y,y_forecast)
strs1='
y='
+str(round(a,2))+'
*ln(x)+'
+str(round(b,2))
strs2='
R^2='
+str(round(R2,4))
plot1=plt.plot(x,y,'
+'
plot2=plt.plot(x,yvals,'
r-'
plt.text(30,10,strs1)
plt.text(30,-30,strs2)
plt.xlabel(u'
移动窗口半径'
fontproperties='
SimHei'
fontsize=14)
plt.ylabel(u'
平均地形起伏度(m)'
Sim-Hei'
#plt.xlabel('
movingwindowR'
plt.legend(loc=4)#指定legend的位置,读者可以自己help它的用法
plt.title(u'
地形起伏度-移动窗口'
plt.savefig('
scatter_curve.png'
4)将序列Xt分为两个序列Xt_1-i,Xt_i-n,分别计算序列的离差平方和S、Si1、Si2。
forjinrange(0,N+1):
X_mean=0;
S=0
foriinrange(0,j):
X_mean=Xt[i]/j+X_mean
S=(Xt[i]-X_mean)*(Xt[i]-X_mean)+S
Si1.append(S)
S_total=S
foriinrange(j,N):
X_mean=Xt[i]/(N-j)+X_mean
Si2.append(S)
5)求两段样本的离差平方和之和Si1+Si2,以及总样本与两段样本离差平方和之和的差值,即:
S-(Si1+Si2),并根据最大差值,确定最佳统计单元。
S_S1_2=[]
x=[]
y=[]
foriinrange(0,N):
Si1_Si2.append([Si1[i],Si2[i]])
#计算离差平方和之和Si=Si1+Si2
S1_plus_S2.append(Si1[i]+Si2[i])
#计算S-Si
S_S1_2.append(S_total-S1_plus_S2[i])
#求序列S-Si的最大值的移动窗口大小
window_size=S_S1_2.index(max(S_S1_2))+2
print'
最佳移动窗口半径R='
window_size
6)第五步结果输出的文件,绘制S-(Si1+Si2)与移动窗口半径散点图。
f_S_Si=file('
\\demclass\\class\\testd\\f_S_Si.txt'
foriinf_S_Si.readlines():
x1=float(i.split('
'
)[0])
x.append(x1)
y1=float(i.split('
)[1])
y.append(y1)
f_S_Si.close()
#drawpicture
x_c=y.index(max(y))+2
y_c=max(y)
xy=[x_c,y_c]
plot1=plt.plot(x,y,'
spot='
R='
+str(R)
plt.annotate(spot,xy,xycoords='
data'
xytext=[20,20],textcoords='
offsetpoints'
arrowprops=dict(arrowstyle='
->
'
))
plt.ylabel('
S-Si'
S_Si.png'
4实例分析
采用SRTM90mDEM数据,数据范围为整个青藏高原,如图1所示。
根据程序首先得到了统计单元与地形起伏度对应关系,见表2。
根据程序,绘制单元面积与平均起伏度对应关系拟合曲线图(见图2),以及确定最佳统计单元。
S-Si图,见图3。
计算地形起伏度最佳统计单元需要对栅格数据进行迭代处理,并完成统计。
ArcPy软件包是基于ArcGIS的Python软件包,对于ESRI格式的空间数据处理就变得异常简单,直接调用已有的FocalStatistics函数即可。
但由于窗口计算从2*2到40*40,因此要迭代39次,如果依靠人工计算,每个过程要3-10min,总耗时4h30min。
对计算结果进行单元面积与平均起伏度关系拟合,手动计算在Excel下进行,则大概需要10min左右才能完成制图和拟合。
进行统计样本的离差(S)计算最繁琐,总样本需要计算一次,而如果总样本分成两个动态样本,进行离差计算,则共有39*2=78个动态样本,每个样本计算一次离差S,S-Si,手动计算在Excel下进行需要花费10min左右,则需耗时78*10=780min。
最终绘图需要10min左右。
因此,进行青藏高原DEM地形起伏度最佳窗口统计计算,人工计算需要耗时大概270+10+780+10=956min,约18h。
编写程序代码3-4h,程序运行4h31min。
以上计算基于如下运行环境:
1)中央处理器CPUIntel(R)Core(TM)i7-4810MQCPU@2.8GHz2.8GHz;
2)内存:
16GB;
3)操作系统:
WIN764bit操作系统;
4)应用软件:
ArcGIS10.1,Python2.7.2,Office2007。
5结论
利用Python及其函数库进行开源GIS程序开发效率尤其明显[21-22],空间处理软件包由ArcPy提供,直接调用即可,统计函数和拟合函数由scipy、numpy软件包提供,绘图函数由matplotlib软件包提供,根据具体要求进行简单的参数调整,即可得到所需结果。
这样能大大减轻工作量,且能保证计算准确。
这些工作往往占了整个GIS数据管理系统的70%的工作量。
在ArcGIS数据处理与分析中正逐步取代AML语言的批处理功能,完成更复杂的工作。
[参考文献]
[1]李钜章.中国地貌形态基本类型数量指标初探[J].地理学报,1982,37
(1):
17-26.
[2]李钜章.中国地貌基本形态划分的探讨[J].地理科学,1987,6
(2):
32-39.
[3]中国科学院地理研究所.中国1:
1000000地貌图制图规范(试行)[M].北京:
科学出版社,1987:
1-44.
[4]DemekJ,EmbletonC.EmbletonC.InternationalGeomorphlogicalMapofEurope(l:
2500000)[M].Cartography,LithographyandPraba,S.P.,1989.
[5]涂汉明,刘振东.中国地势起伏度最佳统计单元的求证[J].湖北大学学报:
自然科学,1990,12(3):
266-271.
[6]CarlosH.Grohmann,EvaluationofTanDEM-XDEMsonselectedBraziliansites:
ComparisonwithSRTM,ASTERGDEMandALOSAW3D30[J].RemoteSensingofEnvironment,212(2018):
121-133.
[7]赵斌滨,程永锋,丁士君,等.基于SRTM-DEM的我国地势起伏度统计单元研究[J].水利学报,2015,46(sup1):
284-290.
[8]陈志明.论中国地貌图的研制原则、内容与方法——以1:
4000000全国地貌图为例[J].地理学报,1993,48
(2):
105-213.
[9]刘新华,杨勤科,汤国安.中国地形起伏度的提取及在水土流失定量评价中的应用[J].水土保持通报,2001,21
(1):
57-60.
[10]高守英,吴泉源,安国强.基于GIS的龙口市泳汶河流域地貌形态定量分析[J].遥感技术与应用,2003,18
(2):
87-90.
[11]唐飞,陈曦,程维明,等.基于DEM的准噶尔盆地及其西北山区地势起伏度研究[J].干旱区地理,2006,29(3):
388-391.
[12]龙恩,程维明,周成虎,等.江基于Srtm-DEM与遥感的长白山基本地貌类型提取方法[J].山地学报,2007,25(5):
557-565.
[13]王玲,同小娟.基于变点分析的地形起伏度研究[J].地理与地理信息科学,2007,23(6):
65-67.
[14]莫申国.基于DEM的秦岭数字地貌格局研究[J].华东师范大学学报:
自然科学版,2008
(2):
8-14.
[15]郎玲玲,程维明,朱启疆,等.多尺度dem提取地势起伏度的对比分析——以福建低山丘陵区为例[J].地球信息科学,2007(6):
135-136.
[16]王岩,刘少峰.基于DEM的青海贵德地区地形起伏度的研究[J].地质通报,2008,27(12):
2117-2121.
[17]程维明,周成虎,柴慧霞,等.中国陆地地貌基本形态类型定量提取与分析[J].地球信息科学学报,2009,11(6):
725-736.
[18]曹伟超,陶和平,孔博,等.青藏高原地貌形态总体特征的GIS识别分析[J].水土保持通报,2011,31(4):
163-168.
[19]张伟,李爱农.基于DEM的中国地形起伏度适宜计算尺度研究[J].地理与地理信息科学,2012,28(4):
9-14.
[20]王让虎,张树文,蒲罗曼,等.基于ASTERGDEM和均值变点分析的中国东北地形起伏度研究[J].干旱区资源与环境,2016,30(6):
49-54.
[21]Python真的比其他语言慢吗[EB/OL].https:
//
[22]Python学习:
给你一个学习Python的理由[EB/OL].http:
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 利用Python ArcPy计算地形起伏度最佳统计单元 利用 Python ArcPy 计算 地形 起伏 最佳 统计 单元