做地形可视化这几年我最大的感触是大多数人不是不会用Python画等高线而是画出来的东西只能算“示意图”离“成图”还差很远。只调一次ax.contour()得到的是一张平铺在线条里的抽象图真正能拿上台面的等高线阴影图是在等高线下面叠了一层由光照模拟出来的山体阴影也就是地图上常见的那种立体地貌渲染效果。这篇文章我用Python从零手搓一张等高线阴影图把从数据准备、Hillshade算法实现、等高线调优到最后的图层叠加完整走一遍。适合已经有Python和NumPy基础、需要做科研图、工程图或地理可视化作品但一直没把成图逻辑理清楚的朋友。1. 一张等高线阴影图里到底有什么1.1 三层信息各管什么先拆解一下成图结构。等高线阴影图不是单一图层而是三层信息的叠加阴影底图Hillshade模拟太阳光从某个方向照射地表用灰度表达山体明暗。它不提供具体的数字信息但能在半秒钟内让你感知整个区域的地形起伏哪里是山脊、哪里是沟谷一目了然。分层设色Hypsometric Tint把高程范围切分成若干区段给不同区段填上不同的颜色。低海拔用绿、中海拔用黄、高海拔用棕或白是经典的地形配色习惯。等高线Contour把相同高程的点连成线附带高程标注是这三层里唯一能够“定量读取”高程信息的要素。三层信息各自独立合在一起却产生明显的一加一大于二效果底图负责立体感设色负责高程整体趋势等高线负责精确细节。这也是为什么地形图出版领域至今仍保留这种组合形式。1.2 为什么一定要加阴影图纯等高线图的问题在于平面感太强。等高线密的地方是陡坡但到底多陡、山形走向如何普通读者很难脑补出来。阴影图恰好把这种脑补过程直接可视化光照方向一致的斜坡是亮的背光面是暗的山脊线通常落在明暗交界处。有了阴影地形就“站起来”了。有人会问用plot_surface画三维曲面不是更有立体感三维图适合交互式展示一旦落到静态图片里透视角度的畸变反而会遮挡重要细节而且无法精确读取高程。等高线阴影图本质上是二维平面图信息密度更高打印输出也不依赖视角这是它在地质、测绘、论文配图场景里经久不衰的原因。1.3 技术选型为什么是Python而不是现成的GIS软件ArcGIS、QGIS都能直接生成阴影图和等高线成品效果很好。但我个人还是坚持用Python做理由很实际第一脚本可复现数据更新后改个路径重跑一遍就行第二可以自由控制每一个视觉参数尤其是配色、等高线间距这些出版细节GUI里的选项反而受限第三后期要叠加气象、地质、水文数据时脚本方案能无缝衔接。当然代价是需要自己处理坐标系、NoData这些麻烦事但看完这篇文章你就有底了。2. 数据准备自己造DEM还是去下载真实地形2.1 自产自足用NumPy合成一个小地形教学场景最忌讳一上来就从网上下数据GDAL装不上、文件格式不对这些坑会直接消磨掉学习热情。我习惯先手工构造一个DEM把流程跑通再换真实数据。构造DEM最朴素的方法是叠加几个二维高斯函数模拟起伏不平的山地import numpy as np y, x np.mgrid[0:100, 0:120] dem ( 120 * np.exp(-((x - 50) ** 2 (y - 40) ** 2) / 800) 90 * np.exp(-((x - 80) ** 2 (y - 75) ** 2) / 500) 60 * np.exp(-((x - 25) ** 2 (y - 80) ** 2) / 300) )这个数组的数值大致在0到200之间可以当作相对高程米。等高线阴影图算法本身不关心高程绝对大小它关心的是高程变化率也就是梯度所以这样的模拟数据完全够用。为什么要用高斯函数因为它在各个方向都平滑既不会出现离谱的尖峰也不会出现大片平地导致梯度为零。多个高斯相加后山与山之间的鞍部、坡向变化都很自然非常适合测试阴影算法。2.2 真实世界公开DEM数据源与rasterio读取如果你要画真实区域优先推荐几个数据源USGS EarthExplorerSRTM 30米、ASTER GDEM等都在这里覆盖全球。OpenTopography按地区直接下载处理好的GeoTIFF对新手最友好。AWS Open DataTerrain Tiles按瓦片读取适合做程序化接入。拿到GeoTIFF后读取通常用rasterio这一行import rasterio with rasterio.open(dem.tif) as src: dem src.read(1).astype(float) transform src.transform res_x transform.a # 东西方向分辨率单位米 res_y transform.e # 南北方向分辨率单位米注意通常是负值 print(dem.shape, abs(res_x), abs(res_y))这里有个容易被忽视的点transform.a是像素宽度transform.e是行方向的分辨率它通常是负数表示地理坐标Y轴向下。计算坡度前记得取绝对值否则梯度方向全反。2.3 数据预处理NoData、单位与裁剪真实DEM几乎都有无效区域常见的是海洋和边缘NoData读出来是-32768这种哨兵值。如果不处理后面算梯度时会算出一圈离谱的数值导致等高线在边界处乱飞。最简单的处理有两种# 方案1替换成NaN这样后续统计都会自动忽略 dem[dem -10000] np.nan # 方案2直接掩膜成布尔数组计算梯度前填0 mask dem -10000 dem np.where(mask, dem, np.nan)需要说明的是np.gradient遇到NaN会把邻近一圈都变NaN所以更好的做法是先对NaN区域做插值填充或者在计算阴影图后处理NaN。如果你只是画图把NaN区域在imshow里直接隐藏也说得过去。裁剪区域我一般提一下不展开大范围数据先用目标经纬度或边界文件切一块减小计算量。rasterio.mask可以做矢量边界裁剪嫌重可以直接用数组切片dem[y0:y1, x0:x1]配合范围计算只要心里清楚切片坐标对应的是行列号。3. 手写Hillshade光照阴影的计算原理与NumPy实现3.1 从高程到坡度为什么先算梯度阴影图的核心是模拟太阳光照到地表后的明暗分布。要算光照先要知道地表每一处有多倾斜、朝哪个方向倾斜这两个量被GIS称作坡度Slope和坡向Aspect。对DEM数组来说坡度本质上是高程在空间上的变化率也就是梯度。NumPy直接给你现成函数gy, gx np.gradient(dem, res_y, res_x)注意返回顺序第一个是沿着数组行方向南北向的梯度第二个是沿着列方向东西向的梯度。如果分辨率是30米那么res_x30, res_y30传给函数的就是每个像素对应的实际距离。这一步特别关键分辨率不对坡度就完全失真。坡度角用梯度模长的反正切计算slope_rad np.arctan(np.sqrt(gx * gx gy * gy))坡向则是两个梯度分量合成后的方向角aspect_rad np.arctan2(-gy, gx)这里有个坐标系约定问题。DEM数组的行方向是南向北但在图像坐标系里y方向向下所以取负号把方向转入地图坐标。如果你发现渲染出来的阴影方向整体反了优先检查的就是这个负号。不同数据投影或处理流程可能导致系统性的方向偏差不要迷信公式要以视觉结果为准。3.2 光照模型太阳方位与入射角的数学表达有了坡度和坡向接下来就是模拟太阳。地图学的标准Hillshade模型只用两个参数控制光照太阳高度角Altitude默认取45度越低影子拉得越长立体感越强但暗部会更多。太阳方位角Azimuth默认取315度也就是西北方向的光源。这个值在传统地图制图中很常用因为阴影落在东南方向不容易在视觉上跟地图文字的排布冲突。单个像素接受的直接光照强度可以写成一个余弦型公式shaded sin(alt) * cos(slope) cos(alt) * sin(slope) * cos(az - aspect)这个公式的含义很直白当太阳正好在头顶alt90°时任何朝向的地表都只受坡度影响不受坡向影响因为sin(alt)等于1、cos(alt)等于0当太阳斜着照过来时地表朝向与光源方向越一致第二项越大就越亮。完全背光的地方整个表达式接近0甚至为负也就是阴影区。3.3 NumPy向量化实现与参数细节把上面三步整合成完整函数核心代码只有几行def hillshade(dem, res_x10.0, res_y10.0, azimuth315.0, altitude45.0): 根据DEM计算山体阴影灰度图值域0~255。 gy, gx np.gradient(dem, res_y, res_x) slope_rad np.arctan(np.sqrt(gx * gx gy * gy)) aspect_rad np.arctan2(-gy, gx) azimuth_rad np.radians(azimuth) altitude_rad np.radians(altitude) shaded ( np.sin(altitude_rad) * np.cos(slope_rad) np.cos(altitude_rad) * np.sin(slope_rad) * np.cos(azimuth_rad - aspect_rad) ) return np.clip(shaded, 0, 1) * 255.0重点看最后一行我先求光照强度的余弦值理论上值域是[-1,1]但负值意味着完全背光。地图制图惯例中背光区直接设为黑色0所以做了一次clip(0,1)。如果你想要更柔和的阴影过渡也可以改成归一化映射return 255.0 * (shaded 1.0) / 2.0这样即使完全背光也只是中灰色不会出现死黑。两者肉眼观感差别很大我建议按图片用途取舍底图输出用前一种视觉冲击强要叠加大量注记用后一种背景不至于太抢。res_x和res_y默认给了10对应合成DEM的理想分辨率。真实数据一定要把rasterio读出来的分辨率传进来宁可不传也不要把单位搞错。如果DEM本身是60米分辨率却填了10坡度会被放大6倍阴影图会看起来极其生硬。4. 等高线生成参数从默认到出版级的调优过程4.1 contour与contourf的分工很多教程把contour和contourf混着用其实职责完全不同。contourf画的是填充色块负责的是“分层设色”那一层视觉效果。contour画的是线条负责的是精确高程边界。画等高线阴影图时两者通常都会用但参数侧重点不同。填充层关心色带范围和透明度线层关心线宽、线色和标注。这里提前提醒一句如果只画线不填充图会显得寡淡如果只填充不画线读者无法快速读出准确高程。两个一起上才能发挥第一章节里讲的三层信息优势。4.2 等高线间距的选择策略等高线间距是整张图最容易出效果也最容易翻车的参数。间距过大地形细节丢失间距过小线条糊成一团。我通常不让间距依赖手动拍脑袋而是用数据范围自动计算vmin np.nanmin(dem) vmax np.nanmax(dem) step 20 levels np.arange(np.floor(vmin / step) * step, np.ceil(vmax / step) * step step, step)这个思路是先把高程下边界取整到20的倍数上边界取整到20的倍数再按20米步长生成等高线列表。好处是输出的等高线都是“20、40、60”这类整洁数字而不是“37.4、41.2”这种没意义的带小数级别。步长到底取多少我的经验是等高线数量控制在10到20条之间视觉效果最好。你可以这样快速估算rough_step (vmax - vmin) / 15 step round(rough_step / 10) * 10 # 取整到10的倍数如果要画陡峭山地等高线会自然堆积在陡坡处这时宁可选稍大的间距否则输出成矢量图后线条过多、文件尺寸也会失控。4.3 标注与线型的细节处理等高线标注重灾区有两个一是标注数字跟线交叉处被压住看不清二是标注字体大小不匹配图片缩放比例。clabel的inlineTrue是最值得开的参数它会在标注数字处把等高线切断留出干净的白色背景阅读性提升非常明显cs ax.contour(dem, levelslevels, colorsblack, linewidths0.7, alpha0.8, extentextent, originupper, zorder3) ax.clabel(cs, inlineTrue, fontsize8, fmt%d)fmt%d是第二个关键细节。默认标注会写成浮点数可能出现50.0这种多余尾巴。用整数格式%d只要级别本身就是整数标注就会是干净的50。如果你的高程有小数位可以改成fmt%.1f但配图审美上最好统一。线型方面常规做法是主线用深色实线遇到一些特殊级别比如每5条线加深一次可以做成叠加方案。这个我在第5章的完整代码里会演示。5. 图层面板叠加顺序、透明度与配色决定成败5.1 图层顺序就是视觉的优先级到这一步相当于做菜的最后组装。图层顺序非常讲究一般是从远到近底图阴影最先画紧跟着是填充色最后是等高线线条和标注。import matplotlib.pyplot as plt extent [0, 120, 0, 100] hs hillshade(dem) fig, ax plt.subplots(figsize(10, 8)) ax.imshow(hs, cmapgray, extentextent, originupper) ax.contourf(dem, levelslevels, cmapterrain, extentextent, originupper, alpha0.6, zorder2) cs ax.contour(dem, levelslevels, colorsblack, linewidths0.7, alpha0.8, extentextent, originupper, zorder3) ax.clabel(cs, inlineTrue, fontsize8, fmt%d) plt.show()这里的originupper和extent必须严格一致否则图像上下翻转或者坐标错位第6章第一个坑就是这个。灰度阴影作为底图用imshow画完后contourf和contour的坐标范围必须与该extent完全对应。5.2 配色方案默认terrain可以但你可以做得更好Matplotlib自带terrain和gist_earth两张地形色带直接用的效果其实不算差但一个通病是低海拔区域颜色偏暗、和灰底图叠在一起不够清爽。我更推荐两种方案一是安装cmocean库它的topo色带就是为地形可视化设计的二是自制分段色带完全控制颜色拐点from matplotlib.colors import LinearSegmentedColormap import matplotlib as mpl colors [#2c7a3e, #8fbf57, #e0c774, #b58a53, #8d5a3a, #f5f5f5] custom_cmap LinearSegmentedColormap.from_list(terrain_custom, colors) norm mpl.colors.Normalize(vminvmin, vmaxvmax) ax.contourf(dem, levelslevels, cmapcustom_cmap, normnorm, extentextent, originupper, alpha0.6, zorder2)从低到高依次是深绿、浅绿、浅黄、黄褐、深褐、白这是地形图很常见的一套印象派配色。透明度alpha0.6是个经验值太低底图阴影盖过颜色层次全无太高颜色盖住灰度阴影立体感又没了。0.5到0.7之间都可以建议输出前多试两个值。5.3 指北针、比例尺与色标的组合一张出版级地图还需要地图要素衬托。色标直接用colorbarcbar fig.colorbar(cs, axax, shrink0.8, labelElevation (m))比例尺优先用matplotlib_scalebar这个第三方控件from matplotlib_scalebar.scalebar import ScaleBar ax.add_artist(ScaleBar(1, unitsm, length_fraction0.2))ScaleBar(1, unitsm)表示一个数据单位对应1米如果你的DEM分辨率是10米这里第一个参数就填10并在文档里说明清楚。指北针不用额外库用annotate几行就能画ax.annotate(N, xy(0.95, 0.92), xycoordsaxes fraction, hacenter, fontsize12, fontweightbold) ax.plot(0.95, 0.90, transformax.transAxes, marker^, colorblack)这段代码把指北针固定在右上角位置用坐标轴比例表示不随数据范围变化。注意箭头在文字下方先后顺序别写反。5.4 完整脚本串一遍把前面所有片段组合成一个可以独立运行的最小脚本方便直接抄作业import numpy as np import matplotlib.pyplot as plt from matplotlib.colors import LinearSegmentedColormap y, x np.mgrid[0:100, 0:120] dem (120 * np.exp(-((x - 50) ** 2 (y - 40) ** 2) / 800) 90 * np.exp(-((x - 80) ** 2 (y - 75) ** 2) / 500) 60 * np.exp(-((x - 25) ** 2 (y - 80) ** 2) / 300)) res_x, res_y 10.0, 10.0 extent [0, 120 * res_x, 0, 100 * res_y] def hillshade(dem, res_x, res_y, azimuth315.0, altitude45.0): gy, gx np.gradient(dem, res_y, res_x) slope_rad np.arctan(np.sqrt(gx * gx gy * gy)) aspect_rad np.arctan2(-gy, gx) azimuth_rad np.radians(azimuth) altitude_rad np.radians(altitude) shaded (np.sin(altitude_rad) * np.cos(slope_rad) np.cos(altitude_rad) * np.sin(slope_rad) * np.cos(azimuth_rad - aspect_rad)) return np.clip(shaded, 0, 1) * 255.0 vmin, vmax np.nanmin(dem), np.nanmax(dem) levels np.arange(0, 200, 10) colors [#2c7a3e, #8fbf57, #e0c774, #b58a53, #8d5a3a, #f5f5f5] custom_cmap LinearSegmentedColormap.from_list(terrain_custom, colors) fig, ax plt.subplots(figsize(10, 8)) ax.imshow(hillshade(dem, res_x, res_y), cmapgray, extentextent, originupper) ax.contourf(dem, levelslevels, cmapcustom_cmap, extentextent, originupper, alpha0.6, zorder2) cs ax.contour(dem, levelslevels, colorsblack, linewidths0.7, alpha0.8, extentextent, originupper, zorder3) ax.clabel(cs, inlineTrue, fontsize8, fmt%d) plt.savefig(contour_hillshade.png, dpi300, bbox_inchestight)把这个脚本跑通一张合格的等高线阴影图就出来了。接下来的内容才是真正让你少走弯路的部分。6. 实测踩坑记录四个高频坑与判断依据6.1 坐标范围不一致导致等高线“起飞”我见过最离谱的一次底图阴影完全正常但等高线全部对不上位置有的线条甚至跑到图外。最后定位到原因是imshow的extent写的是[0, 120, 0, 100]而contour的extent写成了[0, 100, 0, 120]行列顺序颠倒了。排查方法很直接在一个已知的高峰位置打一个点看它是否同时落在底图像元和等高线环上。如果点跟底图对齐但跟等高线错位问题几乎一定在extent或origin。另外imshow默认originuppercontour默认数组坐标也是从上往下但一旦混入其他来源的数据就很容易一个上、一个下视觉表现就是等高线上下镜像错位。统一加originupper能解决大多数问题。6.2 阴影图太暗或过亮阴影图刚出来时很多人觉得山体太黑、细节被吞掉。大多数情况下不是算法错了而是光照参数太激进。高度角越低影子越长45度是经典默认值但对比较平缓的地形45度反而让大片区域处于微光状态。我的经验做法是先用altitude45看整体如果暗部比例超过三分之一就把高度角抬到55~60。另一个思路是把最终灰度做一次线性拉伸把原来0~255的区间压缩到40~255相当于天然加了环境光hs hillshade(dem, res_x, res_y) hs 40 hs * (255 - 40) / 255这样就保证最暗的阴影区也保留一定灰度不会出现死黑。环境光强度你可以按需调整40这个值不算拍脑袋是经过多轮对比得到的平衡点太大会让阴影失去层次太小又回到死黑。顺便说一句如果你用clip(0,1)版本但不做拉伸出来的图通常偏暗这是正常现象不是代码写错了。6.3 大范围数据渲染卡顿真实DEM动辄几千像素见方np.gradient其实非常快真正的瓶颈在等高线算法和imshow渲染。整幅图的内存占用和输出时间会随着像素数线性增长到几百万像素时就明显卡了。两个办法一是按显示用途降采样比如dem[::2, ::2]把行列各抽一半像素一下变成四分之一。降采样后记得把res_x和res_y改成原来的两倍否则坡度会翻倍。二是分块绘制。科研绘图常会遇到这种场景不需要全分辨率渲染只要清晰的等高线。我一般把DEM降到1000像素见方以下再出图肉眼很难分辨细节差异但速度能快一个数量级。注意降采样后extent不能变因为地理范围没变。6.4 中文字体乱码Matplotlib默认字体对中文支持极差图里只要出现中文标注就是方块。这个问题说大不大但每回都能绊倒一些人plt.rcParams[font.sans-serif] [SimHei, Noto Sans CJK SC, WenQuanYi Zen Hei] plt.rcParams[axes.unicode_minus] False第一行设置中文字体按你系统里实际存在的字体选一个第二行是为了让负号正常显示不设置的话坐标轴负号有时会变成乱码。如果是纯英文标注这两行可以完全省略。Linux服务器上渲染图时建议先执行fc-list :langzh看看系统里有没有中文字体没有就拿安装包补一个这是服务器上最常见的坑。最后分享一个我自己的使用习惯做展示用的PNG图把dpi设到300且保存时加bbox_inchestight避免出图四周留白做论文配图或后期还要进AI、CAD处理的建议直接存PDF矢量格式等高线和标注放大不糊。脚本化出图的价值在于数据更新后我只需要改DEM路径参数全部自动重算再也不用在GIS软件里反复调一遍设置。这个工作流我用了很久算是从纯手工画图到自动出图之间性价比最高的一段路。