高程DEM文件操作

1 下载高程DEM文件网站

• https://www.gscloud.cn/

• 注册邮箱，即可登录下载

• 选择对应的经纬度文件下载即可
• 直接高级检索，下载对应区域
  

2 查找安徽池州高程数据

3 高程文件合并与读取

• 高程文件合并，gdal_merge 安装gdal库时自带的脚步程序

xhome@ubuntu:~/dem$ gdal_merge.py -o merged_dem.tif ./dem_chizhou/ASTGTMV003_*_dem.tif
0...10...20...30...40...50...60...70...80...90...100 - done.

#include <iostream>
#include <gdal/gdal_priv.h>
#include <gdal/ogr_spatialref.h>
#include <iomanip>  // 用于控制输出格式


// g++ gdal_demo.cpp -o gdal_demo -lgdal
int main(int argc, char* argv[]) {
    // 检查命令行参数
    if (argc != 2) {
        std::cerr << "用法: " << argv[0] << " <输入TIFF文件路径>" << std::endl;
        return 1;
    }

    const char* tifFilePath = argv[1];

    // 注册GDAL驱动
    GDALAllRegister();

    // 打开TIFF文件
    GDALDataset* dataset = (GDALDataset*)GDALOpen(tifFilePath, GA_ReadOnly);
    if (dataset == nullptr) {
        std::cerr << "错误：无法打开文件 " << tifFilePath << std::endl;
        return 1;
    }

    // 获取地理变换参数
    double geoTransform[6];
    if (dataset->GetGeoTransform(geoTransform) != CE_None) {
        std::cerr << "错误：无法获取地理变换参数" << std::endl;
        GDALClose(dataset);
        return 1;
    }

    // 获取图像尺寸
    int width = dataset->GetRasterXSize();
    int height = dataset->GetRasterYSize();

    // 获取高程数据波段（假设是第一个波段）
    GDALRasterBand* band = dataset->GetRasterBand(1);
    if (band == nullptr) {
        std::cerr << "错误：无法获取高程数据波段" << std::endl;
        GDALClose(dataset);
        return 1;
    }

    // 获取坐标参考系统
    const char* proj = dataset->GetProjectionRef();
    OGRSpatialReference srcSRS;
    if (proj != nullptr && strlen(proj) > 0) {
        if (srcSRS.importFromWkt(&proj) != OGRERR_NONE) {
            std::cerr << "警告：无法解析坐标参考系统，默认使用WGS84" << std::endl;
            srcSRS.SetWellKnownGeogCS("WGS84");
        }
    } else {
        std::cerr << "警告：未找到坐标参考系统，默认使用WGS84" << std::endl;
        srcSRS.SetWellKnownGeogCS("WGS84");
    }

    // 准备目标坐标系（WGS84）
    OGRSpatialReference dstSRS;
    dstSRS.SetWellKnownGeogCS("WGS84");

    // 创建坐标转换对象
    OGRCoordinateTransformation* coordTrans = 
        OGRCreateCoordinateTransformation(&srcSRS, &dstSRS);
    if (coordTrans == nullptr) {
        std::cerr << "警告：无法创建坐标转换对象，将输出原始坐标" << std::endl;
    }

    // 分配行缓冲区
    float* pafScanline = new float[width];
    if (pafScanline == nullptr) {
        std::cerr << "错误：无法分配内存" << std::endl;
        if (coordTrans) OCTDestroyCoordinateTransformation(coordTrans);
        GDALClose(dataset);
        return 1;
    }

    // 遍历所有像素
    for (int y = 0; y < height; y++) {
        // 读取一行高程数据
        if (band->RasterIO(GF_Read, 0, y, width, 1, pafScanline, width, 1, GDT_Float32, 0, 0) != CE_None) {
            std::cerr << "警告：第" << y << "行读取失败" << std::endl;
            continue;
        }

        for (int x = 0; x < width; x++) {
            // 计算像素坐标对应的地理坐标
            double pixelX = geoTransform[0] + x * geoTransform[1] + y * geoTransform[2];
            double pixelY = geoTransform[3] + x * geoTransform[4] + y * geoTransform[5];

            // 如果需要坐标转换
            if (coordTrans) {
                // 创建临时变量用于转换
                double tmpX = pixelX;
                double tmpY = pixelY;
                if (!coordTrans->Transform(1, &tmpX, &tmpY)) {
                    std::cerr << "警告：坐标转换失败，使用原始坐标" << std::endl;
                } else {
                    pixelX = tmpX;
                    pixelY = tmpY;
                }
            }

            // 获取高程值
            float elevation = pafScanline[x];

            // 输出结果（经度、纬度、高程）
            std::cout << std::fixed << std::setprecision(6) 
                      << pixelX << "," << pixelY << "," 
                      << std::setprecision(2) << elevation << std::endl;
        }
    }

    // 释放资源
    delete[] pafScanline;
    if (coordTrans) OCTDestroyCoordinateTransformation(coordTrans);
    GDALClose(dataset);

    return 0;
}