這個方案提供了一個完整的從osgEarth提取高程數據並導出為TIFF文件的工具類。
我們想要從osgEarth中提取某個範圍、某個層級的高程數據，並保存為TIFF文件。
步驟：
- 獲取高程層（ElevationLayer）的數據。
- 確定我們要提取的區域（經緯度範圍）和層級（LOD）。
- 遍歷該區域內的所有TileKey，然後通過ElevationLayer獲取每個TileKey對應的高程數據（HeightField）。
- 將獲取到的高程數據拼接成一個完整的高程網格。
- 使用GDAL將高程網格寫入TIFF文件。
注意：osgEarth中的高程數據可能是多層的，我們可能需要從最底層開始合併，或者直接獲取某一層的數據。
但是，在osgEarth中，高程數據通常以瓦片形式存儲，並且可能有多級細節。我們想要獲取某個層級（LOD）的整個範圍的高程數據，這個層級應該是瓦片金字塔中的某一級。
我們可以通過以下步驟實現：
- 獲取Map中的ElevationLayer。
- 根據給定的範圍（經緯度）和層級，計算出覆蓋該範圍的所有TileKey。
- 對於每個TileKey，使用ElevationLayer::createHeightField方法獲取高程數據。
- 將每個HeightField轉換為GDAL可以識別的數據（例如，浮點型數組）。
- 使用GDAL創建TIFF文件，設置地理參考信息，並將高程數據寫入。
但是，注意：ElevationLayer可能不止一個，我們可能需要合併多個層（例如，有底圖和細節層）。這裏我們假設使用最上面的層（即最詳細的層）或者合併所有層。
在osgEarth中，我們可以通過MapFrame獲取高程層，然後使用ElevationLayerVector的populateHeightFieldAndNormalMap方法來合併多個層。
然而，我們這裏只想要某個層級的數據，所以我們可以直接使用ElevationLayer的createHeightField方法，但是要注意，這個方法可能會從多個層中合併數據（如果設置了多個層的話）。如果我們想要原始數據，可能需要直接訪問底層的TileSource。
但是，為了簡單起見，我們假設使用當前激活的ElevationLayer（可能只有一個）來獲取數據。
另一種方法是使用ElevationQuery，但是ElevationQuery通常用於查詢點的高程，而不是批量提取。
我們也可以使用ElevationPool來獲取數據，但是ElevationPool也是用於查詢的。
因此，我們選擇直接使用ElevationLayer和TileKey來獲取數據。
步驟詳解：
- 獲取Map對象，然後獲取ElevationLayer。
- 根據給定的層級，生成該層級下覆蓋我們指定範圍的所有TileKey。
- 對於每個TileKey，調用ElevationLayer::createHeightField來獲取高程數據。
- 將每個HeightField轉換到同一座標系（例如WGS84）並拼接成一個大的高程網格。
- 使用GDAL寫入TIFF。
但是，注意：每個TileKey對應的高程數據可能是不同的座標系，我們需要將它們轉換到同一個座標系（例如，地理座標系WGS84）下。
另外，我們還需要注意高程數據的垂直基準面（通常我們使用海平面為0）。
由於我們只是保存為TIFF，我們可以將高程數據轉換為WGS84地理座標系，並保存為GeoTIFF。
下面是一個大致的代碼框架：
注意：以下代碼僅為示例，需要根據實際環境調整。
假設我們有一個osgEarth::Map對象，我們想要提取層級為lod，範圍為[minLon, minLat, maxLon, maxLat]的高程數據。
我們假設高程數據在WGS84座標系下。
步驟：
- 獲取ElevationLayer。
- 獲取Profile（通常是全球剖分網格，如WGS84）。
- 根據範圍和層級，計算TileKey。
- 遍歷TileKey，獲取HeightField。
- 將HeightField轉換為GDAL可用的數據。
- 使用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。
步驟：
- 獲取Map的Profile。
- 根據範圍和層級，計算TileKey的範圍。
- 遍歷TileKey，獲取每個TileKey的高程數據。
- 將每個TileKey的高程數據轉換到WGS84，並拼接成一個大的高程網格。
但是，重投影每個TileKey可能會很慢，所以我們也可以先獲取每個TileKey的數據，然後整體重投影。但是，這裏我們為了簡單，假設ElevationLayer的數據已經是WGS84。
實際上，osgEarth的高程層可能使用不同的SRS，我們需要在獲取高程數據時指定SRS，或者將獲取的高程數據轉換到WGS84。
我們可以在創建HeightField時指定SRS，或者使用GeoHeightField類來轉換。
在osgEarth中，GeoHeightField可以處理重投影。
因此，我們使用GeoHeightField來將高程數據轉換到WGS84。
具體步驟：
- 獲取Map的ElevationLayer。
- 獲取Map的Profile，我們使用這個Profile來生成TileKey。
- 計算給定範圍和層級的TileKey列表。
- 創建一個大的高程網格（二維數組），大小由範圍和在層級下的分辨率決定。
- 對於每個TileKey，使用ElevationLayer創建HeightField，然後轉換為GeoHeightField（使用TileKey的Extent和SRS），然後重投影到WGS84。
- 將每個GeoHeightField的數據拷貝到大的高程網格的相應位置。
- 使用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)
這個解決方案的主要特點：
- 範圍提取：根據指定的地理範圍（經緯度）和層級（LOD）提取高程數據
- 多tile合併：自動處理覆蓋指定範圍的多個高程瓦片
- GDAL輸出：使用GDAL庫將數據保存為標準的GeoTIFF格式
- 地理參考：在TIFF文件中包含正確的地理變換和投影信息
注意事項：
- 根據實際數據範圍調整tile合併邏輯
- 處理大範圍數據時注意內存使用
- 可以根據需要調整輸出數據的精度和格式
主要改進點：
- 回退機制：
getHeightFieldWithFallback函數在找不到請求層級數據時，會自動向父層級回退 - 動態尺寸計算：
calculateTileDimensions會根據實際可用的數據源確定tile尺寸 - 重採樣支持：當回退到的層級tile尺寸與預期不符時，會自動進行雙線性插值重採樣
- 靈活配置：可以通過
maxFallbackLevels參數控制最大回退層級數 - 詳細日誌：添加了詳細的日誌輸出，便於調試和數據獲取情況追蹤
這樣修改後，即使請求的高層級數據不存在，系統也會自動使用較低層級的可用數據，大大提高了數據導出的成功率。
1 #include <osgEarth/ElevationLayer>
2 #include <osgEarth/Map>
3 #include <osgEarth/Profile>
4 #include <osgEarth/TileKey>
5 #include <osgEarth/ElevationPool>
6 #include <osg/HeightField>//#include <osg/Shape>
7 #include <gdal_priv.h>
8 #include <cpl_conv.h>
9 #include <vector>
10 #include <algorithm>
11 #include <map>
12 #include <iomanip>
13
14 class ElevationDataExporter {
15 public:
16 struct ExportOptions {
17 unsigned int targetLod = 12;
18 unsigned int maxFallbackLevels = 5;
19 bool useBilinearInterpolation = true;
20 bool fillGaps = true;
21 double fillGapsMaxDistance = 1000.0; // 最大填充距離（米）
22 std::string compression = "DEFLATE";
23 double outputResolution = 0.0; // 輸出分辨率（米/像素），如果為0則自動計算
24 };
25
26 // 導出高程數據為TIFF文件
27 static bool exportElevationToTiff(osgEarth::Map* map,
28 double minX, double minY,
29 double maxX, double maxY,
30 const std::string& outputFilename,
31 const ExportOptions& options = ExportOptions()) {
32 if (!map) {
33 OE_WARN << "Invalid map pointer!" << std::endl;
34 return false;
35 }
36
37 // 獲取高程池 - 更可靠的數據獲取方式
38 osg::ref_ptr<osgEarth::ElevationPool> elevationPool = map->getElevationPool();
39 if (!elevationPool) {
40 OE_WARN << "Failed to get elevation pool!" << std::endl;
41 return false;
42 }
43
44 // 獲取地圖配置文件
45 const osgEarth::Profile* profile = map->getProfile();
46 if (!profile) {
47 OE_WARN << "No profile found!" << std::endl;
48 return false;
49 }
50
51 const osgEarth::SpatialReference* srs = profile->getSRS();
52 if (!srs) {
53 OE_WARN << "No SRS found!" << std::endl;
54 return false;
55 }
56
57 OE_INFO << "Exporting elevation data for extent: "
58 << minX << ", " << minY << " to " << maxX << ", " << maxY
59 << " at LOD " << options.targetLod << std::endl;
60
61 // 計算合適的輸出分辨率
62 unsigned int outputWidth, outputHeight;
63 if (!calculateOutputDimensions(profile, minX, minY, maxX, maxY,
64 options, outputWidth, outputHeight)) {
65 OE_WARN << "Failed to calculate output dimensions!" << std::endl;
66 return false;
67 }
68
69 OE_INFO << "Output dimensions: " << outputWidth << " x " << outputHeight << std::endl;
70
71 // 創建高程數據網格
72 std::vector<float> elevationData(outputWidth * outputHeight, NO_DATA_VALUE);
73
74 // 使用高程池採樣數據（更精確的方法）
75 if (!sampleElevationData(elevationPool, srs, minX, minY, maxX, maxY,
76 outputWidth, outputHeight, options, elevationData)) {
77 OE_WARN << "Failed to sample elevation data!" << std::endl;
78 return false;
79 }
80
81 // 填充數據縫隙
82 if (options.fillGaps) {
83 fillDataGaps(elevationData, outputWidth, outputHeight, options.fillGapsMaxDistance);
84 }
85
86 // 使用GDAL保存為TIFF，包含完整的元數據
87 return saveAsGeoTIFF(outputFilename, elevationData, outputWidth, outputHeight,
88 minX, minY, maxX, maxY, srs, options);
89 }
90
91 private:
92 // 計算輸出尺寸 - 修復版本
93 static bool calculateOutputDimensions(const osgEarth::Profile* profile,
94 double minX, double minY, double maxX, double maxY,
95 const ExportOptions& options,
96 unsigned int& width, unsigned int& height) {
97 try {
98 double extentWidth = maxX - minX;
99 double extentHeight = maxY - minY;
100
101 // 如果用户指定了輸出分辨率，使用它
102 if (options.outputResolution > 0) {
103 width = static_cast<unsigned int>(extentWidth / options.outputResolution) + 1;
104 height = static_cast<unsigned int>(extentHeight / options.outputResolution) + 1;
105 OE_INFO << "Using user-specified resolution: " << options.outputResolution
106 << " meters/pixel" << std::endl;
107 } else {
108 // 基於LOD估算分辨率
109 double resolution = estimateResolutionFromLOD(profile, options.targetLod);
110
111 width = static_cast<unsigned int>(extentWidth / resolution) + 1;
112 height = static_cast<unsigned int>(extentHeight / resolution) + 1;
113
114 OE_INFO << "Using estimated resolution from LOD " << options.targetLod
115 << ": ~" << resolution << " meters/pixel" << std::endl;
116 }
117
118 // 限制最大和最小尺寸
119 width = std::max(256u, std::min(width, 10000u));
120 height = std::max(256u, std::min(height, 10000u));
121
122 OE_INFO << "Calculated dimensions: " << width << " x " << height << std::endl;
123 return true;
124
125 } catch (const std::exception& e) {
126 OE_WARN << "Error calculating dimensions: " << e.what() << std::endl;
127 // 使用保守的默認值
128 width = 1024;
129 height = 1024;
130 return true;
131 }
132 }
133
134 // 根據LOD估算分辨率
135 static double estimateResolutionFromLOD(const osgEarth::Profile* profile, unsigned int lod) {
136 // 地球赤道周長（米）
137 const double EARTH_EQUATORIAL_CIRCUMFERENCE = 40075016.6855785;
138
139 // 基礎瓦片數量（級別0）
140 unsigned int baseTiles = 1;
141
142 // 計算該LOD下的瓦片數量
143 unsigned int numTilesAtLod = baseTiles * (1u << lod); // 2^lod
144
145 // 假設標準瓦片大小
146 unsigned int tileSize = 256;
147
148 // 計算赤道分辨率（米/像素）
149 double equatorialResolution = EARTH_EQUATORIAL_CIRCUMFERENCE / (numTilesAtLod * tileSize);
150
151 // 根據投影類型調整
152 if (profile->getSRS()->isGeographic()) {
153 // 地理座標系，直接使用赤道分辨率
154 return equatorialResolution;
155 } else {
156 // 投影座標系，使用更簡單的估算
157 // 獲取profile的大致範圍來估算
158 osgEarth::GeoExtent profileExtent = profile->getExtent();
159 double profileWidth = profileExtent.width();
160
161 // 假設投影座標單位是米
162 double metersPerDegree = EARTH_EQUATORIAL_CIRCUMFERENCE / 360.0;
163 double profileWidthMeters = profileWidth * metersPerDegree;
164
165 return profileWidthMeters / (numTilesAtLod * tileSize);
166 }
167 }
168
169 // 使用高程池採樣數據（更精確的方法）
170 static bool sampleElevationData(osg::ref_ptr<osgEarth::ElevationPool> elevationPool,
171 const osgEarth::SpatialReference* srs,
172 double minX, double minY, double maxX, double maxY,
173 unsigned int width, unsigned int height,
174 const ExportOptions& options,
175 std::vector<float>& elevationData) {
176 if (!elevationPool) return false;
177
178 // 創建高程包絡線
179 osg::ref_ptr<osgEarth::ElevationEnvelope> envelope =
180 elevationPool->createEnvelope(srs, options.targetLod);
181
182 if (!envelope) {
183 OE_WARN << "Failed to create elevation envelope!" << std::endl;
184 return false;
185 }
186
187 double cellSizeX = (maxX - minX) / (width - 1);
188 double cellSizeY = (maxY - minY) / (height - 1);
189
190 unsigned int totalPoints = width * height;
191 unsigned int pointsPerBatch = 10000; // 分批處理以避免內存問題
192 unsigned int completedPoints = 0;
193
194 OE_INFO << "Sampling " << totalPoints << " elevation points..." << std::endl;
195
196 // 分批處理點
197 for (unsigned int startRow = 0; startRow < height; startRow += pointsPerBatch / width) {
198 unsigned int endRow = std::min(startRow + pointsPerBatch / width, height);
199
200 std::vector<osg::Vec3d> points;
201 std::vector<unsigned int> indices;
202
203 // 準備這一批的點
204 for (unsigned int r = startRow; r < endRow; ++r) {
205 double y = maxY - r * cellSizeY; // 從北向南
206 for (unsigned int c = 0; c < width; ++c) {
207 double x = minX + c * cellSizeX;
208 points.push_back(osg::Vec3d(x, y, 0.0));
209 indices.push_back(r * width + c);
210 }
211 }
212
213 // 批量獲取高程
214 std::vector<float> batchElevations;
215 unsigned int validCount = envelope->getElevations(points, batchElevations);
216
217 // 存儲結果
218 for (size_t i = 0; i < batchElevations.size(); ++i) {
219 elevationData[indices[i]] = batchElevations[i];
220 }
221
222 completedPoints += points.size();
223
224 // 進度報告
225 if (totalPoints > 0) {
226 double progress = (double)completedPoints / totalPoints * 100.0;
227 OE_INFO << std::fixed << std::setprecision(1)
228 << "Progress: " << progress << "% (" << completedPoints
229 << "/" << totalPoints << " points)" << std::endl;
230 }
231
232 if (points.empty()) break;
233 }
234
235 // 統計有效數據點
236 unsigned int validPoints = 0;
237 for (const auto& elev : elevationData) {
238 if (elev != NO_DATA_VALUE) validPoints++;
239 }
240
241 OE_INFO << "Sampling completed: " << validPoints << "/" << totalPoints
242 << " valid elevation points ("
243 << std::fixed << std::setprecision(1)
244 << (double)validPoints / totalPoints * 100.0 << "%)" << std::endl;
245
246 return validPoints > 0;
247 }
248
249 // 填充數據縫隙
250 static void fillDataGaps(std::vector<float>& elevationData,
251 unsigned int width, unsigned int height,
252 double maxFillDistance) {
253 OE_INFO << "Filling data gaps..." << std::endl;
254
255 std::vector<float> filledData = elevationData;
256 unsigned int filledCount = 0;
257
258 // 第一遍：直接鄰居填充
259 for (unsigned int r = 1; r < height - 1; ++r) {
260 for (unsigned int c = 1; c < width - 1; ++c) {
261 unsigned int index = r * width + c;
262
263 if (elevationData[index] == NO_DATA_VALUE) {
264 // 檢查4-鄰域
265 float sum = 0.0f;
266 int count = 0;
267
268 // 上
269 if (elevationData[index - width] != NO_DATA_VALUE) {
270 sum += elevationData[index - width];
271 count++;
272 }
273 // 下
274 if (elevationData[index + width] != NO_DATA_VALUE) {
275 sum += elevationData[index + width];
276 count++;
277 }
278 // 左
279 if (elevationData[index - 1] != NO_DATA_VALUE) {
280 sum += elevationData[index - 1];
281 count++;
282 }
283 // 右
284 if (elevationData[index + 1] != NO_DATA_VALUE) {
285 sum += elevationData[index + 1];
286 count++;
287 }
288
289 if (count > 0) {
290 filledData[index] = sum / count;
291 filledCount++;
292 }
293 }
294 }
295 }
296
297 // 第二遍：擴大搜索範圍填充剩餘的縫隙
298 for (unsigned int r = 0; r < height; ++r) {
299 for (unsigned int c = 0; c < width; ++c) {
300 unsigned int index = r * width + c;
301
302 if (filledData[index] == NO_DATA_VALUE) {
303 float filledValue = interpolateFromNeighbors(elevationData, width, height, c, r, 3);
304 if (filledValue != NO_DATA_VALUE) {
305 filledData[index] = filledValue;
306 filledCount++;
307 }
308 }
309 }
310 }
311
312 elevationData = std::move(filledData);
313 OE_INFO << "Filled " << filledCount << " data gaps" << std::endl;
314 }
315
316 // 從鄰居像素插值
317 static float interpolateFromNeighbors(const std::vector<float>& data,
318 unsigned int width, unsigned int height,
319 unsigned int x, unsigned int y,
320 unsigned int searchRadius) {
321 double sum = 0.0;
322 unsigned int count = 0;
323
324 int startY = std::max(0, static_cast<int>(y) - static_cast<int>(searchRadius));
325 int endY = std::min(static_cast<int>(height) - 1, static_cast<int>(y) + static_cast<int>(searchRadius));
326 int startX = std::max(0, static_cast<int>(x) - static_cast<int>(searchRadius));
327 int endX = std::min(static_cast<int>(width) - 1, static_cast<int>(x) + static_cast<int>(searchRadius));
328
329 for (int ny = startY; ny <= endY; ++ny) {
330 for (int nx = startX; nx <= endX; ++nx) {
331 unsigned int neighborIndex = ny * width + nx;
332 float neighborValue = data[neighborIndex];
333
334 if (neighborValue != NO_DATA_VALUE) {
335 // 使用距離加權
336 double dx = nx - x;
337 double dy = ny - y;
338 double distance = sqrt(dx*dx + dy*dy);
339 double weight = 1.0 / (distance + 1.0); // 避免除零
340
341 sum += neighborValue * weight;
342 count++;
343 }
344 }
345 }
346
347 if (count > 0) {
348 return static_cast<float>(sum / count);
349 }
350
351 return NO_DATA_VALUE;
352 }
353
354 // 使用GDAL保存為GeoTIFF，包含完整元數據
355 static bool saveAsGeoTIFF(const std::string& filename,
356 const std::vector<float>& elevationData,
357 unsigned int width, unsigned int height,
358 double minX, double minY, double maxX, double maxY,
359 const osgEarth::SpatialReference* srs,
360 const ExportOptions& options) {
361 GDALAllRegister();
362
363 // 獲取GTiff驅動
364 GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
365 if (!driver) {
366 OE_WARN << "GTiff driver not available!" << std::endl;
367 return false;
368 }
369
370 // 創建數據集選項
371 char** optionsArray = NULL;
372 optionsArray = CSLSetNameValue(optionsArray, "COMPRESS", options.compression.c_str());
373 optionsArray = CSLSetNameValue(optionsArray, "PREDICTOR", "3"); // 浮點數據用預測器3
374 optionsArray = CSLSetNameValue(optionsArray, "ZLEVEL", "6");
375 optionsArray = CSLSetNameValue(optionsArray, "TILED", "YES");
376 optionsArray = CSLSetNameValue(optionsArray, "BLOCKXSIZE", "256");
377 optionsArray = CSLSetNameValue(optionsArray, "BLOCKYSIZE", "256");
378 optionsArray = CSLSetNameValue(optionsArray, "BIGTIFF", "IF_SAFER");
379
380 // 創建數據集
381 GDALDataset* dataset = driver->Create(filename.c_str(), width, height, 1, GDT_Float32, optionsArray);
382 CSLDestroy(optionsArray);
383
384 if (!dataset) {
385 OE_WARN << "Failed to create TIFF file: " << filename << std::endl;
386 return false;
387 }
388
389 // 設置地理變換參數
390 double pixelSizeX = (maxX - minX) / width;
391 double pixelSizeY = (maxY - minY) / height;
392
393 // 注意：GDAL使用像素中心參考，所以需要調整
394 double geoTransform[6] = {
395 minX + pixelSizeX / 2.0, // 左上角X座標（像素中心）
396 pixelSizeX, // X方向像素大小
397 0, // 旋轉參數（通常為0）
398 maxY - pixelSizeY / 2.0, // 左上角Y座標（像素中心）
399 0, // 旋轉參數（通常為0）
400 -pixelSizeY // Y方向像素大小（負值因為Y從北向南）
401 };
402
403 dataset->SetGeoTransform(geoTransform);
404
405 // 設置投影
406 std::string projWKT = srs->getWKT();
407 if (projWKT.empty()) {
408 // 如果WKT為空，嘗試獲取Proj4字符串
409 std::string proj4 = srs->getHorizInitString();
410 if (!proj4.empty()) {
411 OGRSpatialReference oSRS;
412 if (oSRS.importFromProj4(proj4.c_str()) == OGRERR_NONE) {
413 char* wkt = NULL;
414 oSRS.exportToWkt(&wkt);
415 if (wkt) {
416 projWKT = wkt;
417 CPLFree(wkt);
418 }
419 }
420 }
421 }
422
423 if (!projWKT.empty()) {
424 dataset->SetProjection(projWKT.c_str());
425 } else {
426 OE_WARN << "Could not determine projection WKT!" << std::endl;
427 }
428
429 // 獲取波段並寫入數據
430 GDALRasterBand* band = dataset->GetRasterBand(1);
431
432 // 設置無數據值
433 band->SetNoDataValue(NO_DATA_VALUE);
434
435 // 寫入高程數據
436 CPLErr err = band->RasterIO(GF_Write, 0, 0, width, height,
437 (void*)elevationData.data(), width, height, GDT_Float32,
438 sizeof(float), sizeof(float) * width, NULL);
439
440 if (err != CE_None) {
441 OE_WARN << "Error writing raster data: " << CPLGetLastErrorMsg() << std::endl;
442 GDALClose(dataset);
443 return false;
444 }
445
446 // 設置波段描述
447 band->SetDescription("Elevation");
448 band->SetUnitType("meters");
449
450 // 計算統計信息
451 OE_INFO << "Computing statistics..." << std::endl;
452
453 // 手動計算統計信息，避免GDAL計算時的問題
454 float minVal = FLT_MAX;
455 float maxVal = -FLT_MAX;
456 double sum = 0.0;
457 int count = 0;
458
459 for (const auto& val : elevationData) {
460 if (val != NO_DATA_VALUE) {
461 minVal = std::min(minVal, val);
462 maxVal = std::max(maxVal, val);
463 sum += val;
464 count++;
465 }
466 }
467
468 if (count > 0) {
469 double mean = sum / count;
470 double stddev = 0.0;
471
472 // 計算標準差
473 for (const auto& val : elevationData) {
474 if (val != NO_DATA_VALUE) {
475 double diff = val - mean;
476 stddev += diff * diff;
477 }
478 }
479 stddev = sqrt(stddev / count);
480
481 band->SetStatistics(minVal, maxVal, mean, stddev);
482 OE_INFO << "Elevation statistics - Min: " << minVal << ", Max: " << maxVal
483 << ", Mean: " << mean << ", StdDev: " << stddev << std::endl;
484 } else {
485 OE_WARN << "No valid elevation data for statistics!" << std::endl;
486 }
487
488 // 設置元數據
489 dataset->SetMetadataItem("AREA_OR_POINT", "Point");
490 dataset->SetMetadataItem("TIFFTAG_SOFTWARE", "osgEarth Elevation Exporter");
491 dataset->SetMetadataItem("TIFFTAG_IMAGEDESCRIPTION", "Digital Elevation Model");
492
493 band->SetMetadataItem("ELEVATION_UNITS", "meters");
494 if (count > 0) {
495 band->SetMetadataItem("ELEVATION_MIN", CPLString().Printf("%.3f", minVal));
496 band->SetMetadataItem("ELEVATION_MAX", CPLString().Printf("%.3f", maxVal));
497 }
498 band->SetMetadataItem("NO_DATA_VALUE", CPLString().Printf("%f", NO_DATA_VALUE));
499
500 // 刷新緩存並關閉
501 GDALClose(dataset);
502
503 OE_INFO << "Successfully exported elevation data to: " << filename << std::endl;
504
505 // 驗證輸出文件
506 return validateOutputFile(filename);
507 }
508
509 // 驗證輸出文件
510 static bool validateOutputFile(const std::string& filename) {
511 GDALDataset* ds = (GDALDataset*)GDALOpen(filename.c_str(), GA_ReadOnly);
512 if (!ds) {
513 OE_WARN << "Failed to validate output file: " << filename << std::endl;
514 return false;
515 }
516
517 int width = ds->GetRasterXSize();
518 int height = ds->GetRasterYSize();
519 int bandCount = ds->GetRasterCount();
520
521 OE_INFO << "Output validation - Size: " << width << " x " << height
522 << ", Bands: " << bandCount << std::endl;
523
524 if (bandCount > 0) {
525 GDALRasterBand* band = ds->GetRasterBand(1);
526 int hasNoData = 0;
527 double noDataValue = band->GetNoDataValue(&hasNoData);
528
529 OE_INFO << "NoData value: " << (hasNoData ? std::to_string(noDataValue) : "Not set") << std::endl;
530
531 // 檢查投影
532 const char* projection = ds->GetProjectionRef();
533 if (projection && strlen(projection) > 0) {
534 OE_INFO << "Projection: [Set]" << std::endl;
535 } else {
536 OE_WARN << "Projection: [Missing]" << std::endl;
537 }
538
539 // 檢查地理變換
540 double geoTransform[6];
541 if (ds->GetGeoTransform(geoTransform) == CE_None) {
542 OE_INFO << "GeoTransform: ["
543 << geoTransform[0] << ", " << geoTransform[1] << ", " << geoTransform[2] << ", "
544 << geoTransform[3] << ", " << geoTransform[4] << ", " << geoTransform[5] << "]" << std::endl;
545 } else {
546 OE_WARN << "GeoTransform: [Missing]" << std::endl;
547 }
548 }
549
550 GDALClose(ds);
551 return true;
552 }
553 };
554
555 // 使用示例：
556 void exampleUsage(osgEarth::Map* map) {
557 ElevationDataExporter::ExportOptions options;
558 options.targetLod = 14;
559 options.maxFallbackLevels = 8;
560 options.fillGaps = true;
561 options.fillGapsMaxDistance = 2000.0;
562 options.compression = "DEFLATE";
563 options.outputResolution = 30.0; // 30米/像素
564
565 double minLon = 116.0;
566 double minLat = 39.0;
567 double maxLon = 117.0;
568 double maxLat = 40.0;
569 std::string outputFile = "elevation_export.tif";
570 //DP_ENGINE_PTR->GetKernelData()->GetCurMap()
571 bool success = ElevationDataExporter::exportElevationToTiff(
572 map, minLon, minLat, maxLon, maxLat, outputFile, options);
573
574 if (success) {
575 OE_INFO << "Elevation data exported successfully to: " << outputFile << std::endl;
576 } else {
577 OE_WARN << "Failed to export elevation data!" << std::endl;
578 }
579 }
主要修復和改進：
- 移除不存在的函數：刪除了
getEquatorialResolution函數調用 - 新的分辨率估算方法：實現了
estimateResolutionFromLOD函數，基於地球赤道周長和LOD來估算分辨率 - 支持用户指定分辨率：添加了
outputResolution選項，讓用户可以明確指定輸出分辨率 - 改進的地理變換：修正了GDAL地理變換參數，使用像素中心參考
- 手動統計計算：改進了統計信息計算，避免GDAL內部計算可能的問題
- 更好的縫隙填充：使用兩遍填充策略，先4-鄰域快速填充，再擴大範圍填充剩餘縫隙
- 更詳細的驗證：增加了對地理變換參數的驗證輸出
這個版本應該能夠正確處理高程數據的導出，避免拼接縫隙，並提供完整的高程描述信息。
本文章為轉載內容，我們尊重原作者對文章享有的著作權。如有內容錯誤或侵權問題，歡迎原作者聯繫我們進行內容更正或刪除文章。