使用GDAL库读取SRTM格式的高程数据

文件格式以及原理参考老外的文章：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格式的高程数据相关推荐

python 读取geotiff_科学网—利用python GDAL库读写geotiff格式的遥感影像方法 - 张伟的博文...
(1)利用python GDAL库读写geotiff格式的遥感影像方法,具有很好的参考价值,不错! from osgeo import gdal import numpy as np def read ...
python读取tiff影像_科学网—利用python GDAL库读写geotiff格式的遥感影像方法 - 张伟的博文...
(1)利用python GDAL库读写geotiff格式的遥感影像方法,具有很好的参考价值,不错! from osgeo import gdal import numpy as np def read ...
GDAL库——读取图像并提取基本信息
GDAL库是一个跨平台的栅格地理数据格式库,包括读取.写入.转换.处理各种栅格数据格式(有些特定的格式对一些操作如写入等不支持).它使用了一个单一的抽象数据模型就支持了大多数的栅格数据.这里有GDAL ...
netcdf库读取nc格式文件中的字符串类型的数据
netcdf库读取nc格式文件一.背景二.工具使用三.测试代码四.测试结果一.背景这两天解析数据需要解析nc格式的文件,时间是字符串类型的,最开始还以为这个数据有问题呢,使用panoply ...
MATLAB读取HDF格式的SST数据
利用MATLAB读取HDF格式的SST数据是比较简单的,MATLAB中有专门用于读取HDF的函数hdfread()与hdfinfo()就能很好的读取HDF文件.我们可以在读取文件的时候,通过显示文件中 ...
ajax读取文件数据,Ajax 实现读取 properties 格式资源文件数据
Ajax 的核心是 JavaScript 对象 XmlHttpRequest.该对象在 Internet Explorer 5 中首次引入,它是一种支持异步请求的技术.简而言之,XmlHttpRequ ...
arcgis风向_ArcGIS10.2读取NetCDF格式的气象数据含风向
ArcGIS读取NetCDF格式的气象数据尝试使用Make netCDF File的Tool 参数设置: 变量:PM2_5DRY X维度:XLONG Y维度:XLAT Dimension:南北+东西 ...
编译好的GDAL库,支持ECW格式,支持proj,支持geos
这几天因为工作需要,把GDAL重新编译了一下,现支持ECW格式,即可以用GDALOPEN读取JPEG2000格式的影像数据;编译时也添加了对geos和proj的支持. 内含3个文件夹,分别是gdal, ...
GDAL库读取Envisat ASAR数据
GDAL库本身就可以读取Envisat的图像数据,具体链接为:http://www.gdal.org/frmt_various.html#Envisat. 但是对于ASAR传感器的数据来说,GDAL在 ...
Python读取.edf格式脑电数据文件
MNE-python读取.edf文件 EDF,全称是 European Data Format,是一种标准文件格式,用于交换和存储医疗时间序列. 该格式文件能够存储多通道的数据,允许每个信号拥有不同的 ...

使用GDAL库读取SRTM格式的高程数据

使用GDAL库读取SRTM格式的高程数据相关推荐

最新文章

热门文章