GDAL:DEM渲染和等高线生成

//DEMRenderingontours.h

#include "gdal_priv.h"
#include"gdalwarper.h"
#include "ogr_spatialref.h"
#include "ogrsf_frmts.h"
#include "process.h"
#include<iostream>
using namespace std;
/**
* @brief 颜色插值模式
*/
typedef enum
{
	/*! 线性插值 */
	COLOR_SELECTION_INTERPOLATE = 0,
	/*! 最邻近插值 */
	COLOR_SELECTION_NEAREST_ENTRY = 1,
	/*! 严格查询插值，使用该值时，需要在颜色表文件中列举每个DEM高程值对应的RGB值，否则结果可能是黑色的 */
	COLOR_SELECTION_EXACT_ENTRY = 2
} ColorSelectionMode;
/*!成功执行*/
const int RE_SUCCESS = 0;
/*!文件不存在*/
const int RE_FILENOTEXIST = 1;
/*!文件格式不被支持*/
const int RE_FILENOTSUPPORT = 2;
/*!图像数据类型不正确*/
const int RE_FILETYPEERROR = 3;
/*!创建图像失败*/
const int RE_CREATEFAILED = 4;
/*!输入参数错误*/
const int RE_PARAMERROR = 5;
/*!其他错误*/
const int RE_FAILED = 6;
/*!图像不存在公共区域*/
const int RE_NOSAMEEXTENT = 7;
/*!用户取消操作*/
const int RE_USERCANCEL = 8;
/*!文件已经被使用*/
const int RE_FILEISUESED = 9;
/*!不支持的像素深度*/
const int RE_DEPTHNOTSUPPORT = 10;
/*!波段数量不符合要求*/
const int RE_BANDCOUNTERROR = 11;
/*!文件不存在投影*/
const int RE_NOPROJECTION = 12;
/*!投影不一致*/
const int RE_PROJECTIONDIFF = 13;
class DEMRenderingontours{
public:
	/**
	* @brief DEM Color Relief（DEM颜色渲染）
	* @param pszSrcFilename		输入文件路径
	* @param pszColorTextFile	输入颜色映射文件，如果不指定，自动使用USGS配色方案，GDAL_DATA\\usgs_dem_color_relief.clr
	* 颜色映射文件的格式是：“高程值 red green blue”,如：“1022 255 125 32”
	* @param pszDstFilename		输出文件路径
	* @param sColorMode			颜色插值模式，默认为COLOR_SELECTION_INTERPOLATE
	* @param bAddAlpha			是否添加Alpha通道，默认为不加
	* @param bComputeAtEdges	是否计算边界，默认为不计算
	* @param pszFormat			输出文件格式，详细参考GDAL支持数据类型
	* @param pProcess			进度条指针
	* @return 返回值，表示计算过程中出现的各种错误信息
	*/

	int DemColorRelief(const char* pszSrcFilename, const char* pszColorTextFile, const char* pszDstFilename,
		ColorSelectionMode sColorMode = COLOR_SELECTION_INTERPOLATE, bool bAddAlpha = false,
		bool bComputeAtEdges = false, const char *pszFormat = "GTiff", CProcessBase* pProcess = NULL);
	/**
	* @brief DEM生成等高线
	* @param pszSrcDEM			输入DEM文件路径
	* @param pszDstShp			输出等高线文件路径
	* @param iBandIndex			要处理的波段序号，默认为第一波段
	* @param dInterval			等高距，默认为10米
	* @param pszFormat			输出文件格式，详细参考GDAL支持数据类型
	* @param pProcess			进度条指针
	* @return  返回值，表示计算过程中出现的各种错误信息
	*/
	int DEM2Contour(const char* pszSrcDEM, const char* pszDstShp, int iBandIndex = 1, double dInterval = 10.0,
		const char * pszFormat = "ESRI Shapefile", CProcessBase* pProcess = NULL);
};

//DEMRenderingontours.cpp

