新浪博客

中转博客

2023-05-29 11:20阅读：

http://blog.sina.cn/dpool/blog/u/3115700593

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_Fl
oat32, 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; // 返回程序执行结果
}

我的更多文章

下载客户端阅读体验更佳

APP专享