使用GDAL實作DEM的地貌暈渲圖(三)

繼續(一)(二)的内容，生成DEM彩色的地貌暈渲圖。

1. 原理
- 1) ArcMap生成彩色暈渲圖
- 2) 彩色色帶指派
- 3) 顔色疊加
2. 實作
3. 結語
4. 參考

之前在《使用GDAL實作DEM的地貌暈渲圖(一)》和《使用GDAL實作DEM的地貌暈渲圖(二)》這兩篇文章中詳細介紹了DEM生成地貌暈渲圖的原理與實作。不過之前生成的都是暈渲強度值對應的灰階圖，而實際的應用過程中都會将DEM暈渲成彩色圖。

可以通過ArcMap的做法來參考如何生成彩色暈渲圖（參考[1]），在ArcMap中生成彩色暈渲圖的步驟如下：

通過山體陰影工具生成灰階暈渲圖，這一點與前面文章介紹的相一緻。
然後在原DEM圖的顯示中，選擇最大最小拉伸顯示，然後選擇一個合适的彩色色帶指派。
最後，将步驟一的灰階暈渲圖設定一定的透明度，疊加到步驟二的彩色圖上，就生成了最終具有立體感的彩色暈渲圖。

ArcMap生成的彩色暈渲圖：

不難發現，生成彩色暈渲圖的關鍵是第二步：要選取合适的色帶，讓色帶根據對應的高程指派。查閱了不少的資料，這個色帶應該沒有固定合适通用的模闆，是需要自己根據具體的需要調整的。比如，海平面可以指派成藍色；高山山頂指派成白色；戈壁沙漠指派成黃色；草原森林指派成綠色，這些地貌特征都具有一定的高程相關性，可以根據一定的絕對高程區間指派。

我這裡采取的做法還是跟ArcMap一緻，選取漸變色帶來指派，将漸變色帶限制到DEM的最小最大高程。考慮到地貌的多變性，我這裡生成了藍-綠-黃-紅-紫的多段的漸變色帶。這樣DEM的暈渲效果就是越低越藍，越高越紫。

一般為了保證過渡效果會選擇漸變色帶，漸變色帶的生成也比較簡單，選擇頭尾兩個的顔色的RGB值和一定的漸變範圍，分别讓RGB的值勻速變換就行了。彩色色帶的生成算法如下（可參考第二部分的具體實作來了解）：

//顔色查找表
vector<F_RGB> tableRGB(256);		

//生成漸變色
void Gradient(F_RGB &start, F_RGB &end, vector<F_RGB> &RGBList)
{
	float dr = (end.R - start.R) / RGBList.size();
	float dg = (end.G - start.G) / RGBList.size();
	float db = (end.B - start.B) / RGBList.size();
	for (size_t i = 0; i < RGBList.size(); i++)
	{
		RGBList[i].R = start.R + dr * i;
		RGBList[i].G = start.G + dg * i;
		RGBList[i].B = start.B + db * i;
	}
}

//初始化顔色查找表
void InitColorTable()
{
	F_RGB blue(17, 60, 235);//藍色	
	F_RGB green(17, 235, 86);//綠色
	vector<F_RGB> RGBList(60);
	Gradient(blue, green, RGBList);
	for (int i = 0; i < 60; i++)
	{
		tableRGB[i] = RGBList[i];
	}
	
	F_RGB yellow(235, 173, 17);//黃色	
	RGBList.clear();
	RGBList.resize(60);
	Gradient(green, yellow, RGBList);
	for (int i = 0; i < 60; i++)
	{
		tableRGB[i+60] = RGBList[i];
	}
	
	F_RGB red(235, 60, 17);//紅色
	RGBList.clear();
	RGBList.resize(60);
	Gradient(yellow, red, RGBList);
	for (int i = 0; i < 60; i++)
	{
		tableRGB[i + 120] = RGBList[i];
	}
		
	F_RGB white(235, 17, 235);//紫色
	RGBList.clear();
	RGBList.resize(76);
	Gradient(red, white, RGBList);
	for (int i = 0; i < 76; i++)
	{
		tableRGB[i + 180] = RGBList[i];
	}
}

