告别Cartopy！用Python Basemap + NOAA地形数据，5分钟搞定一张专业全球高程图

当Python地理绘图领域的新贵Cartopy逐渐成为主流时，许多开发者可能已经忘记了那个曾经陪伴我们多年的老朋友——Basemap。尽管官方早已宣布Basemap进入维护模式，但在某些特定场景下，这个"老将"依然展现出令人惊喜的实用价值。特别是当你需要快速绘制一张全球范围的高程图时，Basemap配合NOAA的ETOPO2v2数据源，能够以最简洁的代码实现专业级的地形可视化效果。

1. 为什么选择Basemap而非Cartopy？

在讨论技术实现之前，我们需要明确一个核心问题：在Cartopy已经成为事实标准的今天，为什么还要考虑使用Basemap？答案在于特定场景下的效率与简洁性。

Basemap的三大独特优势 ：

• 安装简单 ：仅需 conda install basemap 或 pip install basemap 即可完成，不像Cartopy可能遇到复杂的依赖问题
• API直观 ：对于简单的地图绘制，Basemap的API更加直接了当
• 数据兼容性好 ：特别是处理像ETOPO2v2这样的传统数据格式时

注意：Basemap虽然不再更新，但其核心功能依然稳定可靠，特别适合快速原型开发和教育用途。

相比之下，Cartopy虽然功能更强大、设计更现代，但在处理全球范围的基础地形图时，其复杂的投影系统和冗长的配置过程反而可能成为负担。

2. 数据准备：获取与理解ETOPO2v2地形数据

NOAA（美国国家海洋和大气管理局）提供的ETOPO2v2全球地形数据集是地理可视化领域的经典资源。这个数据集包含了全球陆地和海洋的地形信息，分辨率达到2弧分（约3.7公里）。

数据集关键特性 ：

属性	值
覆盖范围	全球（-180°到180°经度，-90°到90°纬度）
分辨率	2弧分
数据格式	netCDF
高程单位	米
数据维度	经度(x)、纬度(y)、高程(z)

下载数据只需访问NOAA官方网站：

wget https://www.ngdc.noaa.gov/mgg/global/relief/ETOPO2/ETOPO2v2-2006/ETOPO2v2c/netCDF/ETOPO2v2c_f4_netCDF.zip
unzip ETOPO2v2c_f4_netCDF.zip

使用Python读取这个netCDF文件非常简单，xarray库提供了直观的接口：

import xarray as xr
ds = xr.open_dataset('ETOPO2v2c_f4.nc')
print(ds)

3. Basemap绘图核心流程解析

现在让我们进入核心环节——使用Basemap绘制全球高程图。整个过程可以分为几个清晰的步骤，每个步骤都有其特定的配置选项和技巧。

3.1 基础地图设置

首先需要初始化地图投影。对于全球地图，圆柱投影（'cyl'）是最直接的选择：

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

plt.figure(figsize=(12, 6))
m = Basemap(projection='cyl', resolution='c',
            llcrnrlon=-180, llcrnrlat=-90,
            urcrnrlon=180, urcrnrlat=90)

关键参数说明 ：

• resolution ：控制海岸线精度（'c'粗糙，'l'低，'i'中，'h'高，'f'完整）
• llcrnrlon 等：设置地图边界（左下角经度/纬度，右上角经度/纬度）

3.2 高程数据可视化

ETOPO2v2数据的三维结构（经度、纬度、高程）非常适合用等高线填充图来表现：

import numpy as np

lon = np.linspace(-180, 180, len(ds['x']))
lat = np.linspace(-90, 90, len(ds['y']))
lon, lat = np.meshgrid(lon, lat)
dem = ds['z'].data

levels = [-8000, -4000, -2000, -1000, -500, -200, 0, 
          200, 500, 1000, 2000, 3000, 4000, 5000, 6000]
colors = ['#084594', '#2171b5', '#4292c6', '#6baed6', 
          '#9ecae1', '#c6dbef', '#006837', '#31a354',
          '#78c679', '#addd8e', '#d9f0a3', '#f7fcb9',
          '#c9bc87', '#a69165', '#856b49']

m.contourf(lon, lat, dem, levels=levels, colors=colors, extend='both')

3.3 地图修饰元素

专业的地图离不开精心设计的修饰元素。Basemap提供了丰富的方法来添加这些元素：

添加经纬网格 ：

m.drawmeridians(np.arange(-180, 181, 60), labels=[0,0,0,1], fontsize=8)
m.drawparallels(np.arange(-90, 91, 30), labels=[1,0,0,0], fontsize=8)

