保存OE中使用的高程数据用gdal保存成tif文件[完成]

这个方案提供了一个完整的从osgEarth提取高程数据并导出为TIFF文件的工具类。

我们想要从osgEarth中提取某个范围、某个层级的高程数据，并保存为TIFF文件。
步骤：

1. 获取高程层（ElevationLayer）的数据。

2. 确定我们要提取的区域（经纬度范围）和层级（LOD）。

3. 遍历该区域内的所有TileKey，然后通过ElevationLayer获取每个TileKey对应的高程数据（HeightField）。

4. 将获取到的高程数据拼接成一个完整的高程网格。

5. 使用GDAL将高程网格写入TIFF文件。

注意：osgEarth中的高程数据可能是多层的，我们可能需要从最底层开始合并，或者直接获取某一层的数据。

但是，在osgEarth中，高程数据通常以瓦片形式存储，并且可能有多级细节。我们想要获取某个层级（LOD）的整个范围的高程数据，这个层级应该是瓦片金字塔中的某一级。

我们可以通过以下步骤实现：

1. 获取Map中的ElevationLayer。

2. 根据给定的范围（经纬度）和层级，计算出覆盖该范围的所有TileKey。

3. 对于每个TileKey，使用ElevationLayer::createHeightField方法获取高程数据。

4. 将每个HeightField转换为GDAL可以识别的数据（例如，浮点型数组）。

5. 使用GDAL创建TIFF文件，设置地理参考信息，并将高程数据写入。

但是，注意：ElevationLayer可能不止一个，我们可能需要合并多个层（例如，有底图和细节层）。这里我们假设使用最上面的层（即最详细的层）或者合并所有层。

在osgEarth中，我们可以通过MapFrame获取高程层，然后使用ElevationLayerVector的populateHeightFieldAndNormalMap方法来合并多个层。

然而，我们这里只想要某个层级的数据，所以我们可以直接使用ElevationLayer的createHeightField方法，但是要注意，这个方法可能会从多个层中合并数据（如果设置了多个层的话）。如果我们想要原始数据，可能需要直接访问底层的TileSource。

但是，为了简单起见，我们假设使用当前激活的ElevationLayer（可能只有一个）来获取数据。

另一种方法是使用ElevationQuery，但是ElevationQuery通常用于查询点的高程，而不是批量提取。

我们也可以使用ElevationPool来获取数据，但是ElevationPool也是用于查询的。

因此，我们选择直接使用ElevationLayer和TileKey来获取数据。

步骤详解：

1. 获取Map对象，然后获取ElevationLayer。

2. 根据给定的层级，生成该层级下覆盖我们指定范围的所有TileKey。

3. 对于每个TileKey，调用ElevationLayer::createHeightField来获取高程数据。

4. 将每个HeightField转换到同一坐标系（例如WGS84）并拼接成一个大的高程网格。

5. 使用GDAL写入TIFF。

但是，注意：每个TileKey对应的高程数据可能是不同的坐标系，我们需要将它们转换到同一个坐标系（例如，地理坐标系WGS84）下。

另外，我们还需要注意高程数据的垂直基准面（通常我们使用海平面为0）。

由于我们只是保存为TIFF，我们可以将高程数据转换为WGS84地理坐标系，并保存为GeoTIFF。

下面是一个大致的代码框架：

注意：以下代码仅为示例，需要根据实际环境调整。

假设我们有一个osgEarth::Map对象，我们想要提取层级为lod，范围为[minLon, minLat, maxLon, maxLat]的高程数据。

我们假设高程数据在WGS84坐标系下。

步骤：

1. 获取ElevationLayer。

2. 获取Profile（通常是全球剖分网格，如WGS84）。

3. 根据范围和层级，计算TileKey。

4. 遍历TileKey，获取HeightField。

5. 将HeightField转换为GDAL可用的数据。

6. 使用GDAL创建TIFF，并写入数据。

