使用python将高程（XYZ）数据转换为坡度/坡度图

Question

我有一个包含东距（x）、北距（y）和高程数据（z）的文本文件，如下所示：

   x            y         z
241736.69   3841916.11  132.05
241736.69   3841877.89  138.76
241736.69   3841839.67  142.89
241736.69   3841801.45  148.24
241736.69   3841763.23  157.92
241736.69   3841725.02  165.01
241736.69   3841686.80  171.86
241736.69   3841648.58  178.80
241736.69   3841610.36  185.26
241736.69   3841572.14  189.06
241736.69   3841533.92  191.28
241736.69   3841495.71  193.27
241736.69   3841457.49  193.15
241736.69   3841419.27  194.85
241736.69   3841381.05  192.31
241736.69   3841342.83  188.73
241736.69   3841304.61  183.68
241736.69   3841266.39  176.97
241736.69   3841228.18  160.83
241736.69   3841189.96  145.69
241736.69   3841151.74  129.09
241736.69   3841113.52  120.03
241736.69   3841075.30  111.84
241736.69   3841037.08  104.82
241736.69   3840998.86  101.63
241736.69   3840960.65  97.66
241736.69   3840922.43  93.38
241736.69   3840884.21  88.84
...

我可以用plt.contour和plt.contourf很容易地从上面的数据中得到一张海拔图，如下所示：

但是，我正在尝试获取我所掌握数据的斜率图，如下所示：

我试着用GDAL将XYZ数据转换成DEM，如解释的here，并用richdem加载DEM，如解释的here，但我得到的斜率值不对。在

转换为.tif得到的结果：

这是我用richdem尝试过的代码：

^{pr2}$

我得到的情节是：

colorbar上的值太高而不正确，必须反转绘图以匹配上述绘图（不是我现在的主要问题）。在

当我将python用于GIS时，我不是一个专家（我主要使用python进行数据分析），我希望这不会像我想象的那么复杂。在

Answer 1

eI能够正确地编写一个函数，但首先我需要赞扬这个answer，因为它节省了我编写自己的移动窗口函数的时间（非常好！）公司名称：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from tqdm import trange


def window3x3(arr, shape=(3, 3)):
    r_win = np.floor(shape[0] / 2).astype(int)
    c_win = np.floor(shape[1] / 2).astype(int)
    x, y = arr.shape
    for i in range(x):
        xmin = max(0, i - r_win)
        xmax = min(x, i + r_win + 1)
        for j in range(y):
            ymin = max(0, j - c_win)
            ymax = min(y, j + c_win + 1)
            yield arr[xmin:xmax, ymin:ymax]


