使用 C++ 和 GDAL 库进行基于余弦校正的地形校正的示例代码

在C++中利用GDAL库进行基于余弦校正的地形校正是地理信息系统处理中常见的任务。以下是一个简单的示例，展示了如何读取DEM数据（数字高程模型），进行余弦校正，并保存到新的文件中。请注意，这个示例假设你已经安装了GDAL库并且对基本的C++和GIS编程有所了解。

#include <gdal/gdal.h>
#include <gdal_priv.h>

// 函数声明，用于计算余弦校正因子
double CosineCorrection(double latitude, double demResolution);

int main() {
    // 初始化GDAL
    GDALAllRegister();

    // 打开原始DEM数据
    char* srcDSName = "path_to_your_dsm.tif";
    GDALDataset* srcDS = (GDALDataset*) GDALOpen(srcDSName, GA_ReadOnly);
    if (!srcDS) {
        CPLError(CE_Failure, CPLE_OpenFailed, "Unable to open DSM dataset");
        return -1;
    }

    // 获取数据集信息
    int xSize = srcDS->RasterXSize;
    int ySize = srcDS->RasterYSize;
    double ulX = srcDS->GetGeoTransform()[0];
    double urX = ulX + srcDS->GetGeoTransform()[1] * xSize;
    double lrX = ulX + srcDS->GetGeoTransform()[1] * (xSize - 1);
    double llX = ulX + srcDS->GetGeoTransform()[1] * (ySize - 1);
    double ulY = srcDS->GetGeoTransform()[3] + srcDS->GetGeoTransform()[5] * ySize;
    double urY = ulY;
    double lrY = ulY;
    double llY = ulY + srcDS->GetGeoTransform()[5] * (ySize - 1);

    // 创建目标栅格
    char* dstDSName = "path_to_save_corrected_data.tif";
    GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
    GDALDataset* dstDS = driver->CreateCopy(dstDSName, srcDS, 0);

    // 遍历所有像素并应用余弦校正
    for (int i = 0; i < ySize; ++i) {
        for (int j = 0; j < xSize; ++j) {
            // 纬度转换为弧度
            double lat = (ulY - srcDS->GetY(j)) / srcDS->GetGeoTransform()[5];

            // 读取原点值作为起始高度
            double rawHeight = srcDS->GetRasterBand(1)->ReadAsDouble(i, j);

            // 计算余弦校正因子
            double corrFactor = CosineCorrection(lat, srcDS->GetRasterBand(1)->GetXSize() / (srcDS->GetGeoTransform()[1] * 360));

            // 应用校正
            double correctedHeight = rawHeight / corrFactor;

            // 写入目标文件
            dstDS->GetRasterBand(1)->WriteRaster(j, i, 1, 1, &correctedHeight, 1, GDT_Float32, NULL);
        }
    }

    // 清理资源
    delete srcDS;
    delete dstDS;
    return 0;
}

// 辅助函数，计算余弦校正因子
double CosineCorrection(double latitude, double resolution) {
    // 地球半径（公里）
    const double earthRadiusKm = 6371;
    // 海拔影响修正，一般可以忽略
    const double altitudeCorrection = 0;

    // 将纬度转换为弧度
    double latRad = D2R(latitude);

    // 计算余弦校正因子
    return cos(latRad) * (earthRadiusKm + altitudeCorrection) / resolution;
}