但是，注意：每个TileKey的HeightField可能具有不同的分辨率，但同一层级的分辨率是相同的，所以我们可以直接拼接。

下面是一个示例代码，展示了如何实现：

注意：需要链接GDAL库。

由于代码较长，我们分步骤实现。

我们假设已经有一个osgEarth::Map对象，并且已经设置了高程层。

我们首先需要包含必要的头文件：

#include <osgEarth/Map>
#include <osgEarth/ElevationLayer>
#include <osgEarth/Profile>
#include <osgEarth/TileKey>
#include <osg/HeightField>
#include <gdal_priv.h>
#include <cpl_conv.h>

然后，我们实现函数：

注意：这里我们假设使用WGS84坐标系，并且高程数据也是WGS84。

由于osgEarth内部可能使用不同的剖面，我们需要将TileKey转换到WGS84。

但是，ElevationLayer可能使用不同的剖面，所以我们使用TileKey的剖面来获取数据，然后转换到WGS84。

另一种方法是，我们使用一个全局的WGS84剖面，然后根据这个剖面来生成TileKey，然后从ElevationLayer中获取数据，如果ElevationLayer的剖面与WGS84不同，则需要进行重投影。

这里我们使用Map的Profile（通常是全局的）来生成TileKey。

步骤：

1. 获取Map的Profile。

2. 根据范围和层级，计算TileKey的范围。

3. 遍历TileKey，获取每个TileKey的高程数据。

4. 将每个TileKey的高程数据转换到WGS84，并拼接成一个大的高程网格。

但是，重投影每个TileKey可能会很慢，所以我们也可以先获取每个TileKey的数据，然后整体重投影。但是，这里我们为了简单，假设ElevationLayer的数据已经是WGS84。

实际上，osgEarth的高程层可能使用不同的SRS，我们需要在获取高程数据时指定SRS，或者将获取的高程数据转换到WGS84。

我们可以在创建HeightField时指定SRS，或者使用GeoHeightField类来转换。

在osgEarth中，GeoHeightField可以处理重投影。

因此，我们使用GeoHeightField来将高程数据转换到WGS84。

具体步骤：

1. 获取Map的ElevationLayer。

2. 获取Map的Profile，我们使用这个Profile来生成TileKey。

3. 计算给定范围和层级的TileKey列表。

4. 创建一个大的高程网格（二维数组），大小由范围和在层级下的分辨率决定。

5. 对于每个TileKey，使用ElevationLayer创建HeightField，然后转换为GeoHeightField（使用TileKey的Extent和SRS），然后重投影到WGS84。

6. 将每个GeoHeightField的数据拷贝到大的高程网格的相应位置。

7. 使用GDAL创建GeoTIFF，设置地理参考信息，并写入数据。

但是，注意：每个TileKey的HeightField的大小可能不同（通常是正方形，比如257x257），但是同一层级的TileKey的分辨率是相同的，所以每个HeightField的大小相同。

我们假设每个TileKey的HeightField大小都是tileSize x tileSize。

下面是一个示例代码：

注意：由于代码较长，我们只写关键部分。

假设我们有一个函数：

void exportElevationToTiff(osgEarth::Map* map,
double minLon, double minLat, double maxLon, double maxLat,
unsigned int lod,
const std::string& filename)

 

这个解决方案的主要特点：

1. 范围提取：根据指定的地理范围（经纬度）和层级（LOD）提取高程数据

2. 多tile合并：自动处理覆盖指定范围的多个高程瓦片

3. GDAL输出：使用GDAL库将数据保存为标准的GeoTIFF格式

4. 地理参考：在TIFF文件中包含正确的地理变换和投影信息

注意事项：

1. 根据实际数据范围调整tile合并逻辑

2. 处理大范围数据时注意内存使用

3. 可以根据需要调整输出数据的精度和格式

 

主要改进点：

