GDAL获取指定地理坐标像元值、获取指定地理范围影像数据

标签：gdal   gis   数据   algorithm   

//GdalImage.h
#include "StructDef.h"

#include "gdal1.11.2/gdal_priv.h"
#include "gdal1.11.2/gdal.h"
//#ifdef __cplusplus
//extern "C" {
//#endif


struct stRasterInfo
{
	char fileName[255];
	GDALDataset* pDataset;
	XRECT<rtsDataTypeGeo> rasterRange;
	int nBandCount;
};

struct stBitmapInfo
{
	char fileName[255];
	bool isInters;//是否有相交部分
	BYTE* pBuf;
	int nBandCount;

	int xStart;
	int yStart;
	int xEnd;
	int yEnd;

	int xSize;
	int ySize;
	stBitmapInfo():isInters(false), pBuf(NULL), nBandCount(0), xStart(0),yStart(0),xEnd(0),yEnd(0), xSize(0), ySize(0){}
	~stBitmapInfo()
	{
		if( pBuf )
		{
			free( pBuf );
			pBuf = NULL;
		}			
	}
};


static char g_szRasterPath[255] = {0};	//栅格影像所在文件夹全路径路径
static int g_screenWidth = 0;	//屏幕分辨率，宽
static int g_screenHeight = 0;	//屏幕分辨率，高



//由地理坐标得到图像行列号
bool Projection2ImageRowCol(double *adfGeoTransform, double dProjX, double dProjY, int &iCol, int &iRow);

//根据地理位置获得指定点的像元值2015.3.5
//cszRasterPath[IN]:影像文件绝对路径	GeoX, GeoY[IN]:地理位置	pBuf[OUT]:存放读取到的像元值，为NULL读取失败
//nBandCOunt[OUT]:存放波段数	xSize, ySize[IN]:要获取像元值的范围，宽和高默认为1
void GetPixValByGeoPos( const char *cszRasterPath, const double GeoX, const double GeoY,
	int *&pBuf, int &nBandCount, int xSize, int ySize);// = 1  = 1

//初始化路径及分辨率
void InitVariable( const char *szRasterPath, const int screenWidth, const int screenHeight );

//获取影像对象
void GetBuf( const XRECT<rtsDataTypeGeo> showRange, stBitmapInfo *&pBitmap, int& nSize );

//#ifdef __cplusplus
//}
//#endif

//  GdalImage.cpp
#include <vector>
#include <dirent.h>
//#include <io.h>
#include <string.h>
#include "StructDef.h"
#include "DRECT.h"

#include "GdalImage.h"
#include "gdal1.11.2/gdal_priv.h"
#include "gdal1.11.2/gdal.h"
#include <iostream>
#include "mapWndBase.h"

//const int DirCharMAXLEN = 1024; //定义最大目录长度
static vector<stRasterInfo> g_vecRI;	//初始化时存放指定文件夹下所有tif,bmp,img信息

