Python-提取地形起伏度最佳分析窗口

地形起伏度是指在一定区域范围内的最大高程与最小高程之差. 反映在DEM上, 就是指分析区域内, 栅格最大值与最小值的差异, 表示分析区域的高程起伏情况. 地形起伏度的计算公式为: RF = H_max - H_min .

地形起伏度受分析范围影响较大, 通常随着分析范围的扩大而扩大. 本文采用均值变点法, 此方法多用于地形起伏度的求算, 它对仅一个变化点的检验最为有效.

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()

均值变点法

根据均值变点法的定义，处于最大值时所对应的窗口大小即为最佳分析窗口。而S是固定值，因此只需要算出的最小值，此时的窗口大小就是最佳分析窗口。

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)