1. 原理

之前在《使用GDAL实现DEM的地貌晕渲图(一)》和《使用GDAL实现DEM的地貌晕渲图(二)》这两篇文章中详细介绍了DEM生成地貌晕渲图的原理与实现。不过之前生成的都是晕渲强度值对应的灰度图，而实际的应用过程中都会将DEM晕渲成彩色图。

1) ArcMap生成彩色晕渲图

可以通过ArcMap的做法来参考如何生成彩色晕渲图（参考[1]），在ArcMap中生成彩色晕渲图的步骤如下：

通过山体阴影工具生成灰度晕渲图，这一点与前面文章介绍的相一致。
然后在原DEM图的显示中，选择最大最小拉伸显示，然后选择一个合适的彩色色带赋值。
最后，将步骤一的灰度晕渲图设置一定的透明度，叠加到步骤二的彩色图上，就生成了最终具有立体感的彩色晕渲图。

ArcMap生成的彩色晕渲图：

2) 彩色色带赋值

不难发现，生成彩色晕渲图的关键是第二步：要选取合适的色带，让色带根据对应的高程赋值。查阅了不少的资料，这个色带应该没有固定合适通用的模板，是需要自己根据具体的需要调整的。比如，海平面可以赋值成蓝色；高山山顶赋值成白色；戈壁沙漠赋值成黄色；草原森林赋值成绿色，这些地貌特征都具有一定的高程相关性，可以根据一定的绝对高程区间赋值。

我这里采取的做法还是跟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];}
}

3) 颜色叠加

第一步和第二步分别生成了晕渲强度图和高程彩色色带图，第三步就是将两者的颜色叠加，生成最终的效果图。其实颜色叠加的原理特别简单，对于晕渲强度图的像素值A，令其不透明度为α；对应的高程彩色色带图的像素值B，那么最后叠加的像素值 C=α*A+(1-α)B。可以这么理解:光线到达A，由于透光性只产生了αA的效果，还有(1-α)强度的光线射到B，又产生了(1-α)B的效果，两者叠加就是αA+(1-α)*B。

2. 实现

继续改造之前的代码，最终的实现过程如下：

#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);}//申请bufdepth = 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];           //左上角点坐标Xdouble dx = padfTransform[1];         //X方向的分辨率double startY = padfTransform[3];             //左上角点坐标Ydouble 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);}}//申请bufsize_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 idouble 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基本一致，只是配色效果不同。

3. 结语

关于DEM的地貌晕渲图的实现暂时告一段落了。应该来说还是有些模糊不清的地方，在查阅资料的时候就有所感觉，现在关于GIS的基础原理资料总是不太清晰，原理概念一堆，但就是理解不了。但如果有新的问题或者发现，希望看到这几篇文章的朋友能批评指正下。

4. 参考

[1]. ArcGIS制图手册(3-2)山体阴影和晕渲
[2]. RGB颜色插值渐变原理及算法
[3]. 两个RGBA四通道颜色的叠加计算方法与代码实现

使用GDAL实现DEM的地貌晕渲图(三)相关推荐

ArcGIS中利用DEM制作立体晕渲图的说明
DEM作为4D产品之一,在测绘生产中具有重要作用.利用DEM数据做作晕渲图,并对其进行分层设色,可以很好的反映一个地方的地形.下面以黑龙江省为例,采用90m分辨率的DEM数据,在ArcGIS 10.2 ...
ArcGIS 10.2晕渲图+旋转图制作
晕渲图-通过模拟实际地面本影与落影的方法反映实际地形起伏特征的一种重要的地形图.晕渲图是DEM地表形态表达的一种形式,它通过设置光源的高度角和方位角更形象或者更符合人类视觉的方式展示一个地区的地形.通 ...
OMG！这么优秀的晕渲图，原来四步就可以完成！
晕渲图是表达地形最常见的一种形式,它通过阴影和颜色渐变来展示地表的起伏变化,具有很好的立体感.SuperMap中提供的全球晕渲图地形清晰,配色浅淡,用户可以直接作为底图使用,在上面叠加其它的数据.那么 ...
硬核干货！全球数字高程数据（DEM）详解，还有地形晕渲、等高线等
1 基本概念 DEM是数字高程模型的英文简称(Digital Elevation Model),是研究分析地形.流域.地物识别的重要原始资料.由于DEM 数据能够反映一定分辨率的局部地形特征,因此通过 ...
Matlab下地形图绘图包m_map绘制晕渲（shaded relief）地形图
1. 简介晕渲一词源自绘画,指的是用水墨或颜色渐次浓淡烘染物象,使分出阴阳向背的绘画技法(https://baike.so.com/doc/2078427-2198654.html). 地理学中晕渲 ...
美国高清晕渲地形图分享，每一幅都值得珍藏
以下三维高清晕渲图来源于外网,由王石头下载.整理,并分享给3S技术之家的粉丝朋友! 全都是高清大图,总文件大小达1.4GB. 图片格式有Tiff及jpg,其中Tiff格式带有坐标投影. 原始出处来自于 ...
GDAL使用DEM数据计算山体阴影（Hillshade）
零. 前言说起Hillshade,其实就是模拟太阳光照射地形所引起的明暗对比,然后来对地形图进行渲染,使之看起来具有立体效果的一种方式,常用于地图的渲染,如表1所示,具体的可以参考文献 ...
使用GDAL对DEM进行彩色渲染
在之前的博客中,我们已经看到了gdal对dem数据的强大的处理功能,其中除了坡度坡向等,还有一个比较厉害的,就是使用DEM生成一个彩色的图像.之前关于这方面也翻译了几篇博客,详见<使用GDAL对 ...
GDAL使用DEM数据计算地形指数
零. 前言本文是接上文<GDAL使用DEM数据计算坡度坡向>,还是一篇关于DEM计算地形指数的一篇文章.这里所要计算的地形指数主要包括以下三个指数:地形耐用指数(Terra ...
GDAL使用DEM数据计算坡度坡向
零. 前言之前写过一个3×3的通用模板算子函数的博客<基于GDAL的一个通用的3×3模板函数>,网址:http://blog.csdn.net/liminlu0314/ar ...

使用GDAL实现DEM的地貌晕渲图(三)

文章目录