//由地理坐标得到图像行列号
bool Projection2ImageRowCol(double *adfGeoTransform, double dProjX, double dProjY, int &iCol, int &iRow)
{
	try
	{
		double dTemp = adfGeoTransform[1]*adfGeoTransform[5] - adfGeoTransform[2]*adfGeoTransform[4];
		double dCol = 0.0, dRow = 0.0;
		dCol = (adfGeoTransform[5]*(dProjX - adfGeoTransform[0]) -
			adfGeoTransform[2]*(dProjY - adfGeoTransform[3])) / dTemp + 0.5;
		dRow = (adfGeoTransform[1]*(dProjY - adfGeoTransform[3]) -
			adfGeoTransform[4]*(dProjX - adfGeoTransform[0])) / dTemp + 0.5;

		iCol = int(dCol);
		iRow = int(dRow);
		return true;
	}
	catch(...)
	{
		return false;
	}
}
//根据地理位置获得指定点的像元值2015.3.5
//cszRasterPath[IN]:影像文件绝对路径	GeoX, GeoY[IN]:地理位置	pBuf[OUT]:存放读取到的像元值，为NULL读取失败
//nBandCOunt[OUT]:存放波段数	xSize, ySize[IN]:要获取像元值的范围，宽和高默认为1
void GetPixValByGeoPos( const char *cszRasterPath, const double GeoX, const double GeoY,
	int *&pBuf, int &nBandCount, int xSize = 1, int ySize = 1 )
{
	try
	{
		//GDALAllRegister();
		GDALDataset *pDataset = (GDALDataset*)GDALOpen( cszRasterPath, GA_ReadOnly );
		double GT[6];
		//得到仿射变换模型
		pDataset->GetGeoTransform( GT );
		//得到影像宽和高
		int rasterXSize = pDataset->GetRasterXSize();
		int rasterYSize = pDataset->GetRasterYSize();
		//得到影像波段数
		nBandCount = pDataset->GetRasterCount();
		int iCol, iRow;
		//地理位置转换为影像行列位置
		if( !Projection2ImageRowCol( GT, GeoY, GeoX, iCol, iRow ) )
			throw 0;
		//地理位置转换结果是否超出影像范围
		if( iCol < 0 || iCol > rasterXSize || iRow < 0 || iRow > rasterYSize )
		{
			pBuf = NULL;
		}
		else
		{
			pBuf = new int[nBandCount*xSize*ySize];
			int *panBandMap= new int[nBandCount];
			for( int i = 0; i < nBandCount; ++i )
				panBandMap[i] = i+1;
			//获取指定行列位置每个波段的像元值，存入pBuf
			pDataset->RasterIO( GF_Read, iCol, iRow, xSize, ySize, pBuf, xSize,
				ySize, GDT_CInt32, nBandCount, panBandMap, 4, 0, 0 );

			delete[] panBandMap;
			panBandMap = NULL;
		}
	}
	catch(...)
	{
		if( pBuf )
			delete[] pBuf;
		pBuf = NULL;
	}

	return;
}

//Linux下使用
void GetFileList( const char* szPath, vector<string>& vecFile )
{
	DIR* dp = opendir (szPath);
	if (! dp)
		return;

	if (chdir (szPath) == -1)
		 return;

	struct dirent* de;
	for (de = readdir (dp); de; de = readdir (dp))
	{
		if ( (de -> d_type != DT_DIR) )
		{
			string strName(de->d_name);
			string strFormat = strName.substr(strName.length()-3, 3);
			if( "tif"==strFormat || "img"==strFormat || "bmp"==strFormat )
				vecFile.push_back( strName );
		}
	}
	closedir (dp);
	return ;
}
//windows下使用
//void GetFileList( const char* szPath, vector<string>& vecFile )
//{
//	_finddata_t fdata; //定义文件查找结构对象
//	long done;
//	char tempdir[DirCharMAXLEN] = {0}; //定义一个临时字符数组
//	strcat(tempdir, szPath); //连接字符串
//	strcat(tempdir, "\\*.*"); //连接字符串(搜索所有文件)
//	done = _findfirst(tempdir, &fdata); //开始查找文件
//
//	if(done != -1) //是否查找成功
//	{
//		int ret = 0;
//		while(ret != -1) //定义一个循环
//		{
//			if(fdata.attrib != _A_SUBDIR) //判断文件属性
//			{
//				//cout << fdata.name << endl;
//				if(strcmp(fdata.name, "...")!=0 &&
//					strcmp(fdata.name, "..") != 0 &&
//					strcmp(fdata.name, ".") != 0) //过滤
//				{
//					char dir[DirCharMAXLEN] = {0}; //定义字符数组
//					strcat(dir, szPath); //连接字符串
//					strcat(dir, "\\");
//					strcat(dir, fdata.name);
//					vecFile.push_back(string(dir));
//				}
//			}
//			ret = _findnext(done, &fdata); //查找下一个文件
//		}
//	}
//}

bool ImageRowCol2Projection(double *adfGeoTransform, int iCol, int iRow, double &dProjX, double &dProjY)
{
	//adfGeoTransform[6]  数组adfGeoTransform保存的是仿射变换中的一些参数，分别含义见下
	//adfGeoTransform[0]  左上角x坐标
	//adfGeoTransform[1]  东西方向分辨率
	//adfGeoTransform[2]  旋转角度, 0表示图像 "北方朝上"
	//adfGeoTransform[3]  左上角y坐标
	//adfGeoTransform[4]  旋转角度, 0表示图像 "北方朝上"
	//adfGeoTransform[5]  南北方向分辨率
	try
	{
		dProjX = adfGeoTransform[0] + adfGeoTransform[1] * iCol + adfGeoTransform[2] * iRow;
		dProjY = adfGeoTransform[3] + adfGeoTransform[4] * iCol + adfGeoTransform[5] * iRow;
		return true;
	}
	catch(...)
	{
		return false;
	}
}