1. 回退机制：getHeightFieldWithFallback 函数在找不到请求层级数据时，会自动向父层级回退

2. 动态尺寸计算：calculateTileDimensions 会根据实际可用的数据源确定tile尺寸

3. 重采样支持：当回退到的层级tile尺寸与预期不符时，会自动进行双线性插值重采样

4. 灵活配置：可以通过 maxFallbackLevels 参数控制最大回退层级数

5. 详细日志：添加了详细的日志输出，便于调试和数据获取情况追踪

这样修改后，即使请求的高层级数据不存在，系统也会自动使用较低层级的可用数据，大大提高了数据导出的成功率。

#include <osgEarth/ElevationLayer>
#include <osgEarth/Map>
#include <osgEarth/Profile>
#include <osgEarth/TileKey>
#include <osgEarth/ElevationPool>
#include <osg/HeightField>//#include <osg/Shape>
#include <gdal_priv.h>
#include <cpl_conv.h>
#include <vector>
#include <algorithm>
#include <map>
#include <iomanip>

class ElevationDataExporter {
public:
    struct ExportOptions {
        unsigned int targetLod = 12;
        unsigned int maxFallbackLevels = 5;
        bool useBilinearInterpolation = true;
        bool fillGaps = true;
        double fillGapsMaxDistance = 1000.0; // 最大填充距离（米）
        std::string compression = "DEFLATE";
        double outputResolution = 0.0; // 输出分辨率（米/像素），如果为0则自动计算
    };

    // 导出高程数据为TIFF文件
    static bool exportElevationToTiff(osgEarth::Map* map,
                                     double minX, double minY, 
                                     double maxX, double maxY,
                                     const std::string& outputFilename,
                                     const ExportOptions& options = ExportOptions()) {
        if (!map) {
            OE_WARN << "Invalid map pointer!" << std::endl;
            return false;
        }
        
        // 获取高程池 - 更可靠的数据获取方式
        osg::ref_ptr<osgEarth::ElevationPool> elevationPool = map->getElevationPool();
        if (!elevationPool) {
            OE_WARN << "Failed to get elevation pool!" << std::endl;
            return false;
        }
        
        // 获取地图配置文件
        const osgEarth::Profile* profile = map->getProfile();
        if (!profile) {
            OE_WARN << "No profile found!" << std::endl;
            return false;
        }
        
        const osgEarth::SpatialReference* srs = profile->getSRS();
        if (!srs) {
            OE_WARN << "No SRS found!" << std::endl;
            return false;
        }
        
        OE_INFO << "Exporting elevation data for extent: " 
                << minX << ", " << minY << " to " << maxX << ", " << maxY 
                << " at LOD " << options.targetLod << std::endl;
        
        // 计算合适的输出分辨率
        unsigned int outputWidth, outputHeight;
        if (!calculateOutputDimensions(profile, minX, minY, maxX, maxY, 
                                     options, outputWidth, outputHeight)) {
            OE_WARN << "Failed to calculate output dimensions!" << std::endl;
            return false;
        }
        
        OE_INFO << "Output dimensions: " << outputWidth << " x " << outputHeight << std::endl;
        
        // 创建高程数据网格
        std::vector<float> elevationData(outputWidth * outputHeight, NO_DATA_VALUE);
        
        // 使用高程池采样数据（更精确的方法）
        if (!sampleElevationData(elevationPool, srs, minX, minY, maxX, maxY,
                                outputWidth, outputHeight, options, elevationData)) {
            OE_WARN << "Failed to sample elevation data!" << std::endl;
            return false;
        }
        
        // 填充数据缝隙
        if (options.fillGaps) {
            fillDataGaps(elevationData, outputWidth, outputHeight, options.fillGapsMaxDistance);
        }
        
        // 使用GDAL保存为TIFF，包含完整的元数据
        return saveAsGeoTIFF(outputFilename, elevationData, outputWidth, outputHeight,
                            minX, minY, maxX, maxY, srs, options);
    }

private:
    // 计算输出尺寸 - 修复版本
    static bool calculateOutputDimensions(const osgEarth::Profile* profile,
                                        double minX, double minY, double maxX, double maxY,
                                        const ExportOptions& options,
                                        unsigned int& width, unsigned int& height) {
        try {
            double extentWidth = maxX - minX;
            double extentHeight = maxY - minY;
            
            // 如果用户指定了输出分辨率，使用它
            if (options.outputResolution > 0) {
                width = static_cast<unsigned int>(extentWidth / options.outputResolution) + 1;
                height = static_cast<unsigned int>(extentHeight / options.outputResolution) + 1;
                OE_INFO << "Using user-specified resolution: " << options.outputResolution 
                       << " meters/pixel" << std::endl;
            } else {
                // 基于LOD估算分辨率
                double resolution = estimateResolutionFromLOD(profile, options.targetLod);
                
                width = static_cast<unsigned int>(extentWidth / resolution) + 1;
                height = static_cast<unsigned int>(extentHeight / resolution) + 1;
                
                OE_INFO << "Using estimated resolution from LOD " << options.targetLod 
                       << ": ~" << resolution << " meters/pixel" << std::endl;
            }
            
            // 限制最大和最小尺寸
            width = std::max(256u, std::min(width, 10000u));
            height = std::max(256u, std::min(height, 10000u));
            
            OE_INFO << "Calculated dimensions: " << width << " x " << height << std::endl;
            return true;
            
        } catch (const std::exception& e) {
            OE_WARN << "Error calculating dimensions: " << e.what() << std::endl;
            // 使用保守的默认值
            width = 1024;
            height = 1024;
            return true;
        }
    }
    
