南太平洋热带气旋月气候图pcolormesh热力图不显示及经度标签错误
南太平洋热带气旋路径热力图问题排查与修复方案
问题原因分析
1. pcolormesh热力图空白
- 坐标系映射缺失:轴使用了中心180°的
PlateCarree投影,但pcolormesh未指定transform参数,Cartopy无法将[0,360]范围的经度数据正确映射到投影坐标系,导致热力图无法渲染。 - (注:你的
hist数组维度与网格边缘匹配正确,排除维度不匹配问题)
2. 经度标签显示错误
LONGITUDE_FORMATTER默认适配[-180,180]经度范围,当传入[0,360]的刻度值时,会自动将>180的数值转换为西经(如200°→160°W),<180的转换为东经,完全颠倒了南太平洋的经度标注逻辑。
修复方案
核心修改点
- 给
pcolormesh添加transform=ccrs.PlateCarree()参数,明确数据坐标系。 - 自定义经度格式化函数,将[0,360]经度转换为标准东/西经表示。
- 优化网格线与刻度的设置逻辑,避免重复定义。
修复后完整代码
import numpy as np import pandas as pd from scipy.stats import mode import matplotlib.pyplot as plt from matplotlib.colorbar import ColorbarBase from matplotlib.colors import Normalize import cartopy.crs as ccrs import cartopy.feature as cfeature from cartopy.mpl.gridliner import LATITUDE_FORMATTER # 读取IBTrACS数据 all_TC_tracks = pd.read_csv('C:/Users/Brian/Desktop/ibtracs.SP.list.v04r00.csv', usecols=['SID', 'ISO_TIME', 'LAT', 'LON']) all_TC_tracks['ISO_TIME'] = pd.to_datetime(all_TC_tracks['ISO_TIME']) all_TC_tracks['MONTH'] = all_TC_tracks['ISO_TIME'].dt.month # 经度转换为[0,360] all_TC_tracks['LON'] = np.where(all_TC_tracks['LON'] < 0, all_TC_tracks['LON'] + 360, all_TC_tracks['LON']) # 确定每个TC的主要活跃月份 most_frequent_month = all_TC_tracks.groupby('SID')['MONTH'].agg(lambda x: mode(x)[0]) most_frequent_month = most_frequent_month.map({ 1: 'January', 2: 'February', 3: 'March', 4: 'April', 5: 'May', 6: 'June', 7: 'July', 8: 'August', 9: 'September', 10: 'October', 11: 'November', 12: 'December'}) all_TC_tracks = all_TC_tracks.merge( most_frequent_month, how='left', on='SID').rename( columns={'MONTH_y': 'MOST_FREQUENT_MONTH'}) all_TC_tracks.drop(columns=['ISO_TIME','MONTH_x'], inplace=True) all_TC_tracks.rename(columns={'MOST_FREQUENT_MONTH': 'MONTH'}, inplace=True) # 自定义经度格式化函数:[0,360] → 东/西经 def custom_lon_formatter(x, pos): if x > 180: return f"{360 - x:.0f}°W" else: return f"{x:.0f}°E" def plot_monthly_tracks(month): monthly_tracks = all_TC_tracks[all_TC_tracks['MONTH'] == month] grouped = monthly_tracks.groupby('SID') monthly_TCs = [] for name, group in grouped: tc_dict = {name: list(zip(group['LAT'], group['LON']))} monthly_TCs.append(tc_dict) # 初始化频次直方图 hist = np.zeros((75, 95)) # 70S-5N (75格), 155E-110W (95格) # 统计每个TC经过的网格 for tc in monthly_TCs: for sid, locations in tc.items(): visited_boxes = set() for lat, lon in locations: lat_index = int(lat - (-70)) lon_index = int(lon - 155) if 0 <= lat_index < 75 and 0 <= lon_index < 95: visited_boxes.add((lat_index, lon_index)) for lat_index, lon_index in visited_boxes: hist[lat_index, lon_index] += 1 # 网格边缘 xedges = np.arange(155, 251, 1) # 155到250,共96个边缘点,对应95格 yedges = np.arange(-70, 6, 1) # -70到5,共76个边缘点,对应75格 # 创建绘图对象 fig = plt.figure(figsize=(25., 25.), dpi=250) proj = ccrs.PlateCarree(central_longitude=180) ax = plt.axes(projection=proj) ax.set_extent([155, 250, -70, 5], crs=ccrs.PlateCarree()) # 设置网格线与标签 gridlines = ax.gridlines(draw_labels=True, xlocs=np.arange(155, 251, 5), ylocs=np.arange(-70, 6, 5), color='gray', linestyle='--') gridlines.xformatter = custom_lon_formatter gridlines.yformatter = LATITUDE_FORMATTER gridlines.top_labels = False gridlines.right_labels = False # 添加地理特征 ax.coastlines('10m', edgecolor='black', linewidth=2.5) ax.add_feature(cfeature.BORDERS, edgecolor='black', linewidth=2.5) # 绘制热力图(核心修复:添加transform参数) cax = ax.pcolormesh(xedges, yedges, hist, vmin=0, vmax=50, cmap='PuRd', shading='auto', transform=ccrs.PlateCarree()) # 添加色标 cax_colorbar = fig.add_axes([0.92, 0.25, 0.02, 0.5]) norm = Normalize(vmin=0, vmax=50) colorbar = ColorbarBase(cax_colorbar, cmap='PuRd', norm=norm, extend='max', ticks=np.arange(0, 51, 2)) colorbar.set_label('Frequency of tropical cyclone passage', size=25) cax_colorbar.tick_params(labelsize=15) # 绘制TC轨迹 for tc in monthly_TCs: for sid, locations in tc.items(): lats, lons = zip(*locations) ax.plot(lons, lats, color='black', linewidth=0.2, transform=ccrs.PlateCarree()) # 设置标题 fig.suptitle('IBTrACS | Monthly climatology of Southern Pacific tropical cyclone tracks (1848–2023)', fontsize=29, y=0.851) ax.set_title(month.upper(), fontsize=40, pad=20) plt.show() return hist # 批量绘制月份图 months = ['January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September', 'October', 'November', 'December'] for month in months: hist = plot_monthly_tracks(month)
更优绘图方案建议
1. 用numpy.histogram2d替代手动统计
手动循环统计网格频次效率低且易出错,可直接用numpy.histogram2d一次性完成:
# 替换原手动统计hist的代码 lats_all = [] lons_all = [] for tc in monthly_TCs: for sid, locations in tc.items(): ls, ls_lon = zip(*locations) lats_all.extend(ls) lons_all.extend(ls_lon) hist, yedges, xedges = np.histogram2d(lats_all, lons_all, bins=[np.arange(-70,6,1), np.arange(155,251,1)]) # 注意:histogram2d返回的hist形状是(len(yedges)-1, len(xedges)-1),与原代码一致
2. 用Cartopy原生工具处理跨日界线标签
Cartopy的LongitudeFormatter支持dateline_direction参数,可自动处理跨日界线的经度标签,无需自定义函数:
from cartopy.mpl.ticker import LongitudeFormatter # 替换自定义格式化函数 gridlines.xformatter = LongitudeFormatter(dateline_direction='east', degree_symbol='°')
内容的提问来源于stack exchange,提问作者Brian Añano
相关产品推荐
相关产品推荐

