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

头文件：

//#include "include/gdal.h"  
#include <gdal_priv.h>

#include <math.h>


#ifdef _DEBUG
#define new DEBUG_NEW
#pragma comment(lib, "gdal_i_d.lib")
#else
#pragma comment(lib, "gdal_i.lib")
#endif

函数部分：

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;
}

// 开始计算
void CreadSRTMDlg::OnCompute()
{
	// TODO: Add your command handler code here
	CString filename;
	filename = "srtm_66_21.tif";

	GDALAllRegister();
	GDALDataset *poDataSet;
	GDALRasterBand *pBand;
	int nWidth, nHeight;
	poDataSet = (GDALDataset*)GDALOpen((LPCTSTR)filename, GA_ReadOnly);

	nWidth = poDataSet->GetRasterXSize();//获取图像宽度
	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]);
		}
	}

	
	short int *pData = new short int[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];
	CPLFree(pData);//释放内存
	GDALClose(poDataSet);//关闭数据集
}

我的版本是3.1.3，涉及的库比较多，我使用vcpkg管理，已经导出并放到资源里。