使用Matplotlib和Cartopy绘制台风风场圆/椭圆不显示问题求助

问题 #

我用Matplotlib和Cartopy绘制2023年台风苏拉的路径，数据来自日本气象厅（JMA）的CSV最佳路径文件。路径线条和中心点位能正常显示，但代表大风、暴风风场的圆/椭圆完全看不到。我尝试用matplotlib.patches.Circle和Ellipse类，以台风中心经纬度为圆心，把海里转成公里（乘1.852）作为半径/半轴参数，但地图上只有路径、网格和海岸线，风场图形不显示。

问题根源 #

核心错误是在PlateCarree投影（经纬度坐标）下，直接用公里作为Circle/Ellipse的半径参数无效。Matplotlib的Patch类默认使用轴坐标（此处为经纬度，单位是度），而非真实距离单位。你传入的公里数（比如150海里=277.8公里）对应到经纬度上是277.8度，远远超出地图显示范围（115-130经度，10-25纬度），图形被绘制到可视区域外，自然无法看到。

解决方案 #

要在Cartopy中绘制真实距离的圆/椭圆，需基于大地测量计算，将距离转换为符合投影的坐标。推荐以下方法：

修正后的代码 #

# Program for extracting TC track data from a comma-separated values (CSV) file
# esp. JMA best track data

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from shapely.geometry import Point
from shapely.affinity import rotate

def angle(direction):
    if direction == 'N':
        return 0
    elif direction == 'NE':
        return 45
    elif direction == 'E':
        return 90
    elif direction == 'SE':
        return 135
    elif direction == 'S':
        return 180
    elif direction == 'SW':
        return 225
    elif direction == 'W':
        return 270
    elif direction == 'NW':
        return 315

# Importing TC best track data from a CSV file
saola_besttrack = pd.read_csv('C:/Users/Brian/Desktop/saola_track.csv')

# Setting the latitude and longitude values of the PAR region
PAR_longitudes = [120, 135, 135, 115, 115, 120, 120]
PAR_latitudes = [25, 25, 5, 5, 15, 21, 25]

# 初始化大地测量对象，用于计算真实球面距离
geod = ccrs.Geodetic().globe.get_geod()