#include "DEMRenderingontours.h"
typedef struct
{
	const char *name;
	float r;
	float g;
	float b;
} NamedColor;
static const NamedColor namedColors[] = {
	{ "white", 1.00, 1.00, 1.00 },
	{ "black", 0.00, 0.00, 0.00 },
	{ "red", 1.00, 0.00, 0.00 },
	{ "green", 0.00, 1.00, 0.00 },
	{ "blue", 0.00, 0.00, 1.00 },
	{ "yellow", 1.00, 1.00, 0.00 },
	{ "magenta", 1.00, 0.00, 1.00 },
	{ "cyan", 0.00, 1.00, 1.00 },
	{ "aqua", 0.00, 0.75, 0.75 },
	{ "grey", 0.75, 0.75, 0.75 },
	{ "gray", 0.75, 0.75, 0.75 },
	{ "orange", 1.00, 0.50, 0.00 },
	{ "brown", 0.75, 0.50, 0.25 },
	{ "purple", 0.50, 0.00, 1.00 },
	{ "violet", 0.50, 0.00, 1.00 },
	{ "indigo", 0.00, 0.50, 1.00 },
};
static
int GDALColorReliefFindNamedColor(const char *pszColorName, int *pnR, int *pnG, int *pnB)
{
	unsigned int i;

	*pnR = *pnG = *pnB = 0;
	for (i = 0; i < sizeof(namedColors) / sizeof(namedColors[0]); i++)
	{
		if (EQUAL(pszColorName, namedColors[i].name))
		{
			*pnR = (int)(255. * namedColors[i].r);
			*pnG = (int)(255. * namedColors[i].g);
			*pnB = (int)(255. * namedColors[i].b);
			return TRUE;
		}
	}
	return FALSE;
}
static double GDALColorReliefGetAbsoluteValFromPct(GDALRasterBandH hSrcBand,
	double dfPct)
{
	double dfMin, dfMax;
	int bSuccessMin, bSuccessMax;
	dfMin = GDALGetRasterMinimum(hSrcBand, &bSuccessMin);
	dfMax = GDALGetRasterMaximum(hSrcBand, &bSuccessMax);
	if (!bSuccessMin || !bSuccessMax)
	{
		double dfMean, dfStdDev;
		CPLDebug("DEMAlgortithm", "Computing source raster statistics...\n");
		GDALComputeRasterStatistics(hSrcBand, FALSE, &dfMin, &dfMax,
			&dfMean, &dfStdDev, NULL, NULL);
	}
	return dfMin + dfPct * (dfMax - dfMin);
}
typedef struct
{
	double dfVal;
	int nR;
	int nG;
	int nB;
	int nA;
} ColorAssociation;

static int GDALColorReliefSortColors(const void* pA, const void* pB)
{
	ColorAssociation* pC1 = (ColorAssociation*)pA;
	ColorAssociation* pC2 = (ColorAssociation*)pB;
	return (pC1->dfVal < pC2->dfVal) ? -1 :
		(pC1->dfVal == pC2->dfVal) ? 0 : 1;
}