def gradient(XYZ_file, min=0, max=15, figsize=(15, 10), **kwargs):

    """

    :param XYZ_file: XYZ file in the following format: x,y,z (inlcuding headers)
    :param min: color bar minimum range.
    :param max: color bar maximum range.
    :param figsize: figure size.
    :param kwargs:
           plot: to plot a gradient map. Default is True.
    :return: returns an array with the shape of the grid with the computed slopes


    The algorithm calculates the gradient using a first-order forward or backward difference on the corner points, first
    order central differences at the boarder points, and a 3x3 moving window for every cell with 8 surrounding cells (in
    the middle of the grid) using a third-order finite difference weighted by reciprocal of squared distance

    Assumed 3x3 window:

                                    -
                        |   a   |   b   |   c   |
                                    -
                        |   d   |   e   |   f   |
                                    -
                        |   g   |   h   |   i   |
                                    -


    """

    kwargs.setdefault('plot', True)

    grid = XYZ_file.to_numpy()

    nx = XYZ_file['x'].unique().size
    ny = XYZ_file['y'].unique().size

    xs = grid[:, 0].reshape(ny, nx, order='F')
    ys = grid[:, 1].reshape(ny, nx, order='F')
    zs = grid[:, 2].reshape(ny, nx, order='F')
    dx = abs((xs[:, 1:] - xs[:, :-1]).mean())
    dy = abs((ys[1:, :] - ys[:-1, :]).mean())

    gen = window3x3(zs)
    windows_3x3 = np.asarray(list(gen))
    windows_3x3 = windows_3x3.reshape(ny, nx)

    dzdx = np.empty((ny, nx))
    dzdy = np.empty((ny, nx))
    loc_string = np.empty((ny, nx), dtype="S25")

    for ax_y in trange(ny):
        for ax_x in range(nx):

            # corner points
            if ax_x == 0 and ax_y == 0:  # top left corner
                dzdx[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][0][1] - windows_3x3[ax_y, ax_x][0][0]) / dx
                dzdy[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][1][0] - windows_3x3[ax_y, ax_x][0][0]) / dy
                loc_string[ax_y, ax_x] = 'top left corner'

            elif ax_x == nx - 1 and ax_y == 0:  # top right corner
                dzdx[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][0][1] - windows_3x3[ax_y, ax_x][0][0]) / dx
                dzdy[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][1][1] - windows_3x3[ax_y, ax_x][0][1]) / dy
                loc_string[ax_y, ax_x] = 'top right corner'

            elif ax_x == 0 and ax_y == ny - 1:  # bottom left corner
                dzdx[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][1][1] - windows_3x3[ax_y, ax_x][1][0]) / dx
                dzdy[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][1][0] - windows_3x3[ax_y, ax_x][0][0]) / dy
                loc_string[ax_y, ax_x] = 'bottom left corner'

            elif ax_x == nx - 1 and ax_y == ny - 1:  # bottom right corner
                dzdx[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][1][1] - windows_3x3[ax_y, ax_x][1][0]) / dx
                dzdy[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][1][1] - windows_3x3[ax_y, ax_x][0][1]) / dy
                loc_string[ax_y, ax_x] = 'bottom right corner'

            # top boarder
            elif (ax_y == 0) and (ax_x != 0 and ax_x != nx - 1):
                dzdx[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][0][-1] - windows_3x3[ax_y, ax_x][0][0]) / (2 * dx)
                dzdy[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][1][1] - windows_3x3[ax_y, ax_x][0][1]) / dy
                loc_string[ax_y, ax_x] = 'top boarder'

            # bottom boarder
            elif ax_y == ny - 1 and (ax_x != 0 and ax_x != nx - 1):
                dzdx[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][1][-1] - windows_3x3[ax_y, ax_x][1][0]) / (2 * dx)
                dzdy[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][1][1] - windows_3x3[ax_y, ax_x][0][1]) / dy
                loc_string[ax_y, ax_x] = 'bottom boarder'

            # left boarder
            elif ax_x == 0 and (ax_y != 0 and ax_y != ny - 1):
                dzdx[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][1][1] - windows_3x3[ax_y, ax_x][1][0]) / dx
                dzdy[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][-1][0] - windows_3x3[ax_y, ax_x][0][0]) / (2 * dy)
                loc_string[ax_y, ax_x] = 'left boarder'

            # right boarder
            elif ax_x == nx - 1 and (ax_y != 0 and ax_y != ny - 1):
                dzdx[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][1][1] - windows_3x3[ax_y, ax_x][1][0]) / dx
                dzdy[ax_y, ax_x] = (windows_3x3[ax_y, ax_x][-1][-1] - windows_3x3[ax_y, ax_x][0][-1]) / (2 * dy)
                loc_string[ax_y, ax_x] = 'right boarder'

            # middle grid
            else:
                a = windows_3x3[ax_y, ax_x][0][0]
                b = windows_3x3[ax_y, ax_x][0][1]
                c = windows_3x3[ax_y, ax_x][0][-1]
                d = windows_3x3[ax_y, ax_x][1][0]
                f = windows_3x3[ax_y, ax_x][1][-1]
                g = windows_3x3[ax_y, ax_x][-1][0]
                h = windows_3x3[ax_y, ax_x][-1][1]
                i = windows_3x3[ax_y, ax_x][-1][-1]

                dzdx[ax_y, ax_x] = ((c + 2 * f + i) - (a + 2 * d + g)) / (8 * dx)
                dzdy[ax_y, ax_x] = ((g + 2 * h + i) - (a + 2 * b + c)) / (8 * dy)
                loc_string[ax_y, ax_x] = 'middle grid'

    hpot = np.hypot(abs(dzdy), abs(dzdx))
    slopes_angle = np.degrees(np.arctan(hpot))
    if kwargs['plot']:
        slopes_angle[(slopes_angle < min) | (slopes_angle > max)]

        plt.figure(figsize=figsize)
        plt.pcolormesh(xs, ys, slopes_angle, cmap=plt.cm.gist_yarg, vmax=max, vmin=min)

        plt.colorbar()
        plt.tight_layout()
        plt.show()

    return slopes_angle


if __name__ == '__main__':

    XYZ = pd.read_csv('xyz_file')
    slopes = gradient(XYZ)

最后的情节：

Answer 2

假设数据在一个nx3numpy数组中，首先将elevation列重新解释为矩阵（表示统一网格）：

m=data[:,2].reshape(ny,nx)

然后进行几次切片和减法，得到细胞中心的导数：

^{pr2}$

该系数校正单位（否则为米/单元格，而不是每米），并将总和转换为平均值。（如果每个维度中的间距不同，您将分别将参数缩放为hypot。）请注意，结果数组在每个维度上都比输入小一个；如果大小需要相同，可以使用更复杂的差分方案。numpy.gradient实现了其中的一些，允许

mag=numpy.hypot(*numpy.gradient(m,abs(data[1][1]-data[1][0])))