第一步和第二步分别生成了暈渲強度圖和高程彩色色帶圖，第三步就是将兩者的顔色疊加，生成最終的效果圖。其實顔色疊加的原理特别簡單，對于暈渲強度圖的像素值A，令其不透明度為α；對應的高程彩色色帶圖的像素值B，那麼最後疊加的像素值 C=αA+(1-α)B。可以這麼了解:光線到達A，由于透光性隻産生了αA的效果，還有(1-α)強度的光線射到B，又産生了(1-α)B的效果，兩者疊加就是αA+(1-α)B。

繼續改造之前的代碼，最終的實作過程如下：

#include <iostream>
#include <algorithm>
#include <gdal_priv.h>
#include <osg/Vec3d>
#include <fstream>

using namespace std;
using namespace osg;

//RGB顔色
struct F_RGB
{
	float R;
	float G;
	float B;

	F_RGB()
	{
	}

	F_RGB(float r, float g, float b)
	{
		R = r;
		G = g;
		B = b;
	}

	F_RGB(const F_RGB& rgb) 
	{
		R = rgb.R;
		G = rgb.G;
		B = rgb.B;
	}

	F_RGB& operator= (const F_RGB& rgb)
	{	
		R = rgb.R;
		G = rgb.G;
		B = rgb.B;
		return *this;
	}

};

//顔色查找表
vector<F_RGB> tableRGB(256);		

//生成漸變色
void Gradient(F_RGB &start, F_RGB &end, vector<F_RGB> &RGBList)
{
	float dr = (end.R - start.R) / RGBList.size();
	float dg = (end.G - start.G) / RGBList.size();
	float db = (end.B - start.B) / RGBList.size();
	for (size_t i = 0; i < RGBList.size(); i++)
	{
		RGBList[i].R = start.R + dr * i;
		RGBList[i].G = start.G + dg * i;
		RGBList[i].B = start.B + db * i;
	}
}

//初始化顔色查找表
void InitColorTable()
{
	F_RGB blue(17, 60, 235);//藍色	
	F_RGB green(17, 235, 86);//綠色
	vector<F_RGB> RGBList(60);
	Gradient(blue, green, RGBList);
	for (int i = 0; i < 60; i++)
	{
		tableRGB[i] = RGBList[i];
	}
	
	F_RGB yellow(235, 173, 17);//黃色	
	RGBList.clear();
	RGBList.resize(60);
	Gradient(green, yellow, RGBList);
	for (int i = 0; i < 60; i++)
	{
		tableRGB[i+60] = RGBList[i];
	}
	
	F_RGB red(235, 60, 17);//紅色
	RGBList.clear();
	RGBList.resize(60);
	Gradient(yellow, red, RGBList);
	for (int i = 0; i < 60; i++)
	{
		tableRGB[i + 120] = RGBList[i];
	}
		
	F_RGB white(235, 17, 235);//紫色
	RGBList.clear();
	RGBList.resize(76);
	Gradient(red, white, RGBList);
	for (int i = 0; i < 76; i++)
	{
		tableRGB[i + 180] = RGBList[i];
	}
}

//根據高程選顔色
inline int GetColorIndex(float z, float min_z, float max_z)
{
	int temp = floor((z - min_z) * 255 / (max_z - min_z) + 0.6);
	return temp;
}

// a b c
// d e f
// g h i
double CalHillshade(float *tmpBuf, double Zenith_rad, double Azimuth_rad, double dx, double dy, double z_factor)
{
	double dzdx = ((tmpBuf[2] + 2 * tmpBuf[5] + tmpBuf[8]) - (tmpBuf[0] + 2 * tmpBuf[3] + tmpBuf[6])) / (8 * dx);
	double dzdy = ((tmpBuf[6] + 2 * tmpBuf[7] + tmpBuf[8]) - (tmpBuf[0] + 2 * tmpBuf[1] + tmpBuf[2])) / (8 * dy);

	double Slope_rad = atan(z_factor * sqrt(dzdx*dzdx + dzdy*dzdy));
	double Aspect_rad = 0;
	if (abs(dzdx) > 1e-9)
	{
		Aspect_rad = atan2(dzdy, -dzdx);
		if (Aspect_rad < 0)
		{
			Aspect_rad = 2 * PI + Aspect_rad;
		}
	}
	else
	{
		if (dzdy > 0)
		{
			Aspect_rad = PI / 2;
		}
		else if (dzdy < 0)
		{
			Aspect_rad = 2 * PI - PI / 2;
		}
		else
		{
			Aspect_rad = Aspect_rad;
		}
	}

	double Hillshade = 255.0 * ((cos(Zenith_rad) * cos(Slope_rad)) + (sin(Zenith_rad) * sin(Slope_rad) * cos(Azimuth_rad - Aspect_rad)));
	return Hillshade;
}