static int GDALColorReliefGetRGBA(ColorAssociation* pasColorAssociation,
	int nColorAssociation,
	double dfVal,
	ColorSelectionMode eColorSelectionMode,
	int* pnR,
	int* pnG,
	int* pnB,
	int* pnA)
{
	int i;
	int lower = 0;
	int upper = nColorAssociation - 1;
	int mid;

	while (TRUE)
	{
		mid = (lower + upper) / 2;
		if (upper - lower <= 1)
		{
			if (dfVal < pasColorAssociation[lower].dfVal)
				i = lower;
			else if (dfVal < pasColorAssociation[upper].dfVal)
				i = upper;
			else
				i = upper + 1;
			break;
		}
		else if (pasColorAssociation[mid].dfVal >= dfVal)
		{
			upper = mid;
		}
		else
		{
			lower = mid;
		}
	}

	if (i == 0)
	{
		if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY &&
			pasColorAssociation[0].dfVal != dfVal)
		{
			*pnR = 0;
			*pnG = 0;
			*pnB = 0;
			*pnA = 0;
			return FALSE;
		}
		else
		{
			*pnR = pasColorAssociation[0].nR;
			*pnG = pasColorAssociation[0].nG;
			*pnB = pasColorAssociation[0].nB;
			*pnA = pasColorAssociation[0].nA;
			return TRUE;
		}
	}
	else if (i == nColorAssociation)
	{
		if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY &&
			pasColorAssociation[i - 1].dfVal != dfVal)
		{
			*pnR = 0;
			*pnG = 0;
			*pnB = 0;
			*pnA = 0;
			return FALSE;
		}
		else
		{
			*pnR = pasColorAssociation[i - 1].nR;
			*pnG = pasColorAssociation[i - 1].nG;
			*pnB = pasColorAssociation[i - 1].nB;
			*pnA = pasColorAssociation[i - 1].nA;
			return TRUE;
		}
	}
	else
	{
		if (eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY &&
			pasColorAssociation[i - 1].dfVal != dfVal)
		{
			*pnR = 0;
			*pnG = 0;
			*pnB = 0;
			*pnA = 0;
			return FALSE;
		}

		if (eColorSelectionMode == COLOR_SELECTION_NEAREST_ENTRY &&
			pasColorAssociation[i - 1].dfVal != dfVal)
		{
			int index;
			if (dfVal - pasColorAssociation[i - 1].dfVal <
				pasColorAssociation[i].dfVal - dfVal)
				index = i - 1;
			else
				index = i;

			*pnR = pasColorAssociation[index].nR;
			*pnG = pasColorAssociation[index].nG;
			*pnB = pasColorAssociation[index].nB;
			*pnA = pasColorAssociation[index].nA;
			return TRUE;
		}

		double dfRatio = (dfVal - pasColorAssociation[i - 1].dfVal) /
			(pasColorAssociation[i].dfVal - pasColorAssociation[i - 1].dfVal);
		*pnR = (int)(0.45 + pasColorAssociation[i - 1].nR + dfRatio *
			(pasColorAssociation[i].nR - pasColorAssociation[i - 1].nR));
		if (*pnR < 0) *pnR = 0;
		else if (*pnR > 255) *pnR = 255;
		*pnG = (int)(0.45 + pasColorAssociation[i - 1].nG + dfRatio *
			(pasColorAssociation[i].nG - pasColorAssociation[i - 1].nG));
		if (*pnG < 0) *pnG = 0;
		else if (*pnG > 255) *pnG = 255;
		*pnB = (int)(0.45 + pasColorAssociation[i - 1].nB + dfRatio *
			(pasColorAssociation[i].nB - pasColorAssociation[i - 1].nB));
		if (*pnB < 0) *pnB = 0;
		else if (*pnB > 255) *pnB = 255;
		*pnA = (int)(0.45 + pasColorAssociation[i - 1].nA + dfRatio *
			(pasColorAssociation[i].nA - pasColorAssociation[i - 1].nA));
		if (*pnA < 0) *pnA = 0;
		else if (*pnA > 255) *pnA = 255;

		return TRUE;
	}
}
static
ColorAssociation* GDALColorReliefParseColorFile(GDALRasterBandH hSrcBand,
const char* pszColorFilename,
int* pnColors)
{
	FILE* fpColorFile = VSIFOpenL(pszColorFilename, "rt");
	if (fpColorFile == NULL)
	{
		CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s", pszColorFilename);
		*pnColors = 0;
		return NULL;
	}

	ColorAssociation* pasColorAssociation = NULL;
	int nColorAssociation = 0;

	int bSrcHasNoData = FALSE;
	double dfSrcNoDataValue = GDALGetRasterNoDataValue(hSrcBand, &bSrcHasNoData);

	const char* pszLine;
	while ((pszLine = CPLReadLineL(fpColorFile)) != NULL)
	{
		char** papszFields = CSLTokenizeStringComplex(pszLine, " ,\t:", FALSE, FALSE);
		/* 跳过注释行 */
		int nTokens = CSLCount(papszFields);
		if (nTokens >= 2 &&
			papszFields[0][0] != '#' &&
			papszFields[0][0] != '/')
		{
			pasColorAssociation =
				(ColorAssociation*)CPLRealloc(pasColorAssociation,
				(nColorAssociation + 1) * sizeof(ColorAssociation));
			if (EQUAL(papszFields[0], "nv") && bSrcHasNoData)
				pasColorAssociation[nColorAssociation].dfVal = dfSrcNoDataValue;
			else if (strlen(papszFields[0]) > 1 && papszFields[0][strlen(papszFields[0]) - 1] == '%')
			{
				double dfPct = atof(papszFields[0]) / 100.;
				if (dfPct < 0.0 || dfPct > 1.0)
				{
					CPLError(CE_Failure, CPLE_AppDefined,
						"Wrong value for a percentage : %s", papszFields[0]);
					CSLDestroy(papszFields);
					VSIFCloseL(fpColorFile);
					CPLFree(pasColorAssociation);
					*pnColors = 0;
					return NULL;
				}
				pasColorAssociation[nColorAssociation].dfVal =
					GDALColorReliefGetAbsoluteValFromPct(hSrcBand, dfPct);
			}
			else
				pasColorAssociation[nColorAssociation].dfVal = atof(papszFields[0]);

			if (nTokens >= 4)
			{
				pasColorAssociation[nColorAssociation].nR = atoi(papszFields[1]);
				pasColorAssociation[nColorAssociation].nG = atoi(papszFields[2]);
				pasColorAssociation[nColorAssociation].nB = atoi(papszFields[3]);
				pasColorAssociation[nColorAssociation].nA =
					(CSLCount(papszFields) >= 5) ? atoi(papszFields[4]) : 255;
			}
			else
			{
				int nR, nG, nB;
				if (!GDALColorReliefFindNamedColor(papszFields[1], &nR, &nG, &nB))
				{
					CPLError(CE_Failure, CPLE_AppDefined,
						"Unknown color : %s", papszFields[1]);
					CSLDestroy(papszFields);
					VSIFCloseL(fpColorFile);
					CPLFree(pasColorAssociation);
					*pnColors = 0;
					return NULL;
				}
				pasColorAssociation[nColorAssociation].nR = nR;
				pasColorAssociation[nColorAssociation].nG = nG;
				pasColorAssociation[nColorAssociation].nB = nB;
				pasColorAssociation[nColorAssociation].nA =
					(CSLCount(papszFields) >= 3) ? atoi(papszFields[2]) : 255;
			}

			nColorAssociation++;
		}
		CSLDestroy(papszFields);
	}
	VSIFCloseL(fpColorFile);

	if (nColorAssociation == 0)
	{
		CPLError(CE_Failure, CPLE_AppDefined,
			"No color association found in %s", pszColorFilename);
		*pnColors = 0;
		return NULL;
	}

	qsort(pasColorAssociation, nColorAssociation,
		sizeof(ColorAssociation), GDALColorReliefSortColors);

	*pnColors = nColorAssociation;
	return pasColorAssociation;
}
static
GByte* GDALColorReliefPrecompute(GDALRasterBandH hSrcBand,
ColorAssociation* pasColorAssociation,
int nColorAssociation,
ColorSelectionMode eColorSelectionMode,
int* pnIndexOffset)
{
	GDALDataType eDT = GDALGetRasterDataType(hSrcBand);
	GByte* pabyPrecomputed = NULL;
	int nIndexOffset = (eDT == GDT_Int16) ? 32768 : 0;
	*pnIndexOffset = nIndexOffset;
	int nXSize = GDALGetRasterBandXSize(hSrcBand);
	int nYSize = GDALGetRasterBandXSize(hSrcBand);
	if (eDT == GDT_Byte ||
		((eDT == GDT_Int16 || eDT == GDT_UInt16) && nXSize * nYSize > 65536))
	{
		int iMax = (eDT == GDT_Byte) ? 256 : 65536;
		pabyPrecomputed = (GByte*)VSIMalloc(4 * iMax);
		if (pabyPrecomputed)
		{
			int i;
			for (i = 0; i < iMax; i++)
			{
				int nR, nG, nB, nA;
				GDALColorReliefGetRGBA(pasColorAssociation,
					nColorAssociation,
					i - nIndexOffset,
					eColorSelectionMode,
					&nR, &nG, &nB, &nA);
				pabyPrecomputed[4 * i] = (GByte)nR;
				pabyPrecomputed[4 * i + 1] = (GByte)nG;
				pabyPrecomputed[4 * i + 2] = (GByte)nB;
				pabyPrecomputed[4 * i + 3] = (GByte)nA;
			}
		}
	}
	return pabyPrecomputed;
}
CPLErr GDALColorRelief(GDALRasterBandH hSrcBand,
	GDALRasterBandH hDstBand1,
	GDALRasterBandH hDstBand2,
	GDALRasterBandH hDstBand3,
	GDALRasterBandH hDstBand4,
	const char* pszColorFilename,
	ColorSelectionMode eColorSelectionMode,
	GDALProgressFunc pfnProgress,
	void * pProgressData)
{
	CPLErr eErr;

	if (hSrcBand == NULL || hDstBand1 == NULL || hDstBand2 == NULL ||
		hDstBand3 == NULL)
		return CE_Failure;

	int nColorAssociation = 0;
	ColorAssociation* pasColorAssociation =
		GDALColorReliefParseColorFile(hSrcBand, pszColorFilename,
		&nColorAssociation);
	if (pasColorAssociation == NULL)
		return CE_Failure;

	int nXSize = GDALGetRasterBandXSize(hSrcBand);
	int nYSize = GDALGetRasterBandYSize(hSrcBand);

	if (pfnProgress == NULL)
		pfnProgress = GDALDummyProgress;

	int nR = 0, nG = 0, nB = 0, nA = 0;

	int nIndexOffset = 0;
	GByte* pabyPrecomputed =
		GDALColorReliefPrecompute(hSrcBand,
		pasColorAssociation,
		nColorAssociation,
		eColorSelectionMode,
		&nIndexOffset);

	float* pafSourceBuf = NULL;
	int* panSourceBuf = NULL;
	if (pabyPrecomputed)
		panSourceBuf = (int *)CPLMalloc(sizeof(int)*nXSize);
	else
		pafSourceBuf = (float *)CPLMalloc(sizeof(float)*nXSize);
	GByte* pabyDestBuf1 = (GByte*)CPLMalloc(4 * nXSize);
	GByte* pabyDestBuf2 = pabyDestBuf1 + nXSize;
	GByte* pabyDestBuf3 = pabyDestBuf2 + nXSize;
	GByte* pabyDestBuf4 = pabyDestBuf3 + nXSize;
	int i, j;

	if (!pfnProgress(0.0, NULL, pProgressData))
	{
		CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
		eErr = CE_Failure;
		goto end;
	}

	for (i = 0; i < nYSize; i++)
	{
		eErr = GDALRasterIO(hSrcBand,
			GF_Read,
			0, i,
			nXSize, 1,
			(panSourceBuf) ? (void*)panSourceBuf : (void*)pafSourceBuf,
			nXSize, 1,
			(panSourceBuf) ? GDT_Int32 : GDT_Float32,
			0, 0);
		if (eErr != CE_None)
			goto end;

		if (panSourceBuf)
		{
			for (j = 0; j < nXSize; j++)
			{
				int nIndex = panSourceBuf[j] + nIndexOffset;
				pabyDestBuf1[j] = pabyPrecomputed[4 * nIndex];
				pabyDestBuf2[j] = pabyPrecomputed[4 * nIndex + 1];
				pabyDestBuf3[j] = pabyPrecomputed[4 * nIndex + 2];
				pabyDestBuf4[j] = pabyPrecomputed[4 * nIndex + 3];
			}
		}
		else
		{
			for (j = 0; j < nXSize; j++)
			{
				GDALColorReliefGetRGBA(pasColorAssociation,
					nColorAssociation,
					pafSourceBuf[j],
					eColorSelectionMode,
					&nR,
					&nG,
					&nB,
					&nA);
				pabyDestBuf1[j] = (GByte)nR;
				pabyDestBuf2[j] = (GByte)nG;
				pabyDestBuf3[j] = (GByte)nB;
				pabyDestBuf4[j] = (GByte)nA;
			}
		}

		eErr = GDALRasterIO(hDstBand1,
			GF_Write,
			0, i, nXSize,
			1, pabyDestBuf1, nXSize, 1, GDT_Byte, 0, 0);
		if (eErr != CE_None)
			goto end;

		eErr = GDALRasterIO(hDstBand2,
			GF_Write,
			0, i, nXSize,
			1, pabyDestBuf2, nXSize, 1, GDT_Byte, 0, 0);
		if (eErr != CE_None)
			goto end;

		eErr = GDALRasterIO(hDstBand3,
			GF_Write,
			0, i, nXSize,
			1, pabyDestBuf3, nXSize, 1, GDT_Byte, 0, 0);
		if (eErr != CE_None)
			goto end;

		if (hDstBand4)
		{
			eErr = GDALRasterIO(hDstBand4,
				GF_Write,
				0, i, nXSize,
				1, pabyDestBuf4, nXSize, 1, GDT_Byte, 0, 0);
			if (eErr != CE_None)
				goto end;
		}

		if (!pfnProgress(1.0 * (i + 1) / nYSize, NULL, pProgressData))
		{
			CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
			eErr = CE_Failure;
			goto end;
		}
	}
	pfnProgress(1.0, NULL, pProgressData);
	eErr = CE_None;

end:
	VSIFree(pabyPrecomputed);
	CPLFree(pafSourceBuf);
	CPLFree(panSourceBuf);
	CPLFree(pabyDestBuf1);
	CPLFree(pasColorAssociation);

	return eErr;
}

