Matlab下地形图绘图包m_map绘制晕渲（shaded relief）地形图

1. 简介

晕渲一词源自绘画，指的是用水墨或颜色渐次浓淡烘染物象，使分出阴阳向背的绘画技法（https://baike.so.com/doc/2078427-2198654.html）。

地理学中晕渲图指的是DEM地表形态表达的一种形式，它通过设置光源的高度角和方位角更形象或者更符合人类视觉的方式展示一个地区的地形。通过晕渲图，可以很好的反映地形地势的变化，有很好的立体感，方便用图者的使用（https://baike.baidu.com/item/%E6%99%95%E6%B8%B2%E5%9B%BE）。

具体可参考维基百科：https://en.wikipedia.org/wiki/Terrain_cartography#Shaded_relief。

m_map作为Matlab下一款不错的地形图绘制工具箱，也带有地形晕渲功能。相关函数有：m_etopo2和m_shadedrelief。

2. 函数说明

先看一下m_shadedrelief函数说明。

>> help m_shadedreliefm_shadedrelief Shaded relief topography in an imagem_shadedrelief(X,Y,Z) presents a shaded relief topography, as wouldbe seen if a 3D model was artifically lit from the side. Slopesfacing the light are lightened, and slopes facing away from the light are darkened. X and Y are horizontal and vertical coordinateVECTORS, and these should be in the same units as the height Z (e.g., all in meters), otherwise the slope angle calculations will be in error. Usage notes:(1) m_shadedrelief is a replacement for a low-level call to IMAGE displaying a true-colour image so it MUST be preceded by COLORMAP and CAXIS calls.  (2) m_shadedrelief probably is most useful as a backdrop to maps witha rectangular outline box - either a cylindrical projection, or someother projection with M_PROJ(...'rectbox','on'). (3) Finally, the simplest way of not running into problems:- if your elevation data is in LAT/LON coords (i.e. in a matrix whereeach row has points with the same latitude, and each column has pointswith the same longitude), use M_PROJ('equidistant cylindrical',...)- if your elevation data is in UTM coords (meters E/N), i.e. in a matrixwhere each row has the same UTM northing and the each column has thesame UTM easting, useM_PROJ('utm',....)m_shadedrelief(...,'parameter',value) lets you set various properties.These are:'coords' : Coordinates of X/Y/Z: 'geog' for lat/lon, Z meters,  (default)'map'  for X/Y map coordinates, Z meters'Z' if X/Y/Z are all in same units (e.g., meters)'lightangle' : true direction (degrees) of light source (default -45, i.e. from the north-west)'gradient': Shading effects increase with slope angle until slopes reach this value (in degrees), and are held constant for higher slopes (default 10). Reduce for smoother surfaces.'clipval' : Fractional change in shading for slopes>='gradient'.0 means no change, 1 means saturation to white or blackif slope is facing directly towards or away from light source, (default 0.9).'nancol'  : RGB colour of NaN values (default [1 1 1]);'lakecol' : RGB colour of lakes (flat sections) (default NaN)If set to NaN lakes are ignored.IM=m_shadedrelief(...) returns a handle to the image.[SR,X,Y]=m_SHADEDRELIEF(...) does not create an image but only returnsthe  true-color matrix SR of the shaded relief, as well as the X/Yvectors needed to display it.Example:load toposubplot(2,1,1);  % Example without itimagesc(topolonlim,topolatlim,topo);caxis([-5000 5000]);colormap([m_colmap('water',64);m_colmap('gland',64)]);set(gca,'ydir','normal');subplot(2,1,2);  % Example with itcaxis([-5000 5000]);colormap([m_colmap('water',64);m_colmap('gland',64)]);m_shadedrelief(topolonlim,topolatlim,topo,'gradient',5e2,'coord','Z');axis tight

函数基本格式如下：

m_shadedrelief('coords',OPTION)

其中coords表示输入的经纬度和高程坐标，OPTION表示可选项。
可选项主要用到2个：lightangle和gradient，分别表示光源方位角和坡度（梯度）阈值。

光源方位角lightangle从地理正北方从0°~360°顺时针旋转，或者逆时针旋转，用负角度表示。比如，m_map的m_shadedrelief函数的光源方位角默认值是-45°，即光源位于正西北方向。

理论上，晕渲图应该设置光源的方位角和高度角。但是m_map通过设置地形坡度阈值来间接地展示不同高度角的光源照明。通过对大坡度地形设置阴影效果来突出地形起伏变化。坡度角越大，阴影效果越明显，直到坡度角达到设定值（以度为单位），并且对于较高的坡度保持不变（默认值为10）。

另外，需要说明的是，m_etopo2('shadedrelief')函数和m_shadedrelief在某种程度上等价，但是m_etopo2('shadedrelief')更加简单。

下面我们看一下如何用这两个函数绘制晕渲地形图。

3. 实例说明

3.1 例1

使用m_etopo2('shadedrelief')，默认参数（光源方位角-45°，坡度阈值10）。