添加色标 ：

cb = m.colorbar(location='bottom', pad=0.3, size=0.05)
cb.set_label('Elevation (meters)', fontsize=10)
cb.ax.tick_params(labelsize=8)

4. 完整代码实现与优化技巧

将上述步骤整合起来，我们得到完整的绘图代码。这段代码不仅实现了基本功能，还包含了一些实用优化：

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

# 数据准备
ds = xr.open_dataset('ETOPO2v2c_f4.nc')
lon = np.linspace(-180, 180, len(ds['x']))
lat = np.linspace(-90, 90, len(ds['y']))
lon, lat = np.meshgrid(lon, lat)
dem = ds['z'].data

# 绘图设置
plt.figure(figsize=(16, 8), dpi=300)
plt.rcParams['font.family'] = 'Arial'

# 创建地图
m = Basemap(projection='cyl', resolution='i',
            llcrnrlon=-180, llcrnrlat=-90,
            urcrnrlon=180, urcrnrlat=90)

# 绘制地形
levels = [-8000, -6000, -4000, -2000, -1000, -500, 
          -200, 0, 200, 500, 1000, 1500, 2000, 
          3000, 4000, 5000, 6000]
colors = ['#084594', '#2171b5', '#4292c6', '#6baed6',
          '#9ecae1', '#c6dbef', '#deebf7', '#006837',
          '#31a354', '#78c679', '#addd8e', '#d9f0a3',
          '#f7fcb9', '#c9bc87', '#a69165', '#856b49',
          '#664830']
cs = m.contourf(lon, lat, dem, levels=levels, 
                colors=colors, extend='both')

# 添加地图元素
m.drawcoastlines(linewidth=0.5, color='gray')
m.drawcountries(linewidth=0.25, color='gray')
m.drawmeridians(np.arange(-180, 181, 60), 
                labels=[0,0,0,1], fontsize=8, linewidth=0.5)
m.drawparallels(np.arange(-90, 91, 30), 
                labels=[1,0,0,0], fontsize=8, linewidth=0.5)

# 添加色标
cb = m.colorbar(cs, location='bottom', pad=0.2, size=0.05)
cb.set_label('Elevation (meters)', fontsize=10, labelpad=10)
cb.ax.tick_params(labelsize=8)

# 保存和显示
plt.title('Global Elevation Map (ETOPO2v2)', fontsize=12, pad=20)
plt.savefig('global_elevation.png', dpi=300, bbox_inches='tight')
plt.show()

性能优化技巧 ：

• 降低数据分辨率：对于快速预览，可以跳过部分数据点

dem = ds['z'].data[::2, ::2]  # 每隔一个点取一个

• 使用更简单的投影：如墨卡托投影（'merc'）在某些情况下渲染更快
• 调整分辨率参数： resolution='l' （低）可以显著提高绘制速度

5. 进阶应用：从静态图到动态探索

基础地形图只是起点，Basemap还能支持更丰富的可视化形式。以下是几个值得尝试的进阶方向：

三维地形渲染 ：

from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')

# 降低数据密度以提高性能
dem_small = dem[::10, ::10]
lon_small = lon[::10, ::10]
lat_small = lat[::10, ::10]

ax.plot_surface(lon_small, lat_small, dem_small, 
                cmap='terrain', rstride=1, cstride=1)
ax.set_zlim(-8000, 8000)

区域聚焦分析 ：

# 聚焦在特定区域（如喜马拉雅山脉）
m = Basemap(projection='merc', resolution='h',
            llcrnrlon=70, llcrnrlat=20,
            urcrnrlon=110, urcrnrlat=40)
m.contourf(lon, lat, dem, levels=levels, colors=colors)
m.drawcountries()
m.drawrivers(color='blue')

交互式探索 ： 虽然Basemap本身不支持交互，但可以结合matplotlib的交互模式实现简单的数据探查：

plt.ion()  # 开启交互模式
fig = plt.figure()
# ...（绘图代码）
def onclick(event):
    if event.inaxes:
        print(f"Clicked at: {event.xdata}, {event.ydata}")
fig.canvas.mpl_connect('button_press_event', onclick)

在实际科研工作中，我发现Basemap特别适合快速验证数据质量和制作初步分析图表。当需要更复杂的空间分析时，再转向Cartopy或其他专业GIS工具也不迟。这种"先用Basemap快速验证，再用其他工具深入分析"的工作流程，往往能节省大量时间。