int DEMRenderingontours::DemColorRelief(const char* pszSrcFilename, const char* pszColorTextFile, const char* pszDstFilename,
	ColorSelectionMode sColorMode, bool bAddAlpha, bool bComputeAtEdges,
	const char *pszFormat, CProcessBase* pProcess)
{
	if (pProcess != NULL)
	{
		pProcess->ReSetProcess();
		pProcess->SetMessage("正在计算颜色渲染...");
	}

	int nBand = 1;
	double  adfGeoTransform[6];

	char **papszCreateOptions = NULL;

	GDALDatasetH hSrcDataset = NULL;
	GDALDatasetH hDstDataset = NULL;
	GDALRasterBandH hSrcBand = NULL;
	GDALRasterBandH hDstBand = NULL;
	GDALDriverH hDriver = NULL;

	GDALProgressFunc pfnProgress = ALGTermProgress;

	int nXSize = 0;
	int nYSize = 0;

	if (pszSrcFilename == NULL)
	{
		if (pProcess != NULL)
			pProcess->SetMessage("源文件路径为空！");

		return RE_FILENOTEXIST;
	}

	if (pszDstFilename == NULL)
	{
		if (pProcess != NULL)
			pProcess->SetMessage("输出文件路径为空！");

		return RE_FILENOTEXIST;
	}

	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
	//打开图像
	hSrcDataset = GDALOpen(pszSrcFilename, GA_ReadOnly);
	if (hSrcDataset == NULL)
	{
		if (pProcess != NULL)
			pProcess->SetMessage("输入文件不能打开！");

		return RE_FILENOTEXIST;
	}

	nXSize = GDALGetRasterXSize(hSrcDataset);
	nYSize = GDALGetRasterYSize(hSrcDataset);

	if (nBand <= 0 || nBand > GDALGetRasterCount(hSrcDataset))
	{
		if (pProcess != NULL)
			pProcess->SetMessage("指定计算的波段不正确！");

		return RE_PARAMERROR;
	}

	hSrcBand = GDALGetRasterBand(hSrcDataset, nBand);
	GDALGetGeoTransform(hSrcDataset, adfGeoTransform);

	hDriver = GDALGetDriverByName(pszFormat);
	if (hDriver == NULL)
	{
		if (pProcess != NULL)
			pProcess->SetMessage("不能创建制定类型的文件，请检查该文件类型GDAL是否支持创建！");

		GDALClose(hSrcDataset);

		return RE_FILENOTSUPPORT;
	}

	//设置输出文件数据类型
	GDALDataType eDstDataType = GDT_Byte;

	int nDstBands = (bAddAlpha) ? 4 : 3;

	hDstDataset = GDALCreate(hDriver,
		pszDstFilename,
		nXSize,
		nYSize,
		nDstBands,
		eDstDataType,
		papszCreateOptions);

	if (hDstDataset == NULL)
	{
		if (pProcess != NULL)
			pProcess->SetMessage("创建输出文件失败，请检查文件路径是否正确！");

		GDALClose(hSrcDataset);

		return RE_CREATEFAILED;	//创建图像失败
	}

	hDstBand = GDALGetRasterBand(hDstDataset, 1);

	GDALSetGeoTransform(hDstDataset, adfGeoTransform);
	GDALSetProjection(hDstDataset, GDALGetProjectionRef(hSrcDataset));

	if (GDALColorRelief(hSrcBand,
		GDALGetRasterBand(hDstDataset, 1),
		GDALGetRasterBand(hDstDataset, 2),
		GDALGetRasterBand(hDstDataset, 3),
		(bAddAlpha) ? GDALGetRasterBand(hDstDataset, 4) : NULL,
		pszColorTextFile,
		sColorMode,
		pfnProgress, pProcess) != CE_None)
	{
		if (pProcess != NULL)
		{
			if (!pProcess->m_bIsContinue)
			{
				GDALClose(hSrcDataset);
				GDALClose(hDstDataset);

				GDALDeleteDataset(hDstDataset, pszDstFilename);	//删除结果文件

				CSLDestroy(papszCreateOptions);

				return RE_USERCANCEL;
			}
		}

		GDALClose(hSrcDataset);
		GDALClose(hDstDataset);

		GDALDeleteDataset(hDstDataset, pszDstFilename);	//删除结果文件
		CSLDestroy(papszCreateOptions);
		return RE_FAILED;
	}

	GDALClose(hSrcDataset);
	GDALClose(hDstDataset);

	CSLDestroy(papszCreateOptions);

	return RE_SUCCESS;
}
int DEMRenderingontours::DEM2Contour(const char* pszSrcDEM, const char* pszDstShp, int iBandIndex, double dInterval,
	const char * pszFormat , CProcessBase* pProcess){
	// 下面的代码假设所有的输入参数都正确，不进行错误判断
	if (pProcess != NULL)
	{
		pProcess->ReSetProcess();
		pProcess->SetMessage("开始DEM生成等高线...");
	}

	GDALAllRegister();
	OGRRegisterAll();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
	GDALDatasetH hSrcDS = GDALOpen(pszSrcDEM, GA_ReadOnly);
	GDALRasterBandH hBand = GDALGetRasterBand(hSrcDS, iBandIndex);

	// 获取原始DEM数据中的NODATA值
	int bNoDataSet = FALSE, bIgnoreNoData = FALSE;
	double dfNoData = 0.0;
	if (!bNoDataSet && !bIgnoreNoData)
		dfNoData = GDALGetRasterNoDataValue(hBand, &bNoDataSet);

	// 从输入的DEM数据中获取空间参考信息
	OGRSpatialReferenceH hSRS = NULL;
	const char *pszWKT = GDALGetProjectionRef(hSrcDS);

	if (pszWKT != NULL && strlen(pszWKT) != 0)
		hSRS = OSRNewSpatialReference(pszWKT);

	// 创建输出等高线矢量文件
	OGRSFDriverH hDriver = OGRGetDriverByName(pszFormat);
	OGRDataSourceH hDS = OGR_Dr_CreateDataSource(hDriver, pszDstShp, NULL);
	OGRLayerH hLayer = OGR_DS_CreateLayer(hDS, "Contour", hSRS, wkbLineString, NULL);

	// 创建保存等高线ID的属性字段
	OGRFieldDefnH hFld;
	hFld = OGR_Fld_Create("ID", OFTInteger);
	OGR_Fld_SetWidth(hFld, 8);
	OGR_L_CreateField(hLayer, hFld, FALSE);
	OGR_Fld_Destroy(hFld);

	// 创建保存等高线高程值的属性字段
	hFld = OGR_Fld_Create("Elevation", OFTReal);
	OGR_Fld_SetWidth(hFld, 12);
	OGR_Fld_SetPrecision(hFld, 3);
	OGR_L_CreateField(hLayer, hFld, FALSE);
	OGR_Fld_Destroy(hFld);
	int nElevField = 1;

	// 调用GDAL库中的函数生成等高线
	CPLErr eErr = GDALContourGenerate(hBand, dInterval, 0.0,
		0, NULL,
		bNoDataSet, dfNoData,
		hLayer, 0, nElevField,
		GDALTermProgress, pProcess);

	OGR_DS_Destroy(hDS);
	GDALClose(hSrcDS);

	GDALDestroyDriverManager();
	OGRCleanupAll();

	return RE_SUCCESS;

}