bool GetRasterInfo( vector<stRasterInfo>& vecRI )
{
	if( !g_szRasterPath )
		return false;
	vector<string> vecFile;
	GetFileList( g_szRasterPath, vecFile );

	try
	{
		//GDALAllRegister();
		for( int i = 0; i < vecFile.size(); ++i )
		{
			GDALDataset *pDataset = (GDALDataset*)GDALOpen( vecFile[i].c_str(), GA_ReadOnly );
			double GT[6];
			//得到仿射变换模型
			pDataset->GetGeoTransform( GT );
			//得到影像宽和高
			int rasterXSize = pDataset->GetRasterXSize();
			int rasterYSize = pDataset->GetRasterYSize();
			//得到影像波段数
			int nBandCount = pDataset->GetRasterCount();
			GDALClose(pDataset);
			pDataset = NULL;
			double dProjX, dProjY;
			if( !ImageRowCol2Projection(GT, rasterXSize, rasterYSize, dProjX, dProjY) )
				throw 0;

			stRasterInfo RI;
			RI.pDataset = NULL;
			RI.nBandCount = nBandCount;
			strcpy( RI.fileName, vecFile[i].c_str() );
			RI.rasterRange.m_left = GT[0];
			RI.rasterRange.m_right = dProjX;
			RI.rasterRange.m_top = GT[3];
			RI.rasterRange.m_bottom = dProjY;
			vecRI.push_back(RI);
		}
	}
	catch(...)
	{
		return false;
	}
	return true;
}

void InitVariable( const char *szRasterPath, const int screenWidth, const int screenHeight )
{
	if( !szRasterPath || screenWidth <= 0 || screenHeight <= 0 )
		return;
	strcpy( g_szRasterPath, szRasterPath );
	g_screenWidth = screenWidth;
	g_screenHeight = screenHeight;
	if( !g_vecRI.empty() )
 		g_vecRI.clear();
 	GetRasterInfo( g_vecRI );
	return;
}

