Python气象数据可视化：用Cartopy+frykit绘制中国AQI分布图（附完整代码）

Python气象数据可视化实战：从散点到专业地图的完整工作流

最近在分析全国空气质量数据时，我遇到了一个典型问题：如何将338个城市的AQI监测点数据，转化为一张既专业又美观的全国分布图？这不仅仅是简单的散点图绘制，而是涉及地理坐标转换、空间插值、地图投影、边界裁剪等一系列复杂操作。经过几周的摸索和实践，我总结出了一套基于Cartopy和frykit的高效解决方案，今天就来分享这个完整的工作流程。

如果你也经常需要处理气象、环境或地理空间数据，这篇文章将为你提供一个可以直接复用的代码框架。不同于简单的教程，我会深入讲解每个环节的技术细节和实际应用中可能遇到的坑，特别是如何处理中国地图的特殊需求。

1. 环境配置与工具选择：为什么是Cartopy+frykit组合？

在Python地理可视化领域，有几个主流选择：Basemap（已停止维护）、GeoPandas、Plotly，以及我们今天要重点介绍的Cartopy。Cartopy由英国气象局开发，专门为科学数据可视化设计，与Matplotlib无缝集成，支持多种地图投影和地理特征。

但Cartopy有个“小问题”——它自带的全球地图数据在某些区域边界上不够精确。这就是为什么我们需要frykit这个国产工具包。frykit专门针对中国地图优化，提供了准确的国界、省界数据，以及九段线等必要元素，还能大幅提升绘制速度。

1.1 环境搭建：一步到位的配置方案

我强烈推荐使用Conda进行环境管理，特别是处理地理空间数据时，依赖库的版本兼容性是个大问题。下面是我的标准配置流程：

# 创建专用环境
conda create -n geo_viz python=3.9 -y
conda activate geo_viz

# 核心地理可视化库
conda install -c conda-forge cartopy -y

# 数据处理与分析基础包
conda install numpy pandas xarray scipy -y

# 插值算法库（克里金插值）
conda install -c conda-forge pykrige -y

# 中国地图专用工具包
pip install frykit

# 可选但推荐的辅助工具
pip install openpyxl  # 处理Excel数据
pip install cnmaps    # 另一个中国地图工具，作为备选

注意：如果你之前尝试过从源码编译安装Cartopy，可能会遇到Proj、GEOS等C++库的依赖问题。Conda-forge渠道已经帮你解决了所有底层依赖，这是最省心的方式。

1.2 工具对比：为什么选择这个技术栈？

为了让你更清楚这个技术组合的优势，我整理了一个对比表格：

这个组合的核心思路是：用Cartopy处理地理投影和基础地图，用frykit补充中国特定的地理要素，用PyKrige进行空间插值，形成一个完整的工作流。

2. 数据准备与预处理：从原始数据到可插值格式

在实际项目中，数据往往来自不同的监测站、不同的格式。我最近处理的一个典型数据集包含338个城市的AQI（空气质量指数）小时数据，存储格式是Excel，结构如下：

2.1 数据读取与质量检查

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# 读取数据
df = pd.read_excel('./data/aqi_20240502_1000.xlsx')

# 基础数据检查
print(f"数据形状: {df.shape}")
print(f"城市数量: {df['城市'].nunique()}")
print(f"数据时间范围: {df['时间戳'].min()} 到 {df['时间戳'].max()}")

# 检查缺失值
missing_values = df.isnull().sum()
print("缺失值统计:")
print(missing_values[missing_values > 0])

# 数据分布概览
print("\nAQI统计描述:")
print(df['AQI'].describe())

# 异常值检测（基于3σ原则）
mean_aqi = df['AQI'].mean()
std_aqi = df['AQI'].std()
outliers = df[(df['AQI'] < mean_aqi - 3*std_aqi) | (df['AQI'] > mean_aqi + 3*std_aqi)]
print(f"\n异常值数量: {len(outliers)}")

数据质量是可视化准确性的基础。我通常会进行以下几个检查：

• 完整性检查：确保所有必要字段都有值
• 范围验证：经度应在73°E-135°E，纬度应在3°N-54°N之间
• 一致性检查：同一城市不同时间的坐标是否一致

2.2 空间分布初步探索

