地形起伏度是指在一定区域范围内的最大高程与最小高程之差. 反映在DEM上, 就是指分析区域内, 栅格最大值与最小值的差异, 表示分析区域的高程起伏情况. 地形起伏度的计算公式为: RF = Hmax - Hmin .
地形起伏度受分析范围影响较大, 通常随着分析范围的扩大而扩大. 本文采用均值变点法, 此方法多用于地形起伏度的求算, 它对仅一个变化点的检验最为有效.
arcpy批量计算地形起伏度
计算2×2到51×51共50个网格单元下的平均地形起伏度, 将窗口大小与对应的地形起伏度均值保存在csv格式文件中.
import arcpy
import os
from arcpy import sa
import csv
# DEM路径
DEM = ur'E:\assignment\dem\test\dem_js\dem_js.tif'
# 输出文件夹路径
outputPath = ur'E:\assignment\dem\test\RFs'
# 计算指定格网的地形起伏度
def calc_gridRF(gridLength, DEM, outputFolder=None):
nbr = sa.NbrRectangle(gridLength, gridLength, 'CELL') # 邻域分析的窗口大小
rasterMax = sa.BlockStatistics(DEM, nbr, 'MAXIMUM')
rasterMin = sa.BlockStatistics(DEM, nbr, 'MINIMUM')
RF = rasterMax - rasterMin
if outputFolder:
output = os.path.join(outputPath, ur'RF_%s' % (gridLength))
RF.save(output)
return RF.mean
if __name__ == "__main__":
header = ['cellsize', 'Mean_RF']
rows = []
for i in range(2, 52):
RFMean = calc_gridRF(i, DEM)
row = (i, RFMean)
rows.append(row)
with open(ur"mean_rf.csv", 'w') as f:
writer = csv.writer(f)
writer.writerow(header)
writer.writerows(rows)
对数函数拟合
确定最佳分析窗口的大致范围.
import csv
import numpy as np
import matplotlib.pyplot as plt
if __name__ == "__main__":
with open("mean_rf.csv", 'r') as f:
reader = csv.reader(f)
x = []
y = []
for i in reader:
x.append(i[0])
y.append(i[1])
x = np.array(list(map(lambda a: float(a), x[1:])))
y = np.array(list(map(lambda a: float(a), y[1:])))
log_x = np.log(x)
coefficients = np.polyfit(log_x, y, 1)
if coefficients[1] >= 0:
print("y = %f * ln(x) + %f" % (float(coefficients[0]), abs(float(coefficients[1]))))
else:
print("y = %f * ln(x) - %f" % (float(coefficients[0]), abs(float(coefficients[1]))))
c = coefficients[0] * log_x + coefficients[1]
r2 = 1 - np.sum((y - c)**2) / np.sum((y - np.mean(y))**2)
print("R2:", r2)
plt.figure(figsize=(8, 4))
plt.xlim(0, 55)
plt.ylim(-5, 30)
plt.xlabel('CellSize')
plt.ylabel('Mean of RF/m')
plt.plot(x, y, "o")
plt.plot(x, c, linewidth=2)
plt.savefig("窗口大小与地形起伏度均值.png")
plt.show()
均值变点法
根据均值变点法的定义,![img](file:///C:/Users/Khru/AppData/Local/Temp/msohtmlclip1/01/clip_image002.png)处于最大值时所对应的窗口大小即为最佳分析窗口。而S是固定值,因此只需要算出![img](file:///C:/Users/Khru/AppData/Local/Temp/msohtmlclip1/01/clip_image004.png)的最小值,此时的窗口大小就是最佳分析窗口。
import numpy
import math
import rfs # 代码1
import csv
class calcReliefAmplitude():
def __init__(self, numrange, gridLength):
self.numrange = numrange
self.len = len(self.numrange)
self.gridLength = gridLength
self.bestWindow = self.calcBestWindow()
def _calcVar(self):
return numpy.var(self.numrange) * self.len
def calcBestWindow(self):
S = []
for i in range(2, len(self.numrange) + 1):
s1 = numpy.var(self.numrange[:i]) * (i - 1)
s2 = numpy.var(self.numrange[i:]) * (len(self.numrange) - i + 1)
S.append(s1 + s2)
bestWindow = S.index(min(S)) + 2
return bestWindow
if __name__ == "__main__":
with open("mean_rf.csv", 'r') as f:
reader = csv.reader(f)
rasterList = []
gridLength = []
for i in reader:
rasterList.append(i[1])
gridLength.append(i[0])
rasterList = list(map(lambda a: float(a), rasterList[1:]))
gridLength = list(map(lambda a: float(a), gridLength[1:]))
DEMbestWindow = calcReliefAmplitude(rasterList, gridLength).bestWindow
print(DEMbestWindow)
rfs.calc_gridRF(DEMbestWindow, rfs.DEM, rfs.outputPath)
标签:__,plt,窗口,Python,self,gridLength,import,提取,csv
From: https://www.cnblogs.com/khrushchefox/p/16990445.html