文件格式以及原理參考老外的文章：https://librenepal.com/article/reading-srtm-data-with-python/

?

GDAL包：https://github.com/OSGeo/gdal/releases/tag/v3.1.3

其實文件格式很容易理解，比如是1/3弧度的精度情況下，就是1度分為1/1200份，所以一個文件表示經度和維度各1度的方格，就是切成1200x1200份，存儲為二維矩陣是1201x1201，因為邊界占了1行1列：

其中二維數組是按照地圖來存儲的，所以從低維度和低經度取索引時候需要計算下：

老外是用python寫的，我用c++重寫的，歌詞大意如下：

#pragma once#include <math.h> #include <string.h> #include <stdio.h> #include <string> #include <math.h> #include <iostream> #include <memory> #include <algorithm> #include <map> #include <mutex> #include <thread>//#include "include/gdal.h" #include <gdal_priv.h>#ifdef _DEBUG#pragma comment(lib, "gdal_i_d.lib") #else #pragma comment(lib, "gdal_i.lib") #endifusing namespace std; // 數據結構 class SRTM { public:SRTM(){if (runOnce == false){GDALAllRegister();// windows操作系統使用GBKCPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");runOnce = true;}}~SRTM(){if (pData != nullptr)delete[] pData;//std::cerr << "SRTM析構釋放數據" << endl;};// 分辨率const int sample = 1200;static bool runOnce;private:SRTM(const SRTM & other) = delete;SRTM & operator = (const SRTM & other) = delete;SRTM(SRTM && other) = delete;SRTM & operator = (SRTM && other) = delete;unsigned int nWidth = 1201;unsigned int nHeight = 1201;// c++17以前不支持動態數組使用shared_ptr管理// std::unique_ptr<short[]> data;short int *pData = nullptr;double minLat;double maxLat;double minLon;double maxLon;static double Mercator2Lon(double lon)//墨卡托轉WGS84:經度{return lon / 20037508.34 * 180.0;}static double Mercator2Lat(double lat)//墨卡托轉WGS84:緯度{double result = 0;double mid = lat / 20037508.34 * 180.0;result = 180.0 / M_PI * (2.0 * atan(exp(mid * M_PI / 180.0)) - M_PI / 2.0);return result;}public:// 從二維數組中查詢，但是編號是二維數組的左下角，需要重新計算一下indexbool query(short & height, double lat, double lon){if (pData == nullptr)return false;int lat_row = int(round((lat - int(lat)) * sample));int lon_row = int(round((lon - int(lon)) * sample));//lat_row = int(round((lat - minLat) * sample));//lon_row = int(round((lon - minLon) * sample));lat_row = abs(lat_row);lon_row = abs(lon_row);size_t index = (nHeight - 1 - lat_row) * nWidth + lon_row;if (index >= sample * sample)return false;height = pData[index];return true;}bool load(const std::string fileName){return load(fileName.c_str());}bool load(const char * fileName){GDALDataset *poDataSet;GDALRasterBand *pBand;poDataSet = (GDALDataset*)GDALOpen(fileName, GA_ReadOnly);if (poDataSet == nullptr)return false;this->nWidth = poDataSet->GetRasterXSize();//獲取圖像寬度this->nHeight = poDataSet->GetRasterYSize();//獲取圖像高度// 存儲邊界信息double adfGeoTransform[6];double value[6];if (poDataSet->GetGeoTransform(adfGeoTransform) == CE_None){value[0] = adfGeoTransform[0]; // 起點，左上經度value[1] = adfGeoTransform[3]; // 起點，左上維度value[2] = adfGeoTransform[1] * (double)nWidth + adfGeoTransform[0]; // 右側經度value[3] = adfGeoTransform[5] * (double)nHeight + adfGeoTransform[3]; // 右下if (value[0] > 180 || value[0] < -180)//墨卡托轉WGS84{value[0] = Mercator2Lon(value[0]);value[1] = Mercator2Lat(value[1]);value[2] = Mercator2Lon(value[2]);value[3] = Mercator2Lat(value[3]);}}this->minLon = value[0];this->maxLon = value[2];this->minLat = value[3];this->maxLat = value[1];if (pData != nullptr)delete[] pData;this->pData = new short[nWidth * nHeight];pBand = poDataSet->GetRasterBand(1);pBand->RasterIO(GF_Read, 0, 0, nWidth, nHeight, pData, nWidth, nHeight,pBand->GetRasterDataType(), 0, 0);//int i = pData[1000 * nWidth + 1];GDALClose(poDataSet);//關閉數據集return true;}public:// 通過經緯度計算標準文件名static std::string getFileName(double lat, double lon){char ns;char ew;if (lat >= 0)ns = 'N';elsens = 'S';if (lon >= 0)ew = 'E';elseew = 'W';char buffer[20];int i_lat = abs(int(lat));int i_lon = abs(int(lon));//snprintf(buffer, 20, "%.1s", &ns);snprintf(buffer, 20, "%.1s%02d%.1s%03d.hgt", &ns, i_lat, &ew, i_lon);return buffer;}}; // 初始化 bool SRTM::runOnce = false;// c++ 17 //namespace fs = std::filesystem; // 對緩存進行管理 class SRTM_Cache { public:SRTM_Cache(const char * dir) : rootDir(dir){setRootPath(dir);}~SRTM_Cache(){// 自動釋放資源//dataMap.clear();}void setRootPath(const char * dir){rootDir = dir;// 去掉末尾的\\ size_t off = rootDir.rfind('\\', rootDir.size());if (off > 0 && (off == rootDir.size()-1))rootDir = rootDir.substr(0, off);}bool query(short & height, double lat, double lon){std::string fileName = SRTM::getFileName(lat, lon);bool ret = false;std::shared_ptr<SRTM> ptr = nullptr;{ // 添加作用域，提早解鎖std::lock_guard<std::mutex> autoLock(mapMutex); // 加鎖auto it = dataMap.find(fileName);if (it == dataMap.end()){// 嘗試加載文件std::string filePath = rootDir + "\\" + fileName;ptr = std::make_shared<SRTM>();ret = ptr->load(filePath);if (ret){dataMap.insert(std::pair<std::string, std::shared_ptr<SRTM> >(fileName, ptr));}else{return false;}return ret;}else{std::shared_ptr<SRTM> ptr = it->second;}} // 添加作用域，提早解鎖//std::cout << ptr.use_count() << endl;if (ptr != nullptr){ret = ptr->query(height, lat, lon);return ret;}return false;} // end of queryprivate:SRTM_Cache(const SRTM_Cache & other) = delete;SRTM_Cache& operator =(const SRTM_Cache & other) = delete;// 用智能指針管理數據類實例，因為已經禁止了賦值和拷貝構造，std::map<std::string, std::shared_ptr<SRTM> > dataMap;std::string rootDir;// 添加多線程互斥std::mutex mapMutex; };

使用的方法如下：

// readDEM.cpp : 此文件包含 "main" 函數。程序執行將在此處開始并結束。 //#include <iostream> #include "SRTM.h"int main() {SRTM_Cache manager("D:\\DEM數據\\SRTM3-90米全國DEM\\");double lat = 39.990618;double lon = 116.169644;short heiht;bool ret = manager.query(heiht, lat, lon);if (ret == false){std::cout << "false" << endl;}else{string str;str.resize(100, '\0');snprintf(const_cast<char *>(str.data()), 100, "%.5f, %.5f 高程：%d", lat, lon, heiht);std::cout << str.c_str() << endl;}lat = 41.56;ret = manager.query(heiht, lat, lon);if (ret == false){std::cout << "false" << endl;}else{string str;str.resize(100, '\0');snprintf(const_cast<char *>(str.data()), 100, "%.5f, %.5f 高程：%d", lat, lon, heiht);std::cout << str.c_str() << endl;} }

結果：

39.99062, 116.16964 高程：509 41.56000, 116.16964 高程：1792

備注：

使用vcpkg管理各種開源包真的非常的方便，比自己一個一個找強多了。

總結

以上是生活随笔為你收集整理的使用GDAL库读取SRTM格式的高程数据的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。