int main()
{
	GDALAllRegister();          //GDAL所有操作都需要先注冊格式
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");  //支援中文路徑
	CPLSetConfigOption("GDAL_DATA", "D:/Work/GDALBuild/gdal_build_result/data");  //支援中文路徑

	//const char* demPath = "D:/CloudSpace/我的技術文章/素材/DEM的渲染/dst.tif";
	const char* demPath = "D:/2.tif";
		
	GDALDataset* img = (GDALDataset *)GDALOpen(demPath, GA_ReadOnly);
	if (!img)
	{
		cout << "Can't Open Image!" << endl;
		return 1;
	}

	int imgWidth = img->GetRasterXSize();   //圖像寬度
	int imgHeight = img->GetRasterYSize();  //圖像高度
	int bandNum = img->GetRasterCount();    //波段數
	int depth = GDALGetDataTypeSize(img->GetRasterBand(1)->GetRasterDataType()) / 8;    //圖像深度

	GDALDriver *pDriver = GetGDALDriverManager()->GetDriverByName("GTIFF"); //圖像驅動
	char** ppszOptions = NULL;
	ppszOptions = CSLSetNameValue(ppszOptions, "BIGTIFF", "IF_NEEDED"); //配置圖像資訊
	const char* dstPath = "D:\\dst.tif";
	int bufWidth = imgWidth;
	int bufHeight = imgHeight;
	int dstBand = 3;
	int dstDepth = 1;
	GDALDataset* dst = pDriver->Create(dstPath, bufWidth, bufHeight, dstBand, GDT_Byte, ppszOptions);
	if (!dst)
	{
		printf("Can't Write Image!");
		return false;
	}

	dst->SetProjection(img->GetProjectionRef());
	double padfTransform[6] = { 0 };
	if (CE_None == img->GetGeoTransform(padfTransform))
	{
		dst->SetGeoTransform(padfTransform);
	}

	//申請buf
        depth = 4;
	size_t imgBufNum = (size_t)bufWidth * bufHeight * bandNum;
	float *imgBuf = new float[imgBufNum];
	//讀取
	img->RasterIO(GF_Read, 0, 0, bufWidth, bufHeight, imgBuf, bufWidth, bufHeight,
		GDT_Float32, bandNum, nullptr, bandNum*depth, bufWidth*bandNum*depth, depth);

	if (bandNum != 1)
	{
		return 1;
	}

	//
	double startX = padfTransform[0];			//左上角點坐标X
	double dx = padfTransform[1];			//X方向的分辨率
	double startY = padfTransform[3]; 			//左上角點坐标Y
	double dy = padfTransform[5];			//Y方向的分辨率

	//
	double minZ = DBL_MAX;
	double maxZ = -DBL_MAX;
	double noValue = img->GetRasterBand(1)->GetNoDataValue();
	
	//
	InitColorTable();
	for (int yi = 0; yi < bufHeight; yi++)
	{
		for (int xi = 0; xi < bufWidth; xi++)
		{
			size_t m = (size_t)bufWidth * yi + xi;
			double x = startX + xi * dx;
			double y = startY + yi * dy;
			double z = imgBuf[m];
		
			if (abs(z - noValue) < 0.01 || z<-11034 || z>8844.43)
			{
				continue;
			}

			minZ = (std::min)(minZ, z);
			maxZ = (std::max)(maxZ, z);
		}
	}
		
	//申請buf
	size_t dstBufNum = (size_t)bufWidth * bufHeight * dstBand;
	GByte *dstBuf = new GByte[dstBufNum];
	memset(dstBuf, 0, dstBufNum*sizeof(GByte));

	//設定方向：平行光
	double solarAltitude = 45.0;
	double solarAzimuth = 315.0;
	
	//
	double Zenith_rad = osg::DegreesToRadians(90 - solarAltitude);
	double Azimuth_math = 360.0 - solarAzimuth + 90;
	if (Azimuth_math >= 360.0)
	{
		Azimuth_math = Azimuth_math - 360.0;
	}	
	double Azimuth_rad = osg::DegreesToRadians(Azimuth_math);
		
	//a b c
	//d e f
	//g h i
	double z_factor = 2;
	double alpha = 0.3;		//A不透明度 α*A+(1-α)*B

	//
	for (int yi = 1; yi < bufHeight-1; yi++)
	{
		for (int xi = 1; xi < bufWidth-1; xi++)
		{			
			size_t e = (size_t)bufWidth * yi + xi;
			size_t f = e + 1;
			size_t d = e - 1;

			size_t b = e - bufWidth;
			size_t c = b + 1;
			size_t a = b - 1;

			size_t h = e + bufWidth;
			size_t i = h + 1;
			size_t g = h - 1;
			
			float tmpBuf[9] = { imgBuf[a], imgBuf[b], imgBuf[c], imgBuf[d], imgBuf[e], imgBuf[f], imgBuf[g],imgBuf[h], imgBuf[i] };
			double Hillshade = CalHillshade(tmpBuf, Zenith_rad, Azimuth_rad, dx, -dy, z_factor);
			GByte value = (GByte)(std::min(std::max(Hillshade, 0.0), 255.0));

			int index = GetColorIndex(imgBuf[e], minZ, maxZ);
			GByte rgb[3] = { (GByte)tableRGB[index].R, (GByte)tableRGB[index].G, (GByte)tableRGB[index].B };

			for (int ib = 0; ib < dstBand; ib++)
			{
				size_t n = (size_t)bufWidth * dstBand * yi + dstBand * xi + ib;
				double v = value * alpha + (1 - alpha) * rgb[ib];
				dstBuf[n] = (GByte)(std::min)((std::max)(v, 0.0), 255.0);
				//dstBuf[n] = (GByte)value;
			}	
		}
	}


	//寫入
	dst->RasterIO(GF_Write, 0, 0, bufWidth, bufHeight, dstBuf, bufWidth, bufHeight,
		GDT_Byte, dstBand, nullptr, dstBand*dstDepth, bufWidth*dstBand*dstDepth, dstDepth);
	
	//釋放
	delete[] imgBuf;
	imgBuf = nullptr;

	//釋放
	delete[] dstBuf;
	dstBuf = nullptr;

	//
	GDALClose(dst);
	dst = nullptr;

	GDALClose(img);
	img = nullptr;

	return 0;
}

