include "gdal_priv.h" // GDAL库头文件

include "cpl_conv.h" // for CPLMalloc()

int main() { // Register GDAL drivers 注册GDAL驱动 GDALAllRegister();

复制代码

// Open DEM file 打开DEM文件

GDALDataset *dem_ds = (GDALDataset*) GDALOpen("dem.tif", GA_ReadOnly); // GA_ReadOnly表示只读方式打开

if (dem_ds == NULL) { // 判断是否成功打开

    printf("Failed to open DEM file.

"); // 打印错误信息 exit(1); // 退出程序 }

复制代码

// Get DEM properties 获取DEM文件的属性

int dem_width = dem_ds->GetRasterXSize(); // DEM文件的宽度

int dem_height = dem_ds->GetRasterYSize(); // DEM文件的高度

double dem_transform[6]; // 存储DEM文件的仿射变换系数

dem_ds->GetGeoTransform(dem_transform); // 获取DEM文件的仿射变换系数

double dem_no_data = dem_ds->GetRasterBand(1)->GetNoDataValue(); // 获取DEM文件的无效值

// Create slope output file 创建输出坡度文件

GDALDriver *driver = GetGDALDriverManager()->GetDriverByName("GTiff"); // 获取输出文件驱动

GDALDataset *slope_ds = driver->Create("slope.tif", dem_width, dem_height, 1, GDT_Float32, NULL); // 创建输出文件

double slope_transform[6]; // 存储输出文件的仿射变换系数

memcpy(slope_transform, dem_transform, sizeof(double) * 6); // 将DEM文件的仿射变换系数赋值给输出文件的仿射变换系数

slope_ds->SetGeoTransform(slope_transform); // 设置输出文件的仿射变换系数

slope_ds->SetProjection(dem_ds->GetProjectionRef()); // 设置输出文件的投影信息

// Calculate slope 计算坡度

double dx, dy, dz; // 存储横向、纵向、高程的差值

float slope; // 存储坡度值

GDALRasterBand *dem_band = dem_ds->GetRasterBand(1); // 获取DEM文件的第1个波段

GDALRasterBand *slope_band = slope_ds->GetRasterBand(1); // 获取输出文件的第1个波段

int nXSize = dem_band->GetXSize();//获取宽度

int nYSize = dem_band->GetYSize();//获取高度

int nBlockSize = 1024;//每块的大小

int nXBlocks = (nXSize + nBlockSize - 1) / nBlockSize;//块的个数

int nYBlocks = (nYSize + nBlockSize - 1) / nBlockSize;

float *pafBuffer = new float[nBlockSize * nBlockSize];//缓存数组

for (int iYBlock = 0; iYBlock < nYBlocks; iYBlock++)

{

    for (int iXBlock = 0; iXBlock < nXBlocks; iXBlock++)

    {

        int nXValid = nBlockSize, nYValid = nBlockSize;

        if ((iXBlock + 1) * nBlockSize > nXSize)

            nXValid = nXSize - iXBlock * nBlockSize;

        if ((iYBlock + 1) * nBlockSize > nYSize)

            nYValid = nYSize - iYBlock * nBlockSize;

        dem_band->RasterIO(GF_Read, iXBlock * nBlockSize, iYBlock * nBlockSize, nXValid, nYValid, pafBuffer, nXValid, nYValid, GDT_Float32, 0, 0);

        slope_band->RasterIO(GF_Write, iXBlock * nBlockSize, iYBlock * nBlockSize, nXValid, nYValid, pafBuffer, nXValid, nYValid, GDT_Float32, 0, 0);

        for (int y = iYBlock * nBlockSize + 1; y < (iYBlock + 1) * nBlockSize - 1 && y < dem_height - 1; y++) { // 遍历DEM文件的每一行

            for (int x = iXBlock * nBlockSize + 1; x < (iXBlock + 1) * nBlockSize - 1 && x < dem_width - 1; x++) { // 遍历DEM文件的每一列

                // Calculate slope using central difference method 使用中心差分法计算坡度

                dx = pafBuffer[(y * nXValid + x + 1)] - pafBuffer[(y * nXValid + x - 1)]; // 计算横向差值

                dy = pafBuffer[((y + 1) * nXValid + x)] - pafBuffer[((y - 1) * nXValid + x)]; // 计算纵向差值

                dz = dem_transform[1] * dx + dem_transform[2] * dy; // 计算高程差值

                slope = atan(sqrt(dx*dx + dy*dy) / fabs(dz)) * 180.0 / M_PI; // 计算坡度值

                // Write slope value to slope file 将坡度值写入输出文件

                slope_band->SetPixel(x, y, slope);

            }

        }

    }

}

// Cleanup 清理内存

GDALClose(dem_ds); // 关闭DEM文件

GDALClose(slope_ds); // 关闭输出文件

GDALDestroyDriverManager(); // 销毁驱动管理器

return 0; // 返回程序执行结果

}