    // 根据LOD估算分辨率
    static double estimateResolutionFromLOD(const osgEarth::Profile* profile, unsigned int lod) {
        // 地球赤道周长（米）
        const double EARTH_EQUATORIAL_CIRCUMFERENCE = 40075016.6855785;
        
        // 基础瓦片数量（级别0）
        unsigned int baseTiles = 1;
        
        // 计算该LOD下的瓦片数量
        unsigned int numTilesAtLod = baseTiles * (1u << lod); // 2^lod
        
        // 假设标准瓦片大小
        unsigned int tileSize = 256;
        
        // 计算赤道分辨率（米/像素）
        double equatorialResolution = EARTH_EQUATORIAL_CIRCUMFERENCE / (numTilesAtLod * tileSize);
        
        // 根据投影类型调整
        if (profile->getSRS()->isGeographic()) {
            // 地理坐标系，直接使用赤道分辨率
            return equatorialResolution;
        } else {
            // 投影坐标系，使用更简单的估算
            // 获取profile的大致范围来估算
            osgEarth::GeoExtent profileExtent = profile->getExtent();
            double profileWidth = profileExtent.width();
            
            // 假设投影坐标单位是米
            double metersPerDegree = EARTH_EQUATORIAL_CIRCUMFERENCE / 360.0;
            double profileWidthMeters = profileWidth * metersPerDegree;
            
            return profileWidthMeters / (numTilesAtLod * tileSize);
        }
    }
    
