Python绘图：使用Basemap绘制全球地形图

目的

• 基于NOAA-NGDC的全球地形数据ETOPO2v2，绘制全球地形图。

思路

1. 读取全球地形数据
2. 绘制全球地形图

数据来源

• ETOPO2v2c_f4.nc下载地址：https://www.ngdc.noaa.gov/mgg/global/relief/ETOPO2/ETOPO2v2-2006/ETOPO2v2c/netCDF/ETOPO2v2c_f4_netCDF.zip

ETOPO2v2是美国国家海洋和大气管理局(NOAA)下属的国家地球物理数据中心(NGDC)开发的全球地形模型，包括全球陆地和海洋的地形，分辨率为2分。

步骤

1. 读取全球地形数据

数据格式为netCDF(.nc)，可以考虑使用xarray库来读取nc文件。

import xarray as xr    # 导入xarray库
ds = xr.open_dataset('ETOPO2v2c_f4.nc')    # 读取全球地形数据
ds   # 显示nc文件信息

• 读取ETOPO2v2c_f4.nc文件信息可知，x表示经度，范围为-180～180度；y表示纬度，范围为-90～90度；z表示海拔。

2. 绘制全球地形图

地形图的绘制使用Basemap包。Basemap是Python可视化库Matplotlib下的一个工具包，主要功能是绘制二维地图，是常用的地理数据可视化工具。尽管Basemap逐渐被Cartopy所取代，但个人认为某些地方Basemap使用起来比Cartopy更加方便好用。

• 导入绘图所需的库

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

• 准备绘图数据

# 准备用于绘图的数据
lon = np.linspace(min(ds['x'].data), max(ds['x'].data), len(ds['x'].data))  # 经度
lat = np.linspace(min(ds['y'].data), max(ds['y'].data), len(ds['y'].data))  # 纬度
lon, lat = np.meshgrid(lon, lat)  # 构建经纬网
dem = ds['z'].data  # DEM数据

• 绘图

    # 设置地图全局属性
    plt.rcParams['font.sans-serif'] = ['Times New Roman']  # 设置整体的字体为Times New Roman
    plt.figure(figsize=(10, 6), dpi=600)  # 设置大小和分辨率

    # 创建底图，设置地图投影为World Plate Carrée，分辨率为高分辨率，地图范围为全球
    m = Basemap(projection='cyl', resolution='h', llcrnrlon=-180, llcrnrlat=-90, urcrnrlon=180, urcrnrlat=90)

    # 设置地图经纬线，并只在左端和底端显示
    m.drawmeridians(np.arange(-180, 181, 30), labels=[0, 0, 0, 1], fontsize=10, linewidth=0.8, color='silver')  # 经线
    m.drawparallels(np.arange(-90, 91, 30), labels=[1, 0, 0, 0], fontsize=10, linewidth=0.8, color='silver')  # 纬线

    # 绘制地图
    levels = [-8000, -6000, -4000, -2000, -1000, -200, -50, 0, 50, 200, 500, 1000, 1500, 2000, 3000, 4000, 
              5000, 6000, 7000, 8000]  # 创建分级
    color = ['#084594', '#2171b5', '#4292c6', '#6baed6', '#9ecae1', '#c6dbef', '#deebf7', '#006837', 
             '#31a354', '#78c679', '#addd8e', '#d9f0a3', '#f7fcb9', '#c9bc87', '#a69165', '#856b49', 
             '#664830', '#ad9591', '#d7ccca']  # 设置色带
    m.contourf(lon, lat, dem, levels=levels, extend='both', colors=color)  # 绘图，并设置图例两端显示尖端

    # 设置图例
    cb = m.colorbar(location='bottom', pad=0.35)  # 图例在底端显示
    cb.set_ticks(levels)  # 设置色带刻度
    cb.ax.tick_params(labelsize=10)  # 刻度字号大小
    cb.set_label('Global Elevation (meter)', fontsize=12)  # 设置图例名称和字体大小

    # 保存图片并显示
    plt.savefig('global_elevation.jpg', dpi=600, bbox_inches='tight', pad_inches=0.1)  # 输出地图，并设置边框空白紧密
    plt.show()  # 显示地图

• 绘制结果

全部代码

# -*- encoding: utf-8 -*-
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

# 主函数
if __name__ == '__main__':
    # 读取全球地形数据
    ds = xr.open_dataset('ETOPO2v2c_f4.nc')

    # 准备用于绘图的数据
    lon = np.linspace(min(ds['x'].data), max(ds['x'].data), len(ds['x'].data))  # 经度
    lat = np.linspace(min(ds['y'].data), max(ds['y'].data), len(ds['y'].data))  # 纬度
    lon, lat = np.meshgrid(lon, lat)  # 构建经纬网
    dem = ds['z'].data  # DEM数据

    # 设置地图全局属性
    plt.rcParams['font.sans-serif'] = ['Times New Roman']  # 设置整体的字体为Times New Roman
    plt.figure(figsize=(10, 6), dpi=600)  # 设置大小和分辨率

    # 创建底图，设置地图投影为World Plate Carrée，分辨率为高分辨率，地图范围为全球
    m = Basemap(projection='cyl', resolution='h', llcrnrlon=-180, llcrnrlat=-90, urcrnrlon=180, urcrnrlat=90)

    # 设置地图经纬线，并只在左端和底端显示
    m.drawmeridians(np.arange(-180, 181, 30), labels=[0, 0, 0, 1], fontsize=10, linewidth=0.8, color='silver')  # 经线
    m.drawparallels(np.arange(-90, 91, 30), labels=[1, 0, 0, 0], fontsize=10, linewidth=0.8, color='silver')  # 纬线

    # 绘制地图
    levels = [-8000, -6000, -4000, -2000, -1000, -200, -50, 0, 50, 200, 500, 1000, 1500, 2000, 3000, 4000, 
              5000, 6000, 7000, 8000]  # 创建分级
    color = ['#084594', '#2171b5', '#4292c6', '#6baed6', '#9ecae1', '#c6dbef', '#deebf7', '#006837', 
             '#31a354', '#78c679', '#addd8e', '#d9f0a3', '#f7fcb9', '#c9bc87', '#a69165', '#856b49', 
             '#664830', '#ad9591', '#d7ccca']  # 设置色带
    m.contourf(lon, lat, dem, levels=levels, extend='both', colors=color)  # 绘图，并设置图例两端显示尖端

    # 设置图例
    cb = m.colorbar(location='bottom', pad=0.35)  # 图例在底端显示
    cb.set_ticks(levels)  # 设置色带刻度
    cb.ax.tick_params(labelsize=10)  # 刻度字号大小
    cb.set_label('Global Elevation (meter)', fontsize=12)  # 设置图例名称和字体大小

    # 保存图片并显示
    plt.savefig('global_elevation.jpg', dpi=600, bbox_inches='tight', pad_inches=0.1)  # 输出地图，并设置边框空白紧密
    plt.show()  # 显示地图

最后

• 内容仅供大家学习参考，若有不足之处，敬请大家批评指正！

参考资料

1. xarray官方文档：http://xarray.pydata.org/en/stable/index.html
2. Basemap官方文档：https://basemaptutorial.readthedocs.io/en/latest/index.html
3. Matplotlib官方文档https://matplotlib.org/3.2.0/index.html
4. 色带设置：https://colorbrewer2.org/