在投入复杂的地图绘制前，先用简单的散点图了解数据分布特征：

def plot_aqi_scatter(df, save_path='./output/aqi_scatter.png'):
    """
    绘制AQI空间分布散点图
    """
    plt.figure(figsize=(12, 8))
    
    # 创建颜色映射
    norm = plt.Normalize(df['AQI'].min(), df['AQI'].max())
    cmap = plt.cm.viridis
    
    # 绘制散点
    scatter = plt.scatter(
        df['经度'], 
        df['纬度'], 
        c=df['AQI'], 
        cmap=cmap,
        norm=norm,
        s=50,           # 点大小
        alpha=0.7,      # 透明度
        edgecolors='w',  # 边缘颜色
        linewidths=0.5   # 边缘线宽
    )
    
    # 添加颜色条
    cbar = plt.colorbar(scatter, shrink=0.8)
    cbar.set_label('AQI指数', fontsize=12)
    
    # 设置图形属性
    plt.xlabel('经度 (°E)', fontsize=12)
    plt.ylabel('纬度 (°N)', fontsize=12)
    plt.title('全国338城市AQI分布散点图', fontsize=14, fontweight='bold')
    plt.grid(True, alpha=0.3, linestyle='--')
    
    # 添加中国大致范围框
    plt.axvline(x=73, color='gray', linestyle='--', alpha=0.5)
    plt.axvline(x=135, color='gray', linestyle='--', alpha=0.5)
    plt.axhline(y=3, color='gray', linestyle='--', alpha=0.5)
    plt.axhline(y=54, color='gray', linestyle='--', alpha=0.5)
    
    plt.tight_layout()
    plt.savefig(save_path, dpi=300, bbox_inches='tight')
    plt.show()
    
    return scatter

# 执行绘制
scatter_plot = plot_aqi_scatter(df)

这个散点图虽然简单，但能快速揭示几个关键信息：

1. 数据密度：东部地区监测点明显更密集
2. 异常值：是否有明显偏离主要分布区的点
3. 空间趋势：是否呈现东高西低或南高北低的趋势

3. 克里金插值：从离散点到连续场

散点图只能显示监测点的位置和数值，但我们需要的是全国范围的连续分布图。这就是空间插值的用武之地。在气象和环境领域，克里金（Kriging）插值是最常用的方法之一，因为它能考虑空间自相关性。

3.1 克里金插值原理与参数选择

克里金插值基于地统计学，基本假设是空间上接近的点比远离的点更相似。它不仅要估计未知点的值，还要给出估计的方差（不确定性）。

from pykrige.ok import OrdinaryKriging
import numpy as np

def perform_kriging(df, grid_resolution=0.5):
    """
    执行普通克里金插值
    
    参数:
    df: 包含经度、纬度、AQI的DataFrame
    grid_resolution: 网格分辨率（度）
    
    返回:
    grid_lon: 网格经度
    grid_lat: 网格纬度  
    grid_aqi: 插值后的AQI网格
    """
    
    # 提取数据
    lon = df['经度'].values
    lat = df['纬度'].values
    aqi = df['AQI'].values
    
    # 定义插值网格范围（覆盖中国全境）
    lon_min, lon_max = 73, 135
    lat_min, lat_max = 3, 54
    
    # 创建网格
    grid_lon = np.arange(lon_min, lon_max + grid_resolution, grid_resolution)
    grid_lat = np.arange(lat_min, lat_max + grid_resolution, grid_resolution)
    
    print(f"网格尺寸: {len(grid_lon)} × {len(grid_lat)} = {len(grid_lon)*len(grid_lat)}个格点")
    print(f"经度范围: {grid_lon[0]:.2f}°E 到 {grid_lon[-1]:.2f}°E")
    print(f"纬度范围: {grid_lat[0]:.2f}°N 到 {grid_lat[-1]:.2f}°N")
    
    # 选择变差函数模型
    # linear: 线性模型，适合趋势明显的数据
    # gaussian: 高斯模型，适合平滑变化的数据
    # exponential: 指数模型，通用性较好
    # spherical: 球状模型，地质学常用
    
    variogram_model = 'exponential'  # 根据数据特性选择
    
    # 执行克里金插值
    print("开始克里金插值计算...")
    OK = OrdinaryKriging(
        lon, lat, aqi,
        variogram_model=variogram_model,
        verbose=True,      # 显示计算过程
        enable_plotting=False,  # 不显示变差函数图
        coordinates_type='geographic'  # 指定为地理坐标
    )
    
    # 生成网格数据
    grid_aqi, ss = OK.execute('grid', grid_lon, grid_lat)
    
    print(f"插值完成，AQI范围: {grid_aqi.min():.1f} - {grid_aqi.max():.1f}")
    
    # 创建网格坐标矩阵
    grid_lon_mesh, grid_lat_mesh = np.meshgrid(grid_lon, grid_lat)
    
    return grid_lon_mesh, grid_lat_mesh, grid_aqi