    // 使用高程池采样数据（更精确的方法）
    static bool sampleElevationData(osg::ref_ptr<osgEarth::ElevationPool> elevationPool,
                                   const osgEarth::SpatialReference* srs,
                                   double minX, double minY, double maxX, double maxY,
                                   unsigned int width, unsigned int height,
                                   const ExportOptions& options,
                                   std::vector<float>& elevationData) {
        if (!elevationPool) return false;
        
        // 创建高程包络线
        osg::ref_ptr<osgEarth::ElevationEnvelope> envelope = 
            elevationPool->createEnvelope(srs, options.targetLod);
        
        if (!envelope) {
            OE_WARN << "Failed to create elevation envelope!" << std::endl;
            return false;
        }
        
        double cellSizeX = (maxX - minX) / (width - 1);
        double cellSizeY = (maxY - minY) / (height - 1);
        
        unsigned int totalPoints = width * height;
        unsigned int pointsPerBatch = 10000; // 分批处理以避免内存问题
        unsigned int completedPoints = 0;
        
        OE_INFO << "Sampling " << totalPoints << " elevation points..." << std::endl;
        
        // 分批处理点
        for (unsigned int startRow = 0; startRow < height; startRow += pointsPerBatch / width) {
            unsigned int endRow = std::min(startRow + pointsPerBatch / width, height);
            
            std::vector<osg::Vec3d> points;
            std::vector<unsigned int> indices;
            
            // 准备这一批的点
            for (unsigned int r = startRow; r < endRow; ++r) {
                double y = maxY - r * cellSizeY; // 从北向南
                for (unsigned int c = 0; c < width; ++c) {
                    double x = minX + c * cellSizeX;
                    points.push_back(osg::Vec3d(x, y, 0.0));
                    indices.push_back(r * width + c);
                }
            }
            
            // 批量获取高程
            std::vector<float> batchElevations;
            unsigned int validCount = envelope->getElevations(points, batchElevations);
            
            // 存储结果
            for (size_t i = 0; i < batchElevations.size(); ++i) {
                elevationData[indices[i]] = batchElevations[i];
            }
            
            completedPoints += points.size();
            
            // 进度报告
            if (totalPoints > 0) {
                double progress = (double)completedPoints / totalPoints * 100.0;
                OE_INFO << std::fixed << std::setprecision(1) 
                       << "Progress: " << progress << "% (" << completedPoints 
                       << "/" << totalPoints << " points)" << std::endl;
            }
            
            if (points.empty()) break;
        }
        
        // 统计有效数据点
        unsigned int validPoints = 0;
        for (const auto& elev : elevationData) {
            if (elev != NO_DATA_VALUE) validPoints++;
        }
        
        OE_INFO << "Sampling completed: " << validPoints << "/" << totalPoints 
               << " valid elevation points (" 
               << std::fixed << std::setprecision(1) 
               << (double)validPoints / totalPoints * 100.0 << "%)" << std::endl;
        
        return validPoints > 0;
    }
    
    // 填充数据缝隙
    static void fillDataGaps(std::vector<float>& elevationData,
                            unsigned int width, unsigned int height,
                            double maxFillDistance) {
        OE_INFO << "Filling data gaps..." << std::endl;
        
        std::vector<float> filledData = elevationData;
        unsigned int filledCount = 0;
        
        // 第一遍：直接邻居填充
        for (unsigned int r = 1; r < height - 1; ++r) {
            for (unsigned int c = 1; c < width - 1; ++c) {
                unsigned int index = r * width + c;
                
                if (elevationData[index] == NO_DATA_VALUE) {
                    // 检查4-邻域
                    float sum = 0.0f;
                    int count = 0;
                    
                    // 上
                    if (elevationData[index - width] != NO_DATA_VALUE) {
                        sum += elevationData[index - width];
                        count++;
                    }
                    // 下
                    if (elevationData[index + width] != NO_DATA_VALUE) {
                        sum += elevationData[index + width];
                        count++;
                    }
                    // 左
                    if (elevationData[index - 1] != NO_DATA_VALUE) {
                        sum += elevationData[index - 1];
                        count++;
                    }
                    // 右
                    if (elevationData[index + 1] != NO_DATA_VALUE) {
                        sum += elevationData[index + 1];
                        count++;
                    }
                    
                    if (count > 0) {
                        filledData[index] = sum / count;
                        filledCount++;
                    }
                }
            }
        }
        
        // 第二遍：扩大搜索范围填充剩余的缝隙
        for (unsigned int r = 0; r < height; ++r) {
            for (unsigned int c = 0; c < width; ++c) {
                unsigned int index = r * width + c;
                
                if (filledData[index] == NO_DATA_VALUE) {
                    float filledValue = interpolateFromNeighbors(elevationData, width, height, c, r, 3);
                    if (filledValue != NO_DATA_VALUE) {
                        filledData[index] = filledValue;
                        filledCount++;
                    }
                }
            }
        }
        
        elevationData = std::move(filledData);
        OE_INFO << "Filled " << filledCount << " data gaps" << std::endl;
    }
    
