编写利用DEM计算坡度的python代码

下面是一段利用Python计算DEM坡度的示例代码，其中使用了GDAL库来读取和处理DEM数据：

import gdal
import numpy as np

# 读取DEM数据
dem_file = "dem.tif"
dem_dataset = gdal.Open(dem_file)
dem_band = dem_dataset.GetRasterBand(1)
dem_data = dem_band.ReadAsArray().astype(np.float32)

# 计算坡度
x_scale, y_scale = dem_dataset.GetGeoTransform()[1], -dem_dataset.GetGeoTransform()[5]
dz_dx, dz_dy = np.gradient(dem_data, x_scale, y_scale)
slope = np.arctan(np.sqrt(dz_dx ** 2 + dz_dy ** 2)) * 180 / np.pi

# 输出坡度数据
slope_file = "slope.tif"
driver = gdal.GetDriverByName("GTiff")
slope_dataset = driver.Create(slope_file, dem_dataset.RasterXSize, dem_dataset.RasterYSize, 1, gdal.GDT_Float32)
slope_dataset.SetGeoTransform(dem_dataset.GetGeoTransform())
slope_dataset.SetProjection(dem_dataset.GetProjection())
slope_band = slope_dataset.GetRasterBand(1)
slope_band.WriteArray(slope)
slope_band.FlushCache()

# 释放资源
del dem_data, dem_band, dem_dataset
del slope_band, slope_dataset

该代码首先读取DEM数据，然后使用Numpy库计算DEM数据梯度，最后将计算出的坡度数据写入到GeoTIFF文件中。需要注意的是，该代码只计算了坡度的大小，而没有区分坡度的方向。如果需要计算坡度方向，则需要进一步处理梯度数据。