# 执行插值
grid_lon, grid_lat, grid_aqi = perform_kriging(df, grid_resolution=0.5)

3.2 插值结果的质量评估

插值完成后，我们需要评估结果的质量。一个简单的方法是交叉验证：

def cross_validate_kriging(df, test_ratio=0.2):
    """
    交叉验证克里金插值精度
    """
    from sklearn.model_selection import train_test_split
    from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
    
    # 划分训练集和测试集
    train_df, test_df = train_test_split(
        df, test_size=test_ratio, random_state=42
    )
    
    print(f"训练集: {len(train_df)}个点，测试集: {len(test_df)}个点")
    
    # 用训练集训练克里金模型
    OK = OrdinaryKriging(
        train_df['经度'].values,
        train_df['纬度'].values,
        train_df['AQI'].values,
        variogram_model='exponential',
        verbose=False
    )
    
    # 在测试集位置进行预测
    pred_aqi, pred_var = OK.execute(
        'points',
        test_df['经度'].values,
        test_df['纬度'].values
    )
    
    # 计算评估指标
    mse = mean_squared_error(test_df['AQI'].values, pred_aqi)
    mae = mean_absolute_error(test_df['AQI'].values, pred_aqi)
    r2 = r2_score(test_df['AQI'].values, pred_aqi)
    
    print("\n交叉验证结果:")
    print(f"均方误差 (MSE): {mse:.2f}")
    print(f"平均绝对误差 (MAE): {mae:.2f}")
    print(f"决定系数 (R²): {r2:.3f}")
    
    # 可视化预测 vs 实际
    plt.figure(figsize=(8, 6))
    plt.scatter(test_df['AQI'].values, pred_aqi, alpha=0.6)
    plt.plot([0, 300], [0, 300], 'r--', label='理想线')
    plt.xlabel('实际AQI')
    plt.ylabel('预测AQI')
    plt.title('克里金插值交叉验证')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.savefig('./output/kriging_cv.png', dpi=300)
    plt.show()
    
    return mse, mae, r2

# 执行交叉验证
metrics = cross_validate_kriging(df)

4. 地图绘制核心：Cartopy与frykit的深度集成

有了插值数据，现在进入最核心的部分——地图绘制。这里我会详细介绍如何将Cartopy的基础功能与frykit的中国地图优化结合起来。

4.1 基础地图框架搭建

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import frykit.plot as fplt
import frykit.shp as fshp
import numpy as np

def create_china_map_base(figsize=(12, 10), dpi=150):
    """
    创建中国地图基础框架
    
    返回:
    fig: 图形对象
    ax_main: 主地图坐标轴
    ax_south_china_sea: 南海小图坐标轴
    """
    
    # 设置中文字体（根据系统调整）
    plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'DejaVu Sans']
    plt.rcParams['axes.unicode_minus'] = False
    
    # 创建图形和主地图
    fig = plt.figure(figsize=figsize, dpi=dpi)
    
    # 选择投影 - 兰勃特等角圆锥投影，适合中国区域
    map_crs = ccrs.LambertConformal(
        central_longitude=105,      # 中央经线
        standard_parallels=(25, 47) # 标准纬线
    )
    
    # 数据坐标系统（经纬度）
    data_crs = ccrs.PlateCarree()
    
    # 主地图范围（中国主要区域）
    main_extent = [73, 136, 3, 54]  # 东经73-136°，北纬3-54°
    
    # 创建主地图
    ax_main = fig.add_subplot(1, 1, 1, projection=map_crs)
    ax_main.set_extent(