XYZ数据到python中的斜率梯度图

Question

我有一个带有Easting（x），Northing（y）和Elevation数据（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尝试过的代码：

import richdem as rd

dem_path = 'convertedXYZ.tif'
dem = rd.LoadGDAL(dem_path, no_data=-9999)
slope = rd.TerrainAttribute(dem, attrib='slope_riserun')
rd.rdShow(slope, axes=True, cmap='gist_yarg', figsize=(16, 9))

我得到的情节：

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

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

[可能解决方案]

我到目前为止的代码。我需要找到一种方法来验证这一点，但它对我来说似乎是正确的（也许其他人可以给我他们的反馈）。

XYZ = pd.read_csv('data')
grid = XYZ.to_numpy()

nx = XYZ['x'].unique().size
ny = XYZ['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')

dxs = np.abs(np.gradient(xs)[1])
dys = np.abs(np.gradient(ys)[0])
dzs = np.abs(np.gradient(zs)[1])

hpot = np.hypot(dys, dxs)
slopes = dzs / hpot
slopes_angle = np.degrees(np.arctan2(dzs, hpot))

Answer 1

假设您的数据是n x 3 Numpy数组，首先将高程列重新解释为矩阵（表示统一网格）：

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

然后执行几个切片和减法以获得细胞中心的衍生物：

dx=m[:,1:]-m[:,:-1]
dy=m[1:,:]-m[:-1,:]
mag=numpy.hypot(dx[1:,:]+dx[:-1,:],
                dy[:,1:]+dy[:,:-1])
mag*=abs(data[1][1]-data[1][0])/2

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

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