做地形可視化這幾年我最大的感觸是大多數(shù)人不是不會用Python畫等高線而是畫出來的東西只能算“示意圖”離“成圖”還差很遠。只調(diào)一次ax.contour()得到的是一張平鋪在線條里的抽象圖真正能拿上臺面的等高線陰影圖是在等高線下面疊了一層由光照模擬出來的山體陰影也就是地圖上常見的那種立體地貌渲染效果。這篇文章我用Python從零手搓一張等高線陰影圖把從數(shù)據(jù)準備、Hillshade算法實現(xiàn)、等高線調(diào)優(yōu)到最后的圖層疊加完整走一遍。適合已經(jīng)有Python和NumPy基礎(chǔ)、需要做科研圖、工程圖或地理可視化作品但一直沒把成圖邏輯理清楚的朋友。1. 一張等高線陰影圖里到底有什么1.1 三層信息各管什么先拆解一下成圖結(jié)構(gòu)。等高線陰影圖不是單一圖層而是三層信息的疊加陰影底圖Hillshade模擬太陽光從某個方向照射地表用灰度表達山體明暗。它不提供具體的數(shù)字信息但能在半秒鐘內(nèi)讓你感知整個區(qū)域的地形起伏哪里是山脊、哪里是溝谷一目了然。分層設(shè)色Hypsometric Tint把高程范圍切分成若干區(qū)段給不同區(qū)段填上不同的顏色。低海拔用綠、中海拔用黃、高海拔用棕或白是經(jīng)典的地形配色習慣。等高線Contour把相同高程的點連成線附帶高程標注是這三層里唯一能夠“定量讀取”高程信息的要素。三層信息各自獨立合在一起卻產(chǎn)生明顯的一加一大于二效果底圖負責立體感設(shè)色負責高程整體趨勢等高線負責精確細節(jié)。這也是為什么地形圖出版領(lǐng)域至今仍保留這種組合形式。1.2 為什么一定要加陰影圖純等高線圖的問題在于平面感太強。等高線密的地方是陡坡但到底多陡、山形走向如何普通讀者很難腦補出來。陰影圖恰好把這種腦補過程直接可視化光照方向一致的斜坡是亮的背光面是暗的山脊線通常落在明暗交界處。有了陰影地形就“站起來”了。有人會問用plot_surface畫三維曲面不是更有立體感三維圖適合交互式展示一旦落到靜態(tài)圖片里透視角度的畸變反而會遮擋重要細節(jié)而且無法精確讀取高程。等高線陰影圖本質(zhì)上是二維平面圖信息密度更高打印輸出也不依賴視角這是它在地質(zhì)、測繪、論文配圖場景里經(jīng)久不衰的原因。1.3 技術(shù)選型為什么是Python而不是現(xiàn)成的GIS軟件ArcGIS、QGIS都能直接生成陰影圖和等高線成品效果很好。但我個人還是堅持用Python做理由很實際第一腳本可復現(xiàn)數(shù)據(jù)更新后改個路徑重跑一遍就行第二可以自由控制每一個視覺參數(shù)尤其是配色、等高線間距這些出版細節(jié)GUI里的選項反而受限第三后期要疊加氣象、地質(zhì)、水文數(shù)據(jù)時腳本方案能無縫銜接。當然代價是需要自己處理坐標系、NoData這些麻煩事但看完這篇文章你就有底了。2. 數(shù)據(jù)準備自己造DEM還是去下載真實地形2.1 自產(chǎn)自足用NumPy合成一個小地形教學場景最忌諱一上來就從網(wǎng)上下數(shù)據(jù)GDAL裝不上、文件格式不對這些坑會直接消磨掉學習熱情。我習慣先手工構(gòu)造一個DEM把流程跑通再換真實數(shù)據(jù)。構(gòu)造DEM最樸素的方法是疊加幾個二維高斯函數(shù)模擬起伏不平的山地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) )這個數(shù)組的數(shù)值大致在0到200之間可以當作相對高程米。等高線陰影圖算法本身不關(guān)心高程絕對大小它關(guān)心的是高程變化率也就是梯度所以這樣的模擬數(shù)據(jù)完全夠用。為什么要用高斯函數(shù)因為它在各個方向都平滑既不會出現(xiàn)離譜的尖峰也不會出現(xiàn)大片平地導致梯度為零。多個高斯相加后山與山之間的鞍部、坡向變化都很自然非常適合測試陰影算法。2.2 真實世界公開DEM數(shù)據(jù)源與rasterio讀取如果你要畫真實區(qū)域優(yōu)先推薦幾個數(shù)據(jù)源USGS EarthExplorerSRTM 30米、ASTER GDEM等都在這里覆蓋全球。OpenTopography按地區(qū)直接下載處理好的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是行方向的分辨率它通常是負數(shù)表示地理坐標Y軸向下。計算坡度前記得取絕對值否則梯度方向全反。2.3 數(shù)據(jù)預處理NoData、單位與裁剪真實DEM幾乎都有無效區(qū)域常見的是海洋和邊緣NoData讀出來是-32768這種哨兵值。如果不處理后面算梯度時會算出一圈離譜的數(shù)值導致等高線在邊界處亂飛。最簡單的處理有兩種# 方案1替換成NaN這樣后續(xù)統(tǒng)計都會自動忽略 dem[dem -10000] np.nan # 方案2直接掩膜成布爾數(shù)組計算梯度前填0 mask dem -10000 dem np.where(mask, dem, np.nan)需要說明的是np.gradient遇到NaN會把鄰近一圈都變NaN所以更好的做法是先對NaN區(qū)域做插值填充或者在計算陰影圖后處理NaN。如果你只是畫圖把NaN區(qū)域在imshow里直接隱藏也說得過去。裁剪區(qū)域我一般提一下不展開大范圍數(shù)據(jù)先用目標經(jīng)緯度或邊界文件切一塊減小計算量。rasterio.mask可以做矢量邊界裁剪嫌重可以直接用數(shù)組切片dem[y0:y1, x0:x1]配合范圍計算只要心里清楚切片坐標對應(yīng)的是行列號。3. 手寫Hillshade光照陰影的計算原理與NumPy實現(xiàn)3.1 從高程到坡度為什么先算梯度陰影圖的核心是模擬太陽光照到地表后的明暗分布。要算光照先要知道地表每一處有多傾斜、朝哪個方向傾斜這兩個量被GIS稱作坡度Slope和坡向Aspect。對DEM數(shù)組來說坡度本質(zhì)上是高程在空間上的變化率也就是梯度。NumPy直接給你現(xiàn)成函數(shù)gy, gx np.gradient(dem, res_y, res_x)注意返回順序第一個是沿著數(shù)組行方向南北向的梯度第二個是沿著列方向東西向的梯度。如果分辨率是30米那么res_x30, res_y30傳給函數(shù)的就是每個像素對應(yīng)的實際距離。這一步特別關(guān)鍵分辨率不對坡度就完全失真。坡度角用梯度模長的反正切計算slope_rad np.arctan(np.sqrt(gx * gx gy * gy))坡向則是兩個梯度分量合成后的方向角aspect_rad np.arctan2(-gy, gx)這里有個坐標系約定問題。DEM數(shù)組的行方向是南向北但在圖像坐標系里y方向向下所以取負號把方向轉(zhuǎn)入地圖坐標。如果你發(fā)現(xiàn)渲染出來的陰影方向整體反了優(yōu)先檢查的就是這個負號。不同數(shù)據(jù)投影或處理流程可能導致系統(tǒng)性的方向偏差不要迷信公式要以視覺結(jié)果為準。3.2 光照模型太陽方位與入射角的數(shù)學表達有了坡度和坡向接下來就是模擬太陽。地圖學的標準Hillshade模型只用兩個參數(shù)控制光照太陽高度角Altitude默認取45度越低影子拉得越長立體感越強但暗部會更多。太陽方位角Azimuth默認取315度也就是西北方向的光源。這個值在傳統(tǒng)地圖制圖中很常用因為陰影落在東南方向不容易在視覺上跟地圖文字的排布沖突。單個像素接受的直接光照強度可以寫成一個余弦型公式shaded sin(alt) * cos(slope) cos(alt) * sin(slope) * cos(az - aspect)這個公式的含義很直白當太陽正好在頭頂alt90°時任何朝向的地表都只受坡度影響不受坡向影響因為sin(alt)等于1、cos(alt)等于0當太陽斜著照過來時地表朝向與光源方向越一致第二項越大就越亮。完全背光的地方整個表達式接近0甚至為負也就是陰影區(qū)。3.3 NumPy向量化實現(xiàn)與參數(shù)細節(jié)把上面三步整合成完整函數(shù)核心代碼只有幾行def hillshade(dem, res_x10.0, res_y10.0, azimuth315.0, altitude45.0): 根據(jù)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]但負值意味著完全背光。地圖制圖慣例中背光區(qū)直接設(shè)為黑色0所以做了一次clip(0,1)。如果你想要更柔和的陰影過渡也可以改成歸一化映射return 255.0 * (shaded 1.0) / 2.0這樣即使完全背光也只是中灰色不會出現(xiàn)死黑。兩者肉眼觀感差別很大我建議按圖片用途取舍底圖輸出用前一種視覺沖擊強要疊加大量注記用后一種背景不至于太搶。res_x和res_y默認給了10對應(yīng)合成DEM的理想分辨率。真實數(shù)據(jù)一定要把rasterio讀出來的分辨率傳進來寧可不傳也不要把單位搞錯。如果DEM本身是60米分辨率卻填了10坡度會被放大6倍陰影圖會看起來極其生硬。4. 等高線生成參數(shù)從默認到出版級的調(diào)優(yōu)過程4.1 contour與contourf的分工很多教程把contour和contourf混著用其實職責完全不同。contourf畫的是填充色塊負責的是“分層設(shè)色”那一層視覺效果。contour畫的是線條負責的是精確高程邊界。畫等高線陰影圖時兩者通常都會用但參數(shù)側(cè)重點不同。填充層關(guān)心色帶范圍和透明度線層關(guān)心線寬、線色和標注。這里提前提醒一句如果只畫線不填充圖會顯得寡淡如果只填充不畫線讀者無法快速讀出準確高程。兩個一起上才能發(fā)揮第一章節(jié)里講的三層信息優(yōu)勢。4.2 等高線間距的選擇策略等高線間距是整張圖最容易出效果也最容易翻車的參數(shù)。間距過大地形細節(jié)丟失間距過小線條糊成一團。我通常不讓間距依賴手動拍腦袋而是用數(shù)據(jù)范圍自動計算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的倍數(shù)上邊界取整到20的倍數(shù)再按20米步長生成等高線列表。好處是輸出的等高線都是“20、40、60”這類整潔數(shù)字而不是“37.4、41.2”這種沒意義的帶小數(shù)級別。步長到底取多少我的經(jīng)驗是等高線數(shù)量控制在10到20條之間視覺效果最好。你可以這樣快速估算rough_step (vmax - vmin) / 15 step round(rough_step / 10) * 10 # 取整到10的倍數(shù)如果要畫陡峭山地等高線會自然堆積在陡坡處這時寧可選稍大的間距否則輸出成矢量圖后線條過多、文件尺寸也會失控。4.3 標注與線型的細節(jié)處理等高線標注重災區(qū)有兩個一是標注數(shù)字跟線交叉處被壓住看不清二是標注字體大小不匹配圖片縮放比例。clabel的inlineTrue是最值得開的參數(shù)它會在標注數(shù)字處把等高線切斷留出干凈的白色背景閱讀性提升非常明顯cs ax.contour(dem, levelslevels, colorsblack, linewidths0.7, alpha0.8, extentextent, originupper, zorder3) ax.clabel(cs, inlineTrue, fontsize8, fmt%d)fmt%d是第二個關(guān)鍵細節(jié)。默認標注會寫成浮點數(shù)可能出現(xiàn)50.0這種多余尾巴。用整數(shù)格式%d只要級別本身就是整數(shù)標注就會是干凈的50。如果你的高程有小數(shù)位可以改成fmt%.1f但配圖審美上最好統(tǒng)一。線型方面常規(guī)做法是主線用深色實線遇到一些特殊級別比如每5條線加深一次可以做成疊加方案。這個我在第5章的完整代碼里會演示。5. 圖層面板疊加順序、透明度與配色決定成敗5.1 圖層順序就是視覺的優(yōu)先級到這一步相當于做菜的最后組裝。圖層順序非常講究一般是從遠到近底圖陰影最先畫緊跟著是填充色最后是等高線線條和標注。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必須嚴格一致否則圖像上下翻轉(zhuǎn)或者坐標錯位第6章第一個坑就是這個?；叶汝幱白鳛榈讏D用imshow畫完后contourf和contour的坐標范圍必須與該extent完全對應(yīng)。5.2 配色方案默認terrain可以但你可以做得更好Matplotlib自帶terrain和gist_earth兩張地形色帶直接用的效果其實不算差但一個通病是低海拔區(qū)域顏色偏暗、和灰底圖疊在一起不夠清爽。我更推薦兩種方案一是安裝cmocean庫它的topo色帶就是為地形可視化設(shè)計的二是自制分段色帶完全控制顏色拐點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是個經(jīng)驗值太低底圖陰影蓋過顏色層次全無太高顏色蓋住灰度陰影立體感又沒了。0.5到0.7之間都可以建議輸出前多試兩個值。5.3 指北針、比例尺與色標的組合一張出版級地圖還需要地圖要素襯托。色標直接用colorbarcbar fig.colorbar(cs, axax, shrink0.8, labelElevation (m))比例尺優(yōu)先用matplotlib_scalebar這個第三方控件from matplotlib_scalebar.scalebar import ScaleBar ax.add_artist(ScaleBar(1, unitsm, length_fraction0.2))ScaleBar(1, unitsm)表示一個數(shù)據(jù)單位對應(yīng)1米如果你的DEM分辨率是10米這里第一個參數(shù)就填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)這段代碼把指北針固定在右上角位置用坐標軸比例表示不隨數(shù)據(jù)范圍變化。注意箭頭在文字下方先后順序別寫反。5.4 完整腳本串一遍把前面所有片段組合成一個可以獨立運行的最小腳本方便直接抄作業(yè)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)把這個腳本跑通一張合格的等高線陰影圖就出來了。接下來的內(nèi)容才是真正讓你少走彎路的部分。6. 實測踩坑記錄四個高頻坑與判斷依據(jù)6.1 坐標范圍不一致導致等高線“起飛”我見過最離譜的一次底圖陰影完全正常但等高線全部對不上位置有的線條甚至跑到圖外。最后定位到原因是imshow的extent寫的是[0, 120, 0, 100]而contour的extent寫成了[0, 100, 0, 120]行列順序顛倒了。排查方法很直接在一個已知的高峰位置打一個點看它是否同時落在底圖像元和等高線環(huán)上。如果點跟底圖對齊但跟等高線錯位問題幾乎一定在extent或origin。另外imshow默認originuppercontour默認數(shù)組坐標也是從上往下但一旦混入其他來源的數(shù)據(jù)就很容易一個上、一個下視覺表現(xiàn)就是等高線上下鏡像錯位。統(tǒng)一加originupper能解決大多數(shù)問題。6.2 陰影圖太暗或過亮陰影圖剛出來時很多人覺得山體太黑、細節(jié)被吞掉。大多數(shù)情況下不是算法錯了而是光照參數(shù)太激進。高度角越低影子越長45度是經(jīng)典默認值但對比較平緩的地形45度反而讓大片區(qū)域處于微光狀態(tài)。我的經(jīng)驗做法是先用altitude45看整體如果暗部比例超過三分之一就把高度角抬到55~60。另一個思路是把最終灰度做一次線性拉伸把原來0~255的區(qū)間壓縮到40~255相當于天然加了環(huán)境光hs hillshade(dem, res_x, res_y) hs 40 hs * (255 - 40) / 255這樣就保證最暗的陰影區(qū)也保留一定灰度不會出現(xiàn)死黑。環(huán)境光強度你可以按需調(diào)整40這個值不算拍腦袋是經(jīng)過多輪對比得到的平衡點太大會讓陰影失去層次太小又回到死黑。順便說一句如果你用clip(0,1)版本但不做拉伸出來的圖通常偏暗這是正?，F(xiàn)象不是代碼寫錯了。6.3 大范圍數(shù)據(jù)渲染卡頓真實DEM動輒幾千像素見方np.gradient其實非?？煺嬲钠款i在等高線算法和imshow渲染。整幅圖的內(nèi)存占用和輸出時間會隨著像素數(shù)線性增長到幾百萬像素時就明顯卡了。兩個辦法一是按顯示用途降采樣比如dem[::2, ::2]把行列各抽一半像素一下變成四分之一。降采樣后記得把res_x和res_y改成原來的兩倍否則坡度會翻倍。二是分塊繪制。科研繪圖常會遇到這種場景不需要全分辨率渲染只要清晰的等高線。我一般把DEM降到1000像素見方以下再出圖肉眼很難分辨細節(jié)差異但速度能快一個數(shù)量級。注意降采樣后extent不能變因為地理范圍沒變。6.4 中文字體亂碼Matplotlib默認字體對中文支持極差圖里只要出現(xiàn)中文標注就是方塊。這個問題說大不大但每回都能絆倒一些人plt.rcParams[font.sans-serif] [SimHei, Noto Sans CJK SC, WenQuanYi Zen Hei] plt.rcParams[axes.unicode_minus] False第一行設(shè)置中文字體按你系統(tǒng)里實際存在的字體選一個第二行是為了讓負號正常顯示不設(shè)置的話坐標軸負號有時會變成亂碼。如果是純英文標注這兩行可以完全省略。Linux服務(wù)器上渲染圖時建議先執(zhí)行fc-list :langzh看看系統(tǒng)里有沒有中文字體沒有就拿安裝包補一個這是服務(wù)器上最常見的坑。最后分享一個我自己的使用習慣做展示用的PNG圖把dpi設(shè)到300且保存時加bbox_inchestight避免出圖四周留白做論文配圖或后期還要進AI、CAD處理的建議直接存PDF矢量格式等高線和標注放大不糊。腳本化出圖的價值在于數(shù)據(jù)更新后我只需要改DEM路徑參數(shù)全部自動重算再也不用在GIS軟件里反復調(diào)一遍設(shè)置。這個工作流我用了很久算是從純手工畫圖到自動出圖之間性價比最高的一段路。