眾所周知，Python的matplotlib是一個非常全面的制圖庫，它不僅可以繪制圖表、地圖，還可以繪制3D效果圖，試想一下，如果你在畫圖的時候，可以將立體地形圖作為底圖，那逼格噌一下子就上來了，今天我就來教大家畫一個立體地形圖，啥也不說，咱先上效果圖：
上面這張圖是展示了基于matplotlib+cartopy的山地陰影圖在不同光影參數(shù)下的變化效果。這個變化效果有利于我們理解matplotlib對該效果的設(shè)計(jì)理念。
在我講解之前，我推薦大家讀一下matplotlib官方文檔庫里的這一篇文章：《Topographic hillshading》，該文章已經(jīng)介紹了如何單獨(dú)基于matplotlib繪制山地陰影圖，并給出了不同渲染參數(shù)下的渲染效果圖。我當(dāng)初對山地立體圖的學(xué)習(xí)就是從這篇文章開始的。
本教程代碼所需依賴：
matplotlib
cartopy>=0.19.0
cnmaps
netCDF4
numpy
本教程使用的DEM數(shù)據(jù)：原始數(shù)據(jù)來自公開數(shù)據(jù)集ASTER DEM，已處理成中國區(qū)的NetCDF格式。
另外下文代碼中會出現(xiàn)cnmaps這個新寫的包，如果你對這個包較陌生想要了解這個包的使用方法的請移步我的往期文章：如何用Python優(yōu)雅地繪制中國的地圖
神說：要有光
光，是三維世界最重要的東西，要繪制山地立體圖，首先需要理解matplotlib中的LightSource對象，顧名思義，這個對象就是“光源”，與3D 建模里的光源是同一個東西，它的調(diào)用方法是：
from matplotlib.colors import LightSource
ls = LightSource(azdeg=360, altdeg=30)
其中azdeg是方位角，取值范圍是0~360，altdeg是高度角，取值范圍是0~90，這兩個參數(shù)可以確定一個光源的投射方向，進(jìn)而可以知道被光源投射的物體，哪一部分應(yīng)該是光，哪一部分應(yīng)該是影，而光影便是打開地形立體效果的鑰匙。
在我們創(chuàng)建了光源以后，就需要基于該光源對地形數(shù)據(jù)生成光影對象，通常情況下，對于山地陰影，我們有兩個方法可以選擇，一個是hillshade，另一個是shade，其中hillshade返回的是以0-1的數(shù)字代表的光影明暗特征，你可以把它理解為一個灰度圖，而shade返回的是一個RGBA數(shù)組，也就是彩圖，下面我們使用shade來看一個實(shí)際的例子：
import netCDF4 as nc
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from matplotlib.colors import LightSource
from cnmaps import get_map, draw_map
ds = nc.Dataset('./data/cldasgrid_dem.nc')
_lon = ds.variables['LON'][:]
_lat = ds.variables['LAT'][:]
_dem = ds.variables['elevation'][:]
lon = _lon[4032: 4662]
lat = _lat[1635: 2134]
dem = _dem[1635: 2134, 4032: 4662]
ls = LightSource(azdeg=360, altdeg=30)
rgb = ls.shade(dem[::-1], cmap=plt.cm.gist_earth, blend_mode='overlay',
vert_exag=0.5, dx=10, dy=10, fraction=1.5, vmin=-2300)
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
img = ax.imshow(rgb, extent=(lon.min(),lon.max(),lat.min(),lat.max()), transform=ccrs.PlateCarree())
draw_map(get_map('河南'), color='w', linewidth=2)
ax.set_extent(get_map('河南').get_extent(buffer=0))
plt.show()
這樣，我們的第一張立體地形圖就出來了，如果你調(diào)整azdeg和altdeg的值，陰影的方位就會隨之改變，就像文章開頭那張動圖一樣，它就是通過以10為間隔修改azdeg的值以達(dá)到光線旋轉(zhuǎn)照射的效果的。
光影參數(shù)詳解
接下來，我們需要了解一下ls.shade方法的各個參數(shù)是干什么的，首先第一個位置函數(shù)肯定是我們的dem數(shù)據(jù)，這里需要注意的是，你必須把dem的緯度順序調(diào)整為低緯->高緯的順序，否則渲染出來的圖片是反的。cmap是色標(biāo)這個大家應(yīng)該都知道就不贅述了，你可以使用matplotlib中預(yù)置的任何你喜歡的色標(biāo)，blend_mode這個參數(shù)大家會比較陌生，它是一種渲染模式選擇，預(yù)置選項(xiàng)有：'hsv'，'overlay'，'soft'。它們分別是什么意思呢？官方文檔在這個參數(shù)上的解釋一大堆，但是根據(jù)我的理解，它們就類似于各種“濾鏡”，你可以調(diào)整這三個濾鏡來看看哪個是你喜歡的效果，我們來試一下：
import netCDF4 as nc
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from matplotlib.colors import LightSource
from cnmaps import get_map, draw_map
ds = nc.Dataset('./data/cldasgrid_dem.nc')
_lon = ds.variables['LON'][:]
_lat = ds.variables['LAT'][:]
_dem = ds.variables['elevation'][:]
lon = _lon[4032: 4662]
lat = _lat[1635: 2134]
dem = _dem[1635: 2134, 4032: 4662]
ls = LightSource(azdeg=360, altdeg=30)
fig = plt.figure(figsize=(8*3, 8))
fig.tight_layout()
for i, mode in enumerate(['soft', 'overlay', 'hsv']):
rgb = ls.shade(dem[::-1], cmap=plt.cm.gist_earth, blend_mode=mode,
vert_exag=0.5, dx=10, dy=10, fraction=1.5, vmin=-2300)
ax = fig.add_subplot(1,3,i+1, projection=ccrs.PlateCarree())
img = ax.imshow(rgb, extent=(lon.min(),lon.max(),lat.min(),lat.max()), transform=ccrs.PlateCarree())
draw_map(get_map('河南'), color='w', linewidth=2)
plt.title(mode)
ax.set_extent(get_map('河南').get_extent(buffer=0))
plt.show()
看得出來，還是中間的overlay看起來更均衡一些。
下面我們來看一下vert_exag參數(shù)，這個參數(shù)表征的是頂點(diǎn)的突出程度，這個點(diǎn)的值越大，立體感會越強(qiáng)，但是如果你把它設(shè)得太大，也會有一些副作用，來看示例：
import netCDF4 as nc
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from matplotlib.colors import LightSource
from cnmaps import get_map, draw_map
ds = nc.Dataset('./data/cldasgrid_dem.nc')
_lon = ds.variables['LON'][:]
_lat = ds.variables['LAT'][:]
_dem = ds.variables['elevation'][:]
lon = _lon[4032: 4662]
lat = _lat[1635: 2134]
dem = _dem[1635: 2134, 4032: 4662]
ls = LightSource(azdeg=360, altdeg=30)
fig = plt.figure(figsize=(8*3, 8))
fig.tight_layout()
for i, ve in enumerate([0.5, 1, 10]):
rgb = ls.shade(dem[::-1], cmap=plt.cm.gist_earth, blend_mode='overlay',
vert_exag=ve, dx=10, dy=10, fraction=1.5, vmin=-2300)
ax = fig.add_subplot(1,3,i+1, projection=ccrs.PlateCarree())
img = ax.imshow(rgb, extent=(lon.min(),lon.max(),lat.min(),lat.max()), transform=ccrs.PlateCarree())
draw_map(get_map('河南'), color='w', linewidth=2)
plt.title(ve)
ax.set_extent(get_map('河南').get_extent(buffer=0))
plt.show()
上圖從左到右vert_exag分別等于0.5、1和10，可以看出來，當(dāng)vert_exag==10的時候，平原地圖會有很強(qiáng)的噪點(diǎn)效果，這也是vert_exag值設(shè)置過大的一個副作用，在實(shí)際繪圖的時候，需要綜合考慮，選一個合適的值。當(dāng)然，對于vert_exag參數(shù)，還有另外兩個參數(shù)會與之配合（或者說制衡），那就是dx和dy，這兩個參數(shù)的含義是在平面空間上單個頂點(diǎn)的重采樣間隔，dx和dy的值越小，圖像越能展現(xiàn)原始的數(shù)據(jù)細(xì)節(jié)，dx和dy的值越大，那么最終出來的圖越平滑，一些原始數(shù)據(jù)的細(xì)節(jié)就會被平滑掉了。下面我們來看一下不同的dx，dy取值，對圖像效果有什么影響。
import netCDF4 as nc
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from matplotlib.colors import LightSource
from cnmaps import get_map, draw_map
ds = nc.Dataset('./data/cldasgrid_dem.nc')
_lon = ds.variables['LON'][:]
_lat = ds.variables['LAT'][:]
_dem = ds.variables['elevation'][:]
lon = _lon[4032: 4662]
lat = _lat[1635: 2134]
dem = _dem[1635: 2134, 4032: 4662]
ls = LightSource(azdeg=360, altdeg=30)
fig = plt.figure(figsize=(8*3, 8))
fig.tight_layout()
for i, dd in enumerate([1, 10, 100]):
rgb = ls.shade(dem[::-1], cmap=plt.cm.gist_earth, blend_mode='overlay',
vert_exag=0.5, dx=dd, dy=dd, fraction=1.5, vmin=-2300)
ax = fig.add_subplot(1,3,i+1, projection=ccrs.PlateCarree())
img = ax.imshow(rgb, extent=(lon.min(),lon.max(),lat.min(),lat.max()), transform=ccrs.PlateCarree())
draw_map(get_map('河南'), color='w', linewidth=2)
plt.title(f'dx={dd}, dy={dd}')
ax.set_extent(get_map('河南').get_extent(buffer=0))
plt.show()
可以看到，當(dāng)dx和dy偏小的時候，同樣出現(xiàn)了噪點(diǎn)問題。
最后還有一個很重要的參數(shù)就是fraction，它是一個控制光影效果強(qiáng)度的參數(shù)，這個值越大，明暗的對比度就越大，我們來看一下對比效果：
import netCDF4 as nc
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from matplotlib.colors import LightSource
from cnmaps import get_map, draw_map
ds = nc.Dataset('./data/cldasgrid_dem.nc')
_lon = ds.variables['LON'][:]
_lat = ds.variables['LAT'][:]
_dem = ds.variables['elevation'][:]
lon = _lon[4032: 4662]
lat = _lat[1635: 2134]
dem = _dem[1635: 2134, 4032: 4662]
ls = LightSource(azdeg=360, altdeg=30)
fig = plt.figure(figsize=(8*3, 8))
fig.tight_layout()
for i, fraction in enumerate([0.1, 1, 2]):
rgb = ls.shade(dem[::-1], cmap=plt.cm.gist_earth, blend_mode='overlay',
vert_exag=0.5, dx=10, dy=10, fraction=fraction, vmin=-2300)
ax = fig.add_subplot(1,3,i+1, projection=ccrs.PlateCarree())
img = ax.imshow(rgb, extent=(lon.min(),lon.max(),lat.min(),lat.max()), transform=ccrs.PlateCarree())
draw_map(get_map('河南'), color='w', linewidth=2)
plt.title(f'fraction={fraction}')
ax.set_extent(get_map('河南').get_extent(buffer=0))
plt.show()
可以看到當(dāng)fraction=0.1時，幾乎沒有陰影效果，而fraction=2的時候，陰影效果很強(qiáng)烈，甚至感覺有點(diǎn)耀眼。
上述的山地陰影圖，不僅可以自嗨，還可以與你的其他數(shù)據(jù)結(jié)合起來，一起組成一個多圖層的效果圖，例如：
上圖展示了2021年7月20日鄭州特大暴雨的逐小時降水量在一天中的分布變化，降水?dāng)?shù)據(jù)源是中國氣象局的CMPAS格點(diǎn)降水產(chǎn)品，由于該數(shù)據(jù)非公開數(shù)據(jù)集，在此不提供數(shù)據(jù)下載。上圖的繪制方法就是在前面代碼的基礎(chǔ)上，增加了ax.countourf函數(shù)對降水?dāng)?shù)據(jù)的疊加，在這里就不再贅述。
個人網(wǎng)站：www.clarmy.net
GitHub：https://github.com/Clarmy
歡迎轉(zhuǎn)發(fā)，如有轉(zhuǎn)發(fā)需要，請后臺聯(lián)系我添加長白。
