IT / 地理/GIS / 气象/Meteorology · 2022年2月19日 7

如何用Python绘制炫酷的立体地形图

众所周知,Python的matplotlib是一个非常全面的制图库,它不仅可以绘制图表、地图,还可以绘制3D效果图,试想一下,如果你在写论文的时候,可以将立体地形图作为底图,那逼格噌一下子就上来了,今天我就来教大家画一个带有地理坐标属性的立体地形图,啥也不说,咱先上效果图:

上面这张图是展示了基于matplotlib+cartopy的山地阴影图在不同光影参数下的变化效果。这个变化效果有利于我们理解matplotlib对该效果的设计理念。

在我讲解之前,我推荐大家读一下matplotlib官方文档库里的这一篇文章:Topographic hillshading,该文章已经介绍了如何单独基于matplotlib绘制山地阴影图,并给出了不同渲染参数下的渲染效果图。我当初对山地立体图的学习就是从这篇文章开始的。

本教程代码所需依赖:

matplotlib
cartopy>=0.19.0
cnmaps==0.2.1
netCDF4
numpy

本教程使用的DEM数据:和鲸社区数据集: CHINA_DEM 

另外如果想要了解cnmaps这个包的请移步往期文章:如何用Python优雅地绘制中国的地图

神说:要有光

光,是三维世界最重要的东西,要绘制山地立体图,首先需要理解matplotlib中的LightSource对象,顾名思义,这个对象就是“光源”,与3D 建模里的光源是同一个东西,它的调用方法是:

from matplotlib.colors import LightSource
ls = LightSource(azdeg=360, altdeg=30)

其中azdeg是方位角,altdeg是高度角,这两个参数可以确定一个光源的投射方向,进而可以知道被光源投射的物体,哪一部分应该是光,哪一部分应该是影,而光影便是实现地形立体效果的金钥匙。

在我们创建了光源以后,就需要基于该光源对地形数据生成光影对象,通常情况下,对于山地阴影,我们有两个方法可以选择,一个是hillshade,另一个是shade,其中hillshade返回的是以0-1的数字代表的光影明暗特征,你可以把它理解为一个灰度图,而shade返回的是一个RGBA数组,也就是彩图,下面我们使用shade来看一个实际的例子:

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()

这样,我们的第一张立体地形图就出来了,是不是很炫酷?此外,如果你调整azdegaltdeg的值,阴影的方位就会随之改变,就像文章开头那张动图一样,它就是通过修改azdeg的值以达到光线旋转照射的效果的。

光影参数详解

接下来,我们需要了解一下ls.shade方法的各个参数是干什么的,首先第一个位置函数肯定是我们的dem数据,这里需要注意的是,你必须把dem的纬度顺序调整为低纬->高纬的顺序,否则渲染出来的图片是反的。cmap是色标这个大家应该都知道就不赘述了,blend_mode这个参数大家会比较陌生,它是一种渲染模式选择,预置选项有:'hsv''overlay''soft'。它们分别是什么意思呢?官方文档在这个参数上的解释一大堆,但是根据我的理解,这个参数其实就类似于“滤镜”,你可以调整这三个滤镜来看看哪个是你喜欢的效果,我们来试一下:

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参数,这个参数表征的是顶点的突出程度,这个点的值越大,立体感会越强,但是如果你把它设得太大,也会有一些副作用,来看示例:

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,可以看出来,当vert_exag==10的时候,平原地图会有很强的噪点效果,这也是vert_exag值设置过大的一个副作用,在实际绘图的时候,需要综合考虑,选一个合适的值。当然,对于vert_exag参数,还有另外两个参数会与之配合(或者说制衡),那就是dxdy,这两个参数的含义是在平面空间上单个顶点的重采样间隔,dxdy的值越小,图像越能展现原始的数据细节,dxdy的值越大,那么最终出来的图越平滑,一些原始数据的细节就会被平滑掉了。下面我们来看一下不同的dxdy取值,对图像效果有什么影响。

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()

可以看到,当dxdy偏小的时候,同样出现了噪点问题。

最后还有一个很重要的参数就是fraction,它是一个控制光影效果强度的参数,这个值越大,明暗的对比度就越大,我们来看一下对比效果:

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()

可以看到当fraction=0.1时,几乎没有阴影效果,而fraction=2的时候,阴影效果很强烈,甚至感觉有点耀眼。

最后,让我们在一个动图效果下结束我们今天的讲解:

上图展示了2021年7月20日郑州特大暴雨的逐小时降水量在一天中的分布变化,降水数据源是中国气象局的CMPAS格点降水产品,由于该数据非公开数据集,在此不提供数据下载。