如何读取netcdf数据并在matplotlib Basemap上绘图

发布时间：2021-12-05 14:47:55 作者：小新
来源：亿速云阅读：336

# 如何读取NetCDF数据并在Matplotlib Basemap上绘图

## 1. 引言

NetCDF（Network Common Data Form）是一种广泛应用于科学数据存储的格式，特别适合存储多维数组数据（如气候、海洋和大气数据）。Matplotlib是Python中最流行的绘图库之一，而Basemap是其扩展工具包，专门用于地理空间数据可视化。本文将详细介绍如何读取NetCDF文件，并使用Matplotlib Basemap进行地理空间数据可视化。

---

## 2. 准备工作

### 2.1 安装必要库
首先需要安装以下Python库：
```bash
pip install netCDF4 matplotlib basemap numpy

2.2 获取示例数据

推荐使用公开的气候数据集（如NOAA或CMIP6模型数据）。本文以air.sig995.2012.nc（NOAA地表气温数据）为例。

3. 读取NetCDF数据

3.1 导入库

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

3.2 加载NetCDF文件

# 打开NetCDF文件
dataset = nc.Dataset('air.sig995.2012.nc')

# 查看文件结构
print(dataset)
print(dataset.variables.keys())

# 提取变量
lats = dataset.variables['lat'][:]  # 纬度
lons = dataset.variables['lon'][:]  # 经度
time = dataset.variables['time'][:] # 时间
air = dataset.variables['air'][:]   # 气温数据（假设为3D数组：time x lat x lon）

3.3 数据预处理

# 选择第一个时间点的数据
air_2d = air[0, :, :]

# 处理缺失值（如有）
air_2d = np.ma.masked_where(air_2d > 1e10, air_2d)

4. 使用Basemap绘制地图

4.1 创建Basemap实例

m = Basemap(
    projection='merc',            # 墨卡托投影
    llcrnrlat=lats.min(),        # 左下角纬度
    urcrnrlat=lats.max(),        # 右上角纬度
    llcrnrlon=lons.min(),        # 左下角经度
    urcrnrlon=lons.max(),        # 右上角经度
    resolution='i'               # 中等分辨率
)

4.2 绘制基础地图元素

fig = plt.figure(figsize=(12, 8))
m.drawcoastlines(linewidth=0.5)
m.drawcountries(linewidth=0.5)
m.drawstates(linewidth=0.2)
m.fillcontinents(color='lightgray', lake_color='aqua')
m.drawparallels(np.arange(-90, 90, 30), labels=[1,0,0,0])
m.drawmeridians(np.arange(-180, 180, 60), labels=[0,0,0,1])

4.3 转换坐标网格

# 创建经纬度网格
lon_grid, lat_grid = np.meshgrid(lons, lats)

# 将经纬度转换为地图投影坐标
x, y = m(lon_grid, lat_grid)

4.4 绘制数据

# 绘制等值线填充图
cs = m.contourf(x, y, air_2d, levels=20, cmap=plt.cm.jet)

# 添加颜色条
cbar = m.colorbar(cs, location='bottom', pad="10%")
cbar.set_label('Temperature (K)')

# 添加标题
plt.title('Surface Air Temperature (2012-01-01)')

4.5 完整代码示例

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

# 读取数据
dataset = nc.Dataset('air.sig995.2012.nc')
lats = dataset.variables['lat'][:]
lons = dataset.variables['lon'][:]
air = dataset.variables['air'][0, :, :]

# 创建地图
m = Basemap(projection='merc', 
            llcrnrlat=lats.min(), 
            urcrnrlat=lats.max(),
            llcrnrlon=lons.min(),
            urcrnrlon=lons.max(),
            resolution='i')

fig = plt.figure(figsize=(12, 8))
m.drawcoastlines()
m.drawcountries()
m.fillcontinents(color='lightgray')

# 绘制数据
lon_grid, lat_grid = np.meshgrid(lons, lats)
x, y = m(lon_grid, lat_grid)
cs = m.contourf(x, y, air, levels=20, cmap=plt.cm.jet)
cbar = m.colorbar(cs, location='bottom', pad="10%")
plt.title('Global Surface Temperature')
plt.show()

5. 进阶技巧

5.1 处理时间维度

from netCDF4 import num2date

# 转换时间坐标
times = dataset.variables['time']
dates = num2date(times[:], units=times.units)

# 绘制时间序列动画
for i in range(len(dates)):
    air_slice = air[i, :, :]
    # 更新绘图代码...

5.2 区域地图设置

# 聚焦中国区域
m = Basemap(llcrnrlon=70, llcrnrlat=15,
            urcrnrlon=140, urcrnrlat=55,
            projection='lcc', lat_1=33, lat_2=45, lon_0=100)

5.3 添加矢量数据

# 绘制风场箭头
u = dataset.variables['uwnd'][0, :, :]
v = dataset.variables['vwnd'][0, :, :]
m.quiver(x[::5,::5], y[::5,::5], u[::5,::5], v[::5,::5], scale=500)

6. 常见问题解决

投影变形问题：
- 尝试不同的投影方式（如'ortho'、'robin'）

内存不足：

使用xarray库进行分块读取

import xarray as xr
ds = xr.open_dataset('large_file.nc', chunks={'time': 10})

坐标不匹配：
- 检查数据是否使用-180~180或0~360经度范围
```
lons = np.where(lons > 180, lons - 360, lons)
```

7. 结论

通过结合NetCDF4和Matplotlib Basemap，我们可以高效地处理和可视化地理空间数据。本文展示了从数据读取到地图绘制的完整流程，并提供了进阶技巧和问题解决方案。对于更复杂的应用，建议进一步学习cartopy（Basemap的替代库）和xarray等工具。

注意：Basemap已停止维护，官方推荐迁移至cartopy。但Basemap仍广泛用于现有项目中。 “`