Python气象数据处理与绘图：基于DEM数据计算坡向、坡度和流向

气海无涯 2021-09-11

2824

1、导入模块

1import matplotlib.pyplot as plt
2import richdem as rd

注：需要提前安装好richdem库，使用教程：

1https://richdem.com/

2、读取DEM数据

1dem = rd.LoadGDAL("/home/mw/work/data/NASADEM.nc")
2dem

输出：

 1rdarray([[28.5633104 , 25.2849467 , 21.96061724, ..., 23.35898004,
 2          18.71621789, 21.994     ],
 3         [28.092     , 21.82340848, 19.94186751, ..., 28.57500373,
 4          23.25103111, 21.09174748],
 5         [21.17122643, 16.98091248, 18.10103525, ..., 20.95329856,
 6          20.77973854, 17.99523309],
 7         ...,
 8         [ 0.        ,  0.        ,  0.        , ..., 15.06485127,
 9          19.94248513, 21.77396547],
10         [ 0.        ,  0.        ,  0.        , ..., 16.40997739,
11          16.82554527, 15.00107439],
12         [ 0.        ,  0.        ,  0.        , ..., 20.        ,
13          15.80715168, 16.0361648 ]])

3、绘制DEM数据

1fig = plt.figure(figsize=(14,10))
2plt.imshow(dem, interpolation='none')
3cbar = plt.colorbar()
4cbar.ax.tick_params(labelsize=15) 
5cbar.set_label(label='DEM (m)', size=15)
6plt.tick_params(labelsize=15)
7plt.show()

4、计算坡向（aspect）

1aspect = rd.TerrainAttribute(dem, attrib='aspect')
2rd.rdShow(aspect, axes=False, cmap='jet', figsize=(12,8))

5、计算剖面曲率（profile_curvature）

1profile_curvature = rd.TerrainAttribute(dem, attrib='profile_curvature')
2rd.rdShow(profile_curvature, axes=False, cmap='jet', figsize=(12,8))

6、计算坡度（slope）

1slope = rd.TerrainAttribute(dem, attrib='slope_riserun')
2rd.rdShow(slope, axes=False, cmap='magma', figsize=(12,8))

7、计算流向（Flow）

(1)方法一：D8, O’Callaghan, J.F., Mark, D.M., 1984. The Extraction of Drainage Networks from Digital Elevation Data. Computer vision, graphics, and image processing 28, 323–344.

1accum_d8 = rd.FlowAccumulation(dem, method='D8')
2d8_fig = rd.rdShow(accum_d8, figsize=(12,8), axes=False, cmap='jet')

(2)方法二：Rho8, Fairfield, J., Leymarie, P., 1991. Drainage networks from grid digital elevation models. Water resources research 27, 709–717.

1accum_rho8 = rd.FlowAccumulation(dem, method='Rho8')
2rd.rdShow(accum_rho8, figsize=(12,8), axes=False, cmap='jet')

(3)方法三：Quinn, Quinn, P., Beven, K., Chevallier, P., Planchon, O., 1991. The Prediction Of Hillslope Flow Paths For Distributed Hydrological Modelling Using Digital Terrain Models. Hydrological Processes 5, 59–79.

1accum_quinn = rd.FlowAccumulation(dem, method='Quinn')
2rd.rdShow(accum_quinn, figsize=(12,8), axes=False, cmap='jet')

有问题可以到QQ群里进行讨论，我们在那边等大家。

QQ群号：854684131

python

文章转载自气海无涯，如果涉嫌侵权，请发送邮件至：contact@modb.pro进行举报，并提供相关证据，一经查实，墨天轮将立刻删除相关内容。