使用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")
#endif

using namespace std;
// 数据结构
class SRTM
{
public:
	SRTM()
	{
		if (runOnce == false)
		{
			GDALAllRegister();
			
			// windows操作系统使用GBK
			CPLSetConfigOption("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:
	// 从二维数组中查询，但是编号是二维数组的左下角，需要重新计算一下index
	bool 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';
		else
			ns = 'S';

		if (lon >= 0)
			ew = 'E';
		else
			ew = '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 query

private:
	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管理各种开源包真的非常的方便，比自己一个一个找强多了。