figname='shadedrelief1';m_proj('mercator','long',[114 124],'lat',[20 26]);caxis([-5000 0]); % 设置显示高程范围，只显示海洋部分colormap((m_colmap('blues',200)));  % 设置色标为m_map自带的m_colmap('blues')
hc=colorbar;
set(get(hc,'title'),'string','Elevation(m)')  %  色标单位m_etopo2('shadedrelief');m_gshhs('ic','patch',[.6 .6 .6]); % 设置中等分辨率海岸线填充m_grid('box','fancy','gridlines','no','fontsize',14);  % 轴设置set(gcf,'position',[100 100 800 500]) %图形大小设置
print('-dpng','-r600',[figname,'.png'])  % 导出png图片

3.2 例2

使用m_etopo2('shadedrelief')，光源方位角设置为45°，默认坡度阈值10°。

figname='shadedrelief2';m_proj('mercator','long',[114 124],'lat',[20 26]);caxis([-5000 0]); % 设置显示高程范围，只显示海洋部分colormap((m_colmap('blues',200)));  % 设置色标为m_map自带的m_colmap('blues')
hc=colorbar;
set(get(hc,'title'),'string','Elevation(m)')  %  色标单位m_etopo2('shadedrelief','lightangle',45); % 设置光源方位角为45°m_gshhs('ic','patch',[.6 .6 .6]); % 设置中等分辨率海岸线填充m_grid('box','fancy','gridlines','no','fontsize',14);  % 轴设置set(gcf,'position',[100 100 800 500]) %图形大小设置
print('-dpng','-r600',[figname,'.png'])  % 导出png图片

3.3 例3

使用m_etopo2('shadedrelief')，光源方位角设置为45°，坡角阈值设置为1°。


figname='shadedrelief3';m_proj('mercator','long',[114 124],'lat',[20 26]);caxis([-5000 0]);    % 设置显示高程范围，只显示海洋部分colormap((m_colmap('blues',200)));  % 设置色标为m_map自带的m_colmap('blues')
hc=colorbar;
set(get(hc,'title'),'string','Elevation(m)')  %  色标单位m_etopo2('shadedrelief','lightangle',45,'gradient',1); % 设置光源方位角为45°,坡角阈值为1m_gshhs('ic','patch',[.6 .6 .6]);   % 设置中等分辨率海岸线填充m_grid('box','fancy','gridlines','no','fontsize',14);  % 轴设置set(gcf,'position',[100 100 800 500]) %图形大小设置
print('-dpng','-r600',[figname,'.png'])  % 导出png图片

3.4 例4

使用m_shadedrelief，默认设置（光源方位角设置为-45°，坡角阈值设置为10°）。

figname='shadedrelief4';m_proj('mercator','long',[114 124],'lat',[20 26]);
[ELEV,LON,LAT]=m_etopo2([114 124 20 26]);
caxis([-5000 3000])
% %caxis要放在colormap之前，colormap要放在m_shadedrelief之前
% 使用两种色标区分陆地海洋，范围比例要和caxis匹配
colormap([m_colmap('blues',200);m_colmap('gland',120)])
hc=colorbar;
set(get(hc,'title'),'string','Elevation(m)')  %  色标单位m_shadedrelief(LON(1,:),LAT(:,1),ELEV)
m_gshhs('ic','color','k')m_grid('box','fancy','tickdir','in','gridlines','no','fontsize',14)set(gcf,'position',[100 100 800 500])
print('-dpng','-r600',[figname,'.png'])  % 导出png图片

3.5 例5

使用m_shadedrelief，光源方位角设置为45°，坡角阈值设置为3°。

figname='shadedrelief5';m_proj('mercator','long',[114 124],'lat',[20 26]);
[ELEV,LON,LAT]=m_etopo2([114 124 20 26]);
caxis([-5000 3000])
% %caxis要放在colormap之前，colormap要放在m_shadedrelief之前
% 使用两种色标区分陆地海洋，范围比例要和caxis匹配
colormap([m_colmap('blues',200);m_colmap('gland',120)])
hc=colorbar;
set(get(hc,'title'),'string','Elevation(m)')  %  色标单位m_shadedrelief(LON(1,:),LAT(:,1),ELEV,'lightangle',45,'gradient',3)
m_gshhs('ic','color','k')m_grid('box','fancy','tickdir','in','gridlines','no','fontsize',14)set(gcf,'position',[100 100 800 500])
print('-dpng','-r600',[figname,'.png'])  % 导出png图片

3.6 例6

在例5（使用m_shadedrelief，光源方位角设置为45°，坡角阈值设置为3°。）的基础上，添加河流，顺便改一下轴类型。

figname='shadedrelief6';
figure
m_proj('mercator','long',[114 124],'lat',[20 26]);
[ELEV,LON,LAT]=m_etopo2([114 124 20 26]);
caxis([-5000 3000])
% %caxis要放在colormap之前，colormap要放在m_shadedrelief之前
% 使用两种色标区分陆地海洋，范围比例要和caxis匹配
colormap([m_colmap('blues',200);m_colmap('gland',120)])
hc=colorbar;
set(get(hc,'title'),'string','Elevation(m)')  %  色标单位m_shadedrelief(LON(1,:),LAT(:,1),ELEV,'lightangle',45,'gradient',3)
m_gshhs('ic','color','k')   % 中等分辨率海岸线
m_gshhs('ir4','color','b')  % 中等分辨率河流m_grid('box','fancy','tickdir','out','gridlines','no','fontsize',14)set(gcf,'position',[100 100 800 500])
print('-dpng','-r600',[figname,'.png'])  % 导出png图片

暂时先写这些。

END