bool GetBitmapInfo( stRasterInfo stRI, int* panBandMap, const XRECT<rtsDataTypeGeo>& showRange,
	stBitmapInfo& stBitmap )
{
	int nBandCount;
	if( 7 == stRI.nBandCount )
		nBandCount = 3;
	else
		nBandCount = stRI.nBandCount;
	
	XRECT<rtsDataTypeGeo> rasterRange = stRI.rasterRange;
	try
	{
		//是否有相交范围，如果有取相交部分
		XRECT<rtsDataTypeGeo> interRect;
		if( !interRect.IntersectRect(showRange, rasterRange) )
		{
			strcpy( stBitmap.fileName, stRI.fileName );
			stBitmap.isInters = false;
			if( stRI.pDataset )
			{
				GDALClose( stRI.pDataset );
				stRI.pDataset = NULL;
			}
			return true;
		}
		if( !stRI.pDataset )
			stRI.pDataset = (GDALDataset*)GDALOpen( stRI.fileName, GA_ReadOnly );
		GDALDataset *pDataset = stRI.pDataset;
		//得到仿射变换模型
		double GT[6];
		pDataset->GetGeoTransform( GT );
		int iCol_L=0, iRow_B=0, iCol_R=0, iRow_T=0;//左，下，右，上 像元行列号
		if( !Projection2ImageRowCol(GT, interRect.m_left, interRect.m_bottom, iCol_L, iRow_B) )
			throw 0;
		if( !Projection2ImageRowCol(GT, interRect.m_right, interRect.m_top, iCol_R, iRow_T) )
			throw 0;
		//得到行列宽，高
		int xSize = abs(iCol_R-iCol_L);
		int ySize = abs(iRow_B-iRow_T);
		//缓冲区大小
		int const bufLen = nBandCount*g_screenWidth*g_screenHeight;
		BYTE *pBuf = new BYTE[bufLen];

		pDataset->RasterIO( GF_Read, iCol_L, iRow_T, xSize, ySize, pBuf, g_screenWidth,
			g_screenHeight, GDT_Byte, nBandCount, panBandMap, nBandCount, 0, 1 );

		//整理pBuf
		BYTE *pBufBack;
		if( 1 == nBandCount )//单波段
		{
			//pBufBack = new int[bufLen*4];
			pBufBack = (BYTE*)malloc( sizeof(int)*bufLen*4 );
			for( int i = 0; i < bufLen; ++i )
			{
				pBufBack[i*4] = pBuf[i];	pBufBack[i*4+1] = pBuf[i];
				pBufBack[i*4+2] = pBuf[i];	pBufBack[i*4+3] = 0;
			}
		}
		else if( 3==nBandCount || 7==nBandCount )//3波段
		{
			//pBufBack = new int[bufLen+bufLen/3];
			pBufBack = (BYTE*)malloc(sizeof(BYTE)*(bufLen+bufLen/3));
			for( int i = 0; i < bufLen/3; ++i )
			{
				pBufBack[i*4] = pBuf[i*3];	pBufBack[i*4+1] = pBuf[i*3+1];
				pBufBack[i*4+2] = pBuf[i*3+2];	pBufBack[i*4+3] = 0;
			}
		}
		else if( 4 == nBandCount )//4波段
		{
			//pBufBack = new int[bufLen];
			pBufBack = (BYTE*)malloc(bufLen*sizeof(int));
			for( int i = 0; i < bufLen; ++i )
			{
				if( 3 == (i%4) )
					pBufBack[i] = 0;
			}
		}
		delete[] pBuf;

		XRECT<rtsDataTypeDevice> bmpDevice;
		GeoToDevice(interRect,bmpDevice);

		strcpy( stBitmap.fileName, stRI.fileName );
		stBitmap.isInters = true;
		stBitmap.nBandCount = nBandCount;
		stBitmap.pBuf = pBufBack;

		stBitmap.xStart=bmpDevice.m_left;
		stBitmap.yStart=bmpDevice.m_top;
		stBitmap.xEnd=bmpDevice.m_right;
		stBitmap.yEnd=bmpDevice.m_bottom;

		stBitmap.xSize = g_screenWidth;
		stBitmap.ySize = g_screenHeight;
	}
	catch(...)
	{
		return false;
	}
	return true;
}

void DeleteBmpArray( stBitmapInfo *&pBitmap )
{
	delete[] pBitmap;
}

void GetBuf( const XRECT<rtsDataTypeGeo> showRange, stBitmapInfo *&pBitmap, int& nBitmap )
{
	nBitmap = 0;
	int nBandCount;//波段
	int *panBandMap;
	int nVecRI = g_vecRI.size();
 	pBitmap = new stBitmapInfo[nVecRI];
	for( int i = 0; i < nVecRI; ++i )
	{
		nBandCount = g_vecRI[i].nBandCount;
		//1波段 或 3波段
		if( 1==nBandCount || 3==nBandCount )
		{
			panBandMap = new int[nBandCount];
			for( int j = 0; j < nBandCount; ++j )
				panBandMap[j] = j+1;//nBandCount-j;
		}
		//如果是7波段
		else if( 7==nBandCount )
		{
			nBandCount = 3;
			panBandMap = new int[nBandCount];
			panBandMap[0] = 3;
			panBandMap[1] = 4;
			panBandMap[2] = 2;
		}
		if( GetBitmapInfo( g_vecRI[i], panBandMap, showRange, pBitmap[nBitmap] ) )
		{
			//只打开可显示的前4个
			if( (pBitmap[nBitmap].isInters) && (++nBitmap==4) )
			{
				//扫描其余RasterInfo，关闭打开的pDataset
				for( int j = i+1; j < nVecRI; ++j )
				{
					if( g_vecRI[j].pDataset )
					{
						GDALClose( g_vecRI[j].pDataset );
						g_vecRI[j].pDataset = NULL;
					}
				}
				break;
			}
		}		
		delete[] panBandMap;
	}
	return;
}

GDAL获取指定地理坐标像元值、获取指定地理范围影像数据

标签：gdal   gis   数据   algorithm   

原文地址：http://blog.csdn.net/liyuan_669/article/details/45221759