FH98ViewOnPlaneSoftware/Src/QB_ImgProcessFun/GetElevation.cpp

#include "StdAfx.h"
#include "GetElevation.h"

// tiff文件读取
#include "gdal_priv.h"  
#include "cpl_conv.h" //for CPLMalloc()

// boost-1.56:准C++标准库
#include <boost/lexical_cast.hpp>
#include <boost/format.hpp>

// 全局变量
std::deque<struDemPic*> g_DemPicList;

// 功能：清理DEM记忆内存
void ClearDemMem()
{
	if (g_DemPicList.empty() == false)
	{
		for (size_t i = 0; i < g_DemPicList.size(); i++)
		{
			if (g_DemPicList.at(i) != nullptr && g_DemPicList.at(i)->pDEMData != nullptr)
			{
				delete[] g_DemPicList.at(i)->pDEMData;
				g_DemPicList.at(i)->pDEMData = nullptr;
				delete g_DemPicList.at(i);
				g_DemPicList.at(i) = nullptr;
			}
		}

		g_DemPicList.clear();
	}
}

// 基于经纬度 计算tif文件名称(图像名称对应该图像所在地图的左下角经纬度)
string  GetDemName( const double dLonDeg, const double dLatDeg )
{
	string sLatLon; // 输出

	int iLon = (int) dLonDeg;
	int iLat = (int) dLatDeg;

	//判断经纬度在正常范围之内
	if (abs(dLonDeg) > 180)
	{
		return "";
	}

	if (abs(dLatDeg) > 90)
	{
		return "";
	}

	string sLat, sLon;
	

	// 纬度 南北纬 区别对待
	boost::format fmt("S%02d");
	if (dLatDeg <= 0)
	{
		fmt % ((-iLat) + 1);
		sLat = fmt.str();
	}
	else
	{
		fmt = boost::format("N%02d");
		fmt % iLat;
		sLat = fmt.str();
	}

	// 经度 东西经 区别对待
	if (dLonDeg <= 0)
	{
		fmt = boost::format("W%03d");
		fmt % ((-iLon) + 1);
		sLon = fmt.str();
	}
	else
	{
		fmt = boost::format("E%03d");
		fmt % iLon;
		sLon = fmt.str();
	}

	sLatLon = "ASTGTM_" + sLat + sLon + "_dem.tif";
	return sLatLon;
}


//  搜寻文件路径是否已经存在
int FindDemDataSet( const string strDemName )
{
	for ( unsigned int i = 0; i < g_DemPicList.size(); i++ )
	{
		if ( g_DemPicList[i]->strDemName == strDemName )
		{
			return i;
		}
	}
	return -1;
}

