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