使用GDAL库实现对dem的修改

一、代码

#include "iostream"
#include "string"
#include "gdal_priv.h"
#include "cpl_conv.h"
#include "gdalwarper.h"
#include "stdlib.h"


typedef unsigned short UINT;

int main()
{
	//Opening the File
	GDALAllRegister();

	std::string pszFilename("srtm30plus_stripped.tif");

	//创建dem实体
	GDALDataset *poDataset;
	poDataset = (GDALDataset*)GDALOpen(pszFilename.c_str(), GA_ReadOnly);
	if (poDataset == NULL)
	{
		std::cout << "OPEN GEOTIFF FAIL!" << std::endl;
		return -1;
	}

	//Getting Dataset Information
	double adfGeoTransform[6];
	float width = poDataset->GetRasterXSize();
	float height = poDataset->GetRasterYSize();
	float nBands = poDataset->GetRasterCount();

	std::cout << "Driver: " << poDataset->GetDriver()->GetDescription() << " , "
		<< poDataset->GetDriver()->GetMetadataItem(GDAL_DMD_LONGNAME) << std::endl;

	std::cout << "Size is " << width << " x " <<
		height << " x " << nBands << std::endl;

	if (poDataset->GetProjectionRef() != NULL)
	{
		std::cout << "Projection is" << poDataset->GetProjectionRef() << std::endl;
	}

	if (poDataset->GetGeoTransform(adfGeoTransform) == CE_None)
	{
		std::cout << "Origin = ( " << adfGeoTransform[0] << " , " << adfGeoTransform[3]
			<< " )" << std::endl;
		std::cout << "Pixel Size = ( " << adfGeoTransform[1] << " , " << adfGeoTransform[5]
			<< " )" << std::endl;
	}


	//Fetching a Raster Band
	GDALRasterBand *poBand;
	int nBlockXSize, nBlockYSize;
	int bGotMin, bGotMax;
	double adfMinMax[2];
	poBand = poDataset->GetRasterBand(1);
	poBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
	std::cout << "Block = " << nBlockXSize << " x " << nBlockYSize << " , " <<
		"Type = " << GDALGetDataTypeName(poBand->GetRasterDataType()) << " , " <<
		"ColorInterp = " << GDALGetColorInterpretationName(poBand->GetColorInterpretation()) << std::endl;
	adfMinMax[0] = poBand->GetMinimum(&bGotMin);
	adfMinMax[1] = poBand->GetMaximum(&bGotMax);
	
	if (!(bGotMin&&bGotMax))
	{
		GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);
	}
	std::cout << "Min = " << adfMinMax[0] << " , " << " Max = " << adfMinMax[1] << std::endl;
	if (poBand->GetOverviewCount() > 0)
	{
		std::cout << "Band has " << poBand->GetOverviewCount() << " overviews." << std::endl;
	}
	if (poBand->GetColorTable() != NULL)
	{
		std::cout << "Band has a color table with " << poBand->GetColorTable()->GetColorEntryCount() <<
			" entries." << std::endl;
	}


	//Reading Raster Data And Create File
	float **pafScanline=new float *[nBands];
	if (pafScanline == NULL)
	{
		std::cout << "UINT **pafScanline = new UINT *[nBands]; failed!" << std::endl;
		return -1;
	}
	float **pafScanlineOutput = new float *[nBands];
	if (pafScanlineOutput == NULL)
	{
		std::cout << "UINT **pafScanlineOutput = new UINT *[nBands]; failed!" << std::endl;
		return -1;
	}

	int nXSize = poBand->GetXSize();
	int nYSize = poBand->GetYSize();
	//pafScanline = (float *)CPLMalloc(sizeof(float)*nXSize*nYSize);
	//poBand->RasterIO(GF_Read, 0, 0, nXSize, nYSize,
	//	pafScanline, nXSize, nYSize, /*poBand->GetRasterDataType()*/GDT_Float32, 0, 0);

	const char *pszFormat = "GTiff";
	GDALDriver *poDriver;
	poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);

	GDALDataset *poDstDs;
	const char *pszDstFilename = "srtm_30m_output.tif";
	char **papxzOptions = NULL;

	poDstDs = poDriver->Create(pszDstFilename, nXSize, nYSize, poDataset->GetRasterCount(),
		GDT_Float32, papxzOptions);

	//一次处理xHei=400行数据  分块处理
	int xHei = 400;
	int timePerxH = nYSize / xHei + 1;

	//防止无符号数越界
	int tmp = 0;

	//开始分块处理
	//i表示分块处理需处理次数计数
	for (int i = 0; i < timePerxH - 1; i++)
	{
		for (int j = 0; j < nBands; j++)
		{
			pafScanline[j] = new float[nXSize*xHei];
			if (pafScanline[j] == NULL)
			{
				std::cout << "pafScanline[j] = new UINT[Width*xHei]; failed!" << std::endl;
				return -1;
			}
			pafScanlineOutput[j] = new float[nXSize*xHei];
			if (pafScanlineOutput[j] == NULL)
			{
				std::cout << "pafScanlineOutput[j] = new UINT[Width*xHei]; failed!" << std::endl;
				return -1;
			}
		}
		poBand->RasterIO(GF_Read, 0, i*xHei, nXSize, xHei, pafScanline[0],
			nXSize, xHei, GDT_Float32, 0, 0);

		//pafScanline像素操作
		for (int i = 0; i < nBands; i++)
		{
			for (int j = 0; j < xHei; j++)
			{
				for (int k = 0; k < nXSize; k++)
				{
					if (pafScanline[i][j*nXSize + k] < 0)
					{
						//std::cout << "修改前" << pafScanline[i][j*nXSize + k] << std::endl;
						pafScanline[i][j*nXSize + k] = 0;
						//std::cout << "修改后" << pafScanline[i][j*nXSize + k] << std::endl;
					}

					pafScanlineOutput[i][j*nXSize + k] = pafScanline[i][j*nXSize + k];
				}
			}
		}
		poDstDs->GetRasterBand(1)->RasterIO(GF_Write, 0, xHei*i, nXSize, xHei,
			pafScanlineOutput[0], nXSize, xHei, GDT_Float32, 0, 0);

		for (int j = 0; j < nBands; j++)
		{
			delete []pafScanline[j];
			pafScanline[j] = NULL;

			delete []pafScanlineOutput[j];
			pafScanlineOutput[j] = NULL;
		}

	}

	//剩余行数处理
	int foreH = (timePerxH - 1)*xHei;
	int leftH = nYSize - foreH;
	for (int j = 0; j < nBands; j++)
	{
		pafScanline[j] = new float[nXSize*leftH];
		if (pafScanline[j] == NULL)
		{
			std::cout << "pafScanline[j] = new UINT[Width*leftH]; failed!" << std::endl;
			return -1;
		}
		pafScanlineOutput[j] = new float[nXSize*leftH];
		if (pafScanlineOutput[j] == NULL)
		{
			std::cout << "pafScanlineOutput[j] = new UINT[Width*leftH]; failed!" << std::endl;
			return -1;
		}
	}
	poBand->RasterIO(GF_Read, 0, foreH, nXSize, leftH, pafScanline[0], nXSize, 
		leftH, GDT_Float32, 0, 0);
	//pafScanline像素操作
	for (int i = 0; i < nBands; i++)
	{
		for (int j = 0; j < xHei; j++)
		{
			for (int k = 0; k < nXSize; k++)
			{
				if (pafScanline[i][j*nXSize + k] < 0)
				{
					//std::cout << "修改前" << pafScanline[i][j*nXSize + k] << std::endl;
					pafScanline[i][j*nXSize + k] = 0;
					//std::cout << "修改后" << pafScanline[i][j*nXSize + k] << std::endl;
				}

				pafScanlineOutput[i][j*nXSize + k] = pafScanline[i][j*nXSize + k];
				//std::cout << "pafScanlineOutput的数据 ：" << pafScanlineOutput[i][j*nXSize + k] << std::endl;
			}
		}
	}
	poDstDs->GetRasterBand(1)->RasterIO(GF_Write, 0, foreH, nXSize, leftH, 
		pafScanlineOutput[0], nXSize, leftH, GDT_Float32, 0, 0);
	for (int j = 0; j < nBands; j++)
	{
		delete[]pafScanline[j];
		pafScanline[j] = NULL;

		delete[]pafScanlineOutput[j];
		pafScanlineOutput[j] = NULL;
	}

	delete pafScanline;
	pafScanline = NULL;
	delete pafScanlineOutput;
	pafScanlineOutput = NULL;

	poDstDs->SetGeoTransform(adfGeoTransform);
	poDstDs->SetProjection(poDataset->GetProjectionRef());



	GDALClose(poDataset);

	system("pause");
	return 0;
}

二、心得

要注意保存像素数据的数组需要和图像文件本身存储的数据类型相匹配，否则出来的数据是错误的。一般来说dem文件都十分巨大，在读入时需要注意分段。

参考文献

[原][osg][gdal]两种方式修改tiff高程
GDAL documentation