//调用

#include "DEMRenderingontours.h"
int main(){
	/********************************************************************************************************/
	//------------------------------------------DEM渲染-----------------------------------------------------//
	/********************************************************************************************************/
	const char* pszSrcFile3 = "C:/Users/Administrator/Documents/Tencent Files/328302992/FileRecv/GDAL书籍代码及数据/GDAL书籍数据/7-Data/dem.tif";
	const char* pszDstFile3 = "C:/Users/Administrator/Documents/Tencent Files/328302992/FileRecv/GDAL书籍代码及数据/GDAL书籍数据/7-Data/dem_渲染后1.tif";
	const char* pszcolorFile3 = "C:/Users/Administrator/Documents/Tencent Files/328302992/FileRecv/GDAL书籍代码及数据/GDAL书籍数据/7-Data/color1.txt";
	CConsoleProcess *pProdem = new CConsoleProcess;
	DEMRenderingontours ImageCorrectbyDEM;
	ImageCorrectbyDEM.DemColorRelief(pszSrcFile3, pszcolorFile3, pszDstFile3, COLOR_SELECTION_INTERPOLATE, false, false, "GTiff", pProdem);
	delete pProdem;
	/********************************************************************************************************/
	//------------------------------------------DEM生成等高线-----------------------------------------------//
	/********************************************************************************************************/
	const char* pszSrcFile4 = "C:/Users/Administrator/Documents/Tencent Files/328302992/FileRecv/GDAL书籍代码及数据/GDAL书籍数据/7-Data/dem.tif";
	const char* pszDstFile4 = "C:/Users/Administrator/Documents/Tencent Files/328302992/FileRecv/GDAL书籍代码及数据/GDAL书籍数据/7-Data/dem_counter.shp";
	CConsoleProcess *pProdemcounter = new CConsoleProcess;
	DEMRenderingontours ImageCorrectbyDEMcounter;
	ImageCorrectbyDEMcounter.DEM2Contour(pszSrcFile4, pszDstFile4, 1, 10.0, "ESRI Shapefile", pProdemcounter);
	delete pProdemcounter;
	return 0;
}

颜色映射表

两个映射表生成的渲染效果为：

等高线生成：