for i in range(len(saola_besttrack)):
    subset = saola_besttrack.iloc[:i+1]
    
    # 按风力等级拆分路径数据
    TD_track = subset[subset['Wind (kt)'] <= 33]
    TS_track = subset[(subset['Wind (kt)'] >= 34) & (subset['Wind (kt)'] <= 47)]
    STS_track = subset[(subset['Wind (kt)'] >= 48) & (subset['Wind (kt)'] <= 63)]
    TY_track = subset[(subset['Wind (kt)'] >= 64) & (subset['Wind (kt)'] <= 99)]
    STY_track = subset[subset['Wind (kt)'] >= 100]
    
    # 初始化地图
    fig = plt.figure(figsize=(25., 25.), dpi=250)
    ax = plt.axes(projection=ccrs.PlateCarree())
    ax.set_extent([115, 130, 10, 25], ccrs.PlateCarree())
    
    # 添加网格线
    gridlines = ax.gridlines(draw_labels=True, xlocs=np.arange(-180, 181, 2),
                             ylocs=np.arange(-90, 91, 2), color='gray',
                             linestyle='--')
    gridlines.top_labels = False
    gridlines.right_labels = False
    gridlines.xlabel_style = {'size': 20}
    gridlines.ylabel_style = {'size': 20}
    
    # 添加海岸线
    ax.coastlines('10m', edgecolor='black', linewidth=2.5)
    
    # 绘制PAR边界
    ax.plot(PAR_longitudes, PAR_latitudes, color='black', linewidth=6,
            transform=ccrs.PlateCarree(), linestyle='-.')
    
    # 获取当前时刻的台风数据
    current_data = subset.iloc[i]
    lon, lat = current_data['Long.'], current_data['Lat.']
    
    # 绘制大风风场
    if not pd.isna(current_data['Direc. of Major Gale Axis']):
        # 转换海里为米（1海里=1852米）
        major_radius_m = current_data['Radius of Major Gale Axis (nm)'] * 1852
        minor_radius_m = current_data['Radius of Minor Gale Axis (nm)'] * 1852
        
        if current_data['Direc. of Major Gale Axis'] == 'symmetric':
            # 创建真实球面距离的圆
            circle = Point(lon, lat).buffer(major_radius_m, resolution=32, metric=True, geod=geod)
            ax.add_geometries([circle], crs=ccrs.PlateCarree(), edgecolor='yellow', facecolor='none', linewidth=2)
        else:
            # 创建椭圆：先按长半轴建圆，缩放短半轴后旋转
            ellipse = Point(lon, lat).buffer(major_radius_m, resolution=32, metric=True, geod=geod)
            ellipse = ellipse.scale(xfact=minor_radius_m/major_radius_m, yfact=1.0, origin=(lon, lat))
            rot_angle = angle(current_data['Direc. of Major Gale Axis'])
            ellipse = rotate(ellipse, rot_angle, origin=(lon, lat), use_radians=False)
            ax.add_geometries([ellipse], crs=ccrs.PlateCarree(), edgecolor='yellow', facecolor='none', linewidth=2)
    
    # 绘制暴风风场
    if not pd.isna(current_data['Direc. of Major Storm Axis']):
        major_radius_m = current_data['Radius of Major Storm Axis (nm)'] * 1852
        minor_radius_m = current_data['Radius of Minor Storm Axis (nm)'] * 1852
        
        if current_data['Direc. of Major Storm Axis'] == 'symmetric':
            circle = Point(lon, lat).buffer(major_radius_m, resolution=32, metric=True, geod=geod)
            ax.add_geometries([circle], crs=ccrs.PlateCarree(), edgecolor='red', facecolor='none', linewidth=2)
        else:
            ellipse = Point(lon, lat).buffer(major_radius_m, resolution=32, metric=True, geod=geod)
            ellipse = ellipse.scale(xfact=minor_radius_m/major_radius_m, yfact=1.0, origin=(lon, lat))
            rot_angle = angle(current_data['Direc. of Major Storm Axis'])
            ellipse = rotate(ellipse, rot_angle, origin=(lon, lat), use_radians=False)
            ax.add_geometries([ellipse], crs=ccrs.PlateCarree(), edgecolor='red', facecolor='none', linewidth=2)
    
    # 绘制台风路径
    ax.plot(subset['Long.'], subset['Lat.'], '-k', transform=ccrs.PlateCarree(),
            linewidth=5)
    ax.plot(TD_track['Long.'], TD_track['Lat.'], 'o', color='blue',
            transform=ccrs.PlateCarree(), markersize=25)
    ax.plot(TS_track['Long.'], TS_track['Lat.'], 'o', color='green',
            transform=ccrs.PlateCarree(), markersize=25)
    ax.plot(STS_track['Long.'], STS_track['Lat.'], 'o', color='orange',
            transform=ccrs.PlateCarree(), markersize=25)
    ax.plot(TY_track['Long.'], TY_track['Lat.'], 'o', color='red',
            transform=ccrs.PlateCarree(), markersize=25)
    ax.plot(STY_track['Long.'], STY_track['Lat.'], 'o', color='purple',
            transform=ccrs.PlateCarree(), markersize=25)
    
    plt.show()

# End of program

关键修改说明 #

1. 引入Shapely库：用Point.buffer()创建基于真实球面距离的圆，通过geod参数确保距离计算符合地球曲面特性。
2. 单位转换：直接将海里转为米（1海里=1852米），匹配Shapely的距离单位要求。
3. 椭圆绘制逻辑：先创建长半轴对应的圆，缩放得到短半轴比例，再按指定方向旋转，保证椭圆的方向和大小与数据一致。
4. 代码简化：用subset.iloc[i]直接获取当前时刻数据，减少重复调用，提升可读性。

内容的提问来源于stack exchange，提问作者Brian Año