bool GetDem( float& fElevator, double dLonDeg, double dLatDeg, struDemPic* pDemPic, bool bNearst)
{
	// 空值判断
	if (pDemPic == NULL)
	{
		return false;
	}

	// 分辨率 判断
	if (pDemPic->adfGeoTransform[1] <= 0 || pDemPic->adfGeoTransform[5] >= 0)
	{
		return false;
	}

	// 边界
	double leftLon   = pDemPic->adfGeoTransform[0];
	double rightLon  = leftLon + 1;
	double topLat    = pDemPic->adfGeoTransform[3];
	double bottomLat = topLat - 1;

	// 计算像素坐标
	double px = (dLonDeg - pDemPic->adfGeoTransform[0]) / pDemPic->adfGeoTransform[1];
	double py = (dLatDeg - pDemPic->adfGeoTransform[3]) / pDemPic->adfGeoTransform[5];

	int ipx = int(px);
	int ipy = int(py);

	// 在图像内部，不含下边界和右边界
	if (dLonDeg >= leftLon && dLonDeg < rightLon && dLatDeg > bottomLat && dLatDeg <= topLat)
	{
		// 最近邻法
		if (bNearst == true)
		{
			int x = int(px + 0.5);
			int y = int(py + 0.5);
			fElevator = *(pDemPic->pDEMData + y * pDemPic->iXSize + x);
		}
		else
		{
			// 规则:双线性内插 计算高程值
			int yt = 0, yb = 0; // y 上下
			int	xl = 0, xr = 0; // x 左右

			yt = int(py);
			yb = yt + 1;
			xl = int(px);
			xr = xl + 1;

			// 相邻四个像素的值
			signed short int DNlt = *(pDemPic->pDEMData + yt * pDemPic->iXSize + xl);
			signed short int DNrt = *(pDemPic->pDEMData + yt * pDemPic->iXSize + xr);
			signed short int DNlb = *(pDemPic->pDEMData + yb * pDemPic->iXSize + xl);
			signed short int DNrb = *(pDemPic->pDEMData + yb * pDemPic->iXSize + xr);

			// 间距
			double h, v;
			h = py - yt;
			v = px - xl;

			fElevator = float(DNlt * (1 - h) * (1 - v) + DNlb * h * (1 - v)+ DNrt * (1 - h) * v + DNrb * h * v);
		}
	}
	else
	{
		// 边缘用最近邻法取值
		if ((px < 0) || (px > pDemPic->iXSize - 1) || 
			(py < 0) || ( py > pDemPic->iYSize - 1))
		{
			signed short DNlt = *(pDemPic->pDEMData + ipy * pDemPic->iXSize + ipx);
			fElevator = DNlt;
		}
		else
		{
			return false;
		}
	}

	return true;
}



// 获取高程
bool GetDem( float& fElevator, double dLonDeg, double dLatDeg, string sDir, bool bNearst)
{
	// tif文件 名称
	string strName = GetDemName(dLonDeg, dLatDeg);
	if (strName == "")
	{
		return false;
	}

	// tif文件 路径名称
	string strPathName; 

	if (sDir != "")
	{
		strPathName = sDir + "\\" + strName;
	}
	else
	{
		return false;
	}

	// 是否能在记录中找到
	int idx = FindDemDataSet(strPathName);
	if (idx >= 0)
	{
		return GetDem(fElevator, dLonDeg, dLatDeg, g_DemPicList[idx], bNearst);;
	}

	//如果在数据集数组中找不到 读取文件并添加到数据集
	// 注册
	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");

	GDALDataset* tempdata = NULL;
	tempdata = (GDALDataset*) GDALOpen(strPathName.c_str(), GA_ReadOnly);
	if ( tempdata != NULL )
	{
		GDALDataType dataType = tempdata->GetRasterBand(1)->GetRasterDataType();
		int iXSzie = tempdata->GetRasterBand(1)->GetXSize();
		int iYSzie = tempdata->GetRasterBand(1)->GetYSize();
		struDemPic* temp = new struDemPic();
		temp->iXSize = iXSzie;
		temp->iYSize = iYSzie;
		tempdata->GetGeoTransform(temp->adfGeoTransform);

		temp->pDEMData = new signed short int[iXSzie * iYSzie];

		memset((void*) temp->pDEMData, 0, iXSzie * iYSzie * sizeof(signed short int));
		tempdata->GetRasterBand(1)->RasterIO( GF_Read, 0, 0, iXSzie, iYSzie, temp->pDEMData, iXSzie, iYSzie, GDT_Int16, 0, 0 );
		GDALClose(tempdata);
		tempdata = NULL;

		// 入队 存储
		temp->strDemName = strPathName;
		g_DemPicList.push_back(temp);

		if (g_DemPicList.size() > 9) // 通常一次分析区域不会超过3*3度
		{
			delete [] g_DemPicList[0]->pDEMData;
			g_DemPicList[0]->pDEMData = nullptr;
			delete g_DemPicList[0];
			g_DemPicList[0] = nullptr;
			g_DemPicList.pop_front();
		}

		return GetDem(fElevator, dLonDeg, dLatDeg, g_DemPicList[g_DemPicList.size() - 1], bNearst);;
	}
	else
	{
		return false;
	}
}