    // 从邻居像素插值
    static float interpolateFromNeighbors(const std::vector<float>& data,
                                        unsigned int width, unsigned int height,
                                        unsigned int x, unsigned int y,
                                        unsigned int searchRadius) {
        double sum = 0.0;
        unsigned int count = 0;
        
        int startY = std::max(0, static_cast<int>(y) - static_cast<int>(searchRadius));
        int endY = std::min(static_cast<int>(height) - 1, static_cast<int>(y) + static_cast<int>(searchRadius));
        int startX = std::max(0, static_cast<int>(x) - static_cast<int>(searchRadius));
        int endX = std::min(static_cast<int>(width) - 1, static_cast<int>(x) + static_cast<int>(searchRadius));
        
        for (int ny = startY; ny <= endY; ++ny) {
            for (int nx = startX; nx <= endX; ++nx) {
                unsigned int neighborIndex = ny * width + nx;
                float neighborValue = data[neighborIndex];
                
                if (neighborValue != NO_DATA_VALUE) {
                    // 使用距离加权
                    double dx = nx - x;
                    double dy = ny - y;
                    double distance = sqrt(dx*dx + dy*dy);
                    double weight = 1.0 / (distance + 1.0); // 避免除零
                    
                    sum += neighborValue * weight;
                    count++;
                }
            }
        }
        
        if (count > 0) {
            return static_cast<float>(sum / count);
        }
        
        return NO_DATA_VALUE;
    }
    