最終實作的效果圖如下，可以看到确實實作了高程越低越藍，越高越紫的暈渲效果，同時具有深度感，能看清山脊地貌，效果與ArcMap基本一緻，隻是配色效果不同。

關于DEM的地貌暈渲圖的實作暫時告一段落了。應該來說還是有些模糊不清的地方，在查閱資料的時候就有所感覺，現在關于GIS的基礎原理資料總是不太清晰，原理概念一堆，但就是了解不了。但如果有新的問題或者發現，希望看到這幾篇文章的朋友能批評指正下。

[1]. ArcGIS制圖手冊(3-2)山體陰影和暈渲

[2]. RGB顔色插值漸變原理及算法

[3]. 兩個RGBA四通道顔色的疊加計算方法與代碼實作

使用GDAL實作DEM的地貌暈渲圖(三)

繼續閱讀

geotrellis-栅格資料切片前言步驟總結

geotrellis入門01——跑例子geotrellis-landsat-tutorial

geotrellis入門03——專業術語

矢量切片

python 利用arcpy批量裁剪tif檔案，并生成tif檔案

從招聘資料看GIS就業

GIS—手把手用ArcGIS打開osm檔案

地圖資料服務GIS相關

空間分析判斷以及坐标轉換工具

proj4測繪坐标系轉換

玩轉地圖投影

GIS | 坐标系統與地圖投影坐标系統三個半概念地圖投影投影坐标系統

GlobalMapper20如何導入文本類（txt，csv，excel）的坐标資料序：一、如何導入離散的坐标點二、離散的坐标點如何生産互相連接配接的線三、離散的坐标點如何生成範圍面

GlobalMapper20正方形的資料變成了矩形，變形了？一、問題描述二、問題分析：三、顯示方式調整

創業,不能兼職(69)---空間資料和屬性資料

坐标系統和投影變換在桌面産品中的應用