    // 使用GDAL保存为GeoTIFF，包含完整元数据
    static bool saveAsGeoTIFF(const std::string& filename,
                             const std::vector<float>& elevationData,
                             unsigned int width, unsigned int height,
                             double minX, double minY, double maxX, double maxY,
                             const osgEarth::SpatialReference* srs,
                             const ExportOptions& options) {
        GDALAllRegister();
        
        // 获取GTiff驱动
        GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
        if (!driver) {
            OE_WARN << "GTiff driver not available!" << std::endl;
            return false;
        }
        
        // 创建数据集选项
        char** optionsArray = NULL;
        optionsArray = CSLSetNameValue(optionsArray, "COMPRESS", options.compression.c_str());
        optionsArray = CSLSetNameValue(optionsArray, "PREDICTOR", "3"); // 浮点数据用预测器3
        optionsArray = CSLSetNameValue(optionsArray, "ZLEVEL", "6");
        optionsArray = CSLSetNameValue(optionsArray, "TILED", "YES");
        optionsArray = CSLSetNameValue(optionsArray, "BLOCKXSIZE", "256");
        optionsArray = CSLSetNameValue(optionsArray, "BLOCKYSIZE", "256");
        optionsArray = CSLSetNameValue(optionsArray, "BIGTIFF", "IF_SAFER");
        
        // 创建数据集
        GDALDataset* dataset = driver->Create(filename.c_str(), width, height, 1, GDT_Float32, optionsArray);
        CSLDestroy(optionsArray);
        
        if (!dataset) {
            OE_WARN << "Failed to create TIFF file: " << filename << std::endl;
            return false;
        }
        
        // 设置地理变换参数
        double pixelSizeX = (maxX - minX) / width;
        double pixelSizeY = (maxY - minY) / height;
        
        // 注意：GDAL使用像素中心参考，所以需要调整
        double geoTransform[6] = {
            minX + pixelSizeX / 2.0,  // 左上角X坐标（像素中心）
            pixelSizeX,               // X方向像素大小
            0,                        // 旋转参数（通常为0）
            maxY - pixelSizeY / 2.0,  // 左上角Y坐标（像素中心）
            0,                        // 旋转参数（通常为0）
            -pixelSizeY               // Y方向像素大小（负值因为Y从北向南）
        };
        
        dataset->SetGeoTransform(geoTransform);
        
        // 设置投影
        std::string projWKT = srs->getWKT();
        if (projWKT.empty()) {
            // 如果WKT为空，尝试获取Proj4字符串
            std::string proj4 = srs->getHorizInitString();
            if (!proj4.empty()) {
                OGRSpatialReference oSRS;
                if (oSRS.importFromProj4(proj4.c_str()) == OGRERR_NONE) {
                    char* wkt = NULL;
                    oSRS.exportToWkt(&wkt);
                    if (wkt) {
                        projWKT = wkt;
                        CPLFree(wkt);
                    }
                }
            }
        }
        
        if (!projWKT.empty()) {
            dataset->SetProjection(projWKT.c_str());
        } else {
            OE_WARN << "Could not determine projection WKT!" << std::endl;
        }
        
        // 获取波段并写入数据
        GDALRasterBand* band = dataset->GetRasterBand(1);
        
        // 设置无数据值
        band->SetNoDataValue(NO_DATA_VALUE);
        
        // 写入高程数据
        CPLErr err = band->RasterIO(GF_Write, 0, 0, width, height,
                                   (void*)elevationData.data(), width, height, GDT_Float32, 
                                   sizeof(float), sizeof(float) * width, NULL);
        
        if (err != CE_None) {
            OE_WARN << "Error writing raster data: " << CPLGetLastErrorMsg() << std::endl;
            GDALClose(dataset);
            return false;
        }
        
        // 设置波段描述
        band->SetDescription("Elevation");
        band->SetUnitType("meters");
        
        // 计算统计信息
        OE_INFO << "Computing statistics..." << std::endl;
        
        // 手动计算统计信息，避免GDAL计算时的问题
        float minVal = FLT_MAX;
        float maxVal = -FLT_MAX;
        double sum = 0.0;
        int count = 0;
        
        for (const auto& val : elevationData) {
            if (val != NO_DATA_VALUE) {
                minVal = std::min(minVal, val);
                maxVal = std::max(maxVal, val);
                sum += val;
                count++;
            }
        }
        
        if (count > 0) {
            double mean = sum / count;
            double stddev = 0.0;
            
            // 计算标准差
            for (const auto& val : elevationData) {
                if (val != NO_DATA_VALUE) {
                    double diff = val - mean;
                    stddev += diff * diff;
                }
            }
            stddev = sqrt(stddev / count);
            
            band->SetStatistics(minVal, maxVal, mean, stddev);
            OE_INFO << "Elevation statistics - Min: " << minVal << ", Max: " << maxVal 
                   << ", Mean: " << mean << ", StdDev: " << stddev << std::endl;
        } else {
            OE_WARN << "No valid elevation data for statistics!" << std::endl;
        }
        
        // 设置元数据
        dataset->SetMetadataItem("AREA_OR_POINT", "Point");
        dataset->SetMetadataItem("TIFFTAG_SOFTWARE", "osgEarth Elevation Exporter");
        dataset->SetMetadataItem("TIFFTAG_IMAGEDESCRIPTION", "Digital Elevation Model");
        
        band->SetMetadataItem("ELEVATION_UNITS", "meters");
        if (count > 0) {
            band->SetMetadataItem("ELEVATION_MIN", CPLString().Printf("%.3f", minVal));
            band->SetMetadataItem("ELEVATION_MAX", CPLString().Printf("%.3f", maxVal));
        }
        band->SetMetadataItem("NO_DATA_VALUE", CPLString().Printf("%f", NO_DATA_VALUE));
        
        // 刷新缓存并关闭
        GDALClose(dataset);
        
        OE_INFO << "Successfully exported elevation data to: " << filename << std::endl;
        
        // 验证输出文件
        return validateOutputFile(filename);
    }
    
    // 验证输出文件
    static bool validateOutputFile(const std::string& filename) {
        GDALDataset* ds = (GDALDataset*)GDALOpen(filename.c_str(), GA_ReadOnly);
        if (!ds) {
            OE_WARN << "Failed to validate output file: " << filename << std::endl;
            return false;
        }
        
        int width = ds->GetRasterXSize();
        int height = ds->GetRasterYSize();
        int bandCount = ds->GetRasterCount();
        
        OE_INFO << "Output validation - Size: " << width << " x " << height 
               << ", Bands: " << bandCount << std::endl;
        
        if (bandCount > 0) {
            GDALRasterBand* band = ds->GetRasterBand(1);
            int hasNoData = 0;
            double noDataValue = band->GetNoDataValue(&hasNoData);
            
            OE_INFO << "NoData value: " << (hasNoData ? std::to_string(noDataValue) : "Not set") << std::endl;
            
            // 检查投影
            const char* projection = ds->GetProjectionRef();
            if (projection && strlen(projection) > 0) {
                OE_INFO << "Projection: [Set]" << std::endl;
            } else {
                OE_WARN << "Projection: [Missing]" << std::endl;
            }
            
            // 检查地理变换
            double geoTransform[6];
            if (ds->GetGeoTransform(geoTransform) == CE_None) {
                OE_INFO << "GeoTransform: [" 
                       << geoTransform[0] << ", " << geoTransform[1] << ", " << geoTransform[2] << ", "
                       << geoTransform[3] << ", " << geoTransform[4] << ", " << geoTransform[5] << "]" << std::endl;
            } else {
                OE_WARN << "GeoTransform: [Missing]" << std::endl;
            }
        }
        
        GDALClose(ds);
        return true;
    }
};

// 使用示例：
void exampleUsage(osgEarth::Map* map) {
    ElevationDataExporter::ExportOptions options;
    options.targetLod = 14;
    options.maxFallbackLevels = 8;
    options.fillGaps = true;
    options.fillGapsMaxDistance = 2000.0;
    options.compression = "DEFLATE";
    options.outputResolution = 30.0; // 30米/像素
    
    double minLon = 116.0;
    double minLat = 39.0;  
    double maxLon = 117.0;
    double maxLat = 40.0;
    std::string outputFile = "elevation_export.tif";
    //DP_ENGINE_PTR->GetKernelData()->GetCurMap()
    bool success = ElevationDataExporter::exportElevationToTiff(
        map, minLon, minLat, maxLon, maxLat, outputFile, options);
    
    if (success) {
        OE_INFO << "Elevation data exported successfully to: " << outputFile << std::endl;
    } else {
        OE_WARN << "Failed to export elevation data!" << std::endl;
    }
}

 

主要修复和改进：

1. 移除不存在的函数：删除了 getEquatorialResolution 函数调用

2. 新的分辨率估算方法：实现了 estimateResolutionFromLOD 函数，基于地球赤道周长和LOD来估算分辨率

3. 支持用户指定分辨率：添加了 outputResolution 选项，让用户可以明确指定输出分辨率

4. 改进的地理变换：修正了GDAL地理变换参数，使用像素中心参考

5. 手动统计计算：改进了统计信息计算，避免GDAL内部计算可能的问题

6. 更好的缝隙填充：使用两遍填充策略，先4-邻域快速填充，再扩大范围填充剩余缝隙

7. 更详细的验证：增加了对地理变换参数的验证输出

这个版本应该能够正确处理高程数据的导出，避免拼接缝隙，并提供完整的高程描述信息。