Python pcolormesh 热图在绘图上不可见（南太平洋每月热带气旋路径气候学）

Question

我正在尝试制作南太平洋盆地每月热带气旋 (TC) 路径气候学的 Python 图。但是，我本应使用 pcolormesh 制作的每个网格 TC 轨道通过频率的网格热图并未出现在绘图上，从而导致地图具有纯白色背景。

显然，这个问题与处理跨越反子午线（180°）的数据和地图的问题有关，其中我不太熟悉所需的解决方法。我在其他地区（例如西北太平洋、印度洋、北大西洋、北太平洋东部）做了同样的绘图，没有任何问题，因为这些地区都没有穿过反子午线。我正在使用来自 IBTrACS（国际气候管理最佳跟踪档案） 数据集的数据。

这是我到目前为止所做的：

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 LONGITUDE_FORMATTER, LATITUDE_FORMATTER

# Extracting TC best track data from the IBTrACS CSV files
all_TC_tracks = pd.read_csv('C:/Users/Brian/Desktop/ibtracs.SP.list.v04r00.csv',
                      usecols=['SID', 'ISO_TIME', 'LAT', 'LON'])

# Converting the string times into datetime type
all_TC_tracks['ISO_TIME'] = pd.to_datetime(all_TC_tracks['ISO_TIME'])

# Extracting the month value from the datetimes
all_TC_tracks['MONTH'] = all_TC_tracks['ISO_TIME'].dt.month

# Converting the longitude values from [-180, 180] notation to [0, 360] notation
all_TC_tracks['LON'] = np.where(all_TC_tracks['LON'] < 0, all_TC_tracks['LON'] + 360, all_TC_tracks['LON'])
# This is necessary since the Southern Pacific basin crosses the antimeridian (180°).

# Determining the month of occurence of each TC
# by identifying the mode of all month values per TC best track
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)

def plot_monthly_tracks(month):
    # Filter the DataFrame to only include rows for a specific month
    monthly_tracks = all_TC_tracks[all_TC_tracks['MONTH'] == month]

    # Group the DataFrame by 'SID'
    grouped = monthly_tracks.groupby('SID')
    
    # Initialize an empty list to store the dictionaries
    monthly_TCs = []
    
    # Iterate over the grouped DataFrame
    for name, group in grouped:
        # Create a dictionary where the key is the SID and the value is a list of lat-lon tuples
        tc_dict = {name: list(zip(group['LAT'], group['LON']))}
        # Append the dictionary to the list
        monthly_TCs.append(tc_dict)
    
    # Create a 2D array for the histogram
    hist = np.zeros((75, 95)) #70S-5N, 155E-110W
    print(hist)
    
    # Iterate over the tropical cyclones
    for tc in monthly_TCs:
        for sid, locations in tc.items():
            # Create a set to store the grid boxes that this tropical cyclone passes through
            visited_boxes = set()
            for lat, lon in locations:
                # Calculate the indices of the grid box that this location falls into
                lat_index = int(lat - (-70))
                lon_index = int(lon - 155)
                # Check if the indices are within the expected range
                if 0 <= lat_index < 75 and 0 <= lon_index < 95:
                    # Add the grid box to the set of visited boxes
                    visited_boxes.add((lat_index, lon_index))
            # Increment the count in the histogram for each visited box
            for lat_index, lon_index in visited_boxes:
                hist[lat_index, lon_index] += 1
                # print(hist[lat_index, lon_index])
    
    # Define the edges of the bins for the histogram
    xedges = np.arange(155, 251, 1)  # 95 elements
    yedges = np.arange(-70, 6, 1)  # 75 elements
    
    print(xedges)
    print(yedges)
    
    fig = plt.figure(figsize=(25., 25.), dpi=250)
    revised_projection = ccrs.PlateCarree(central_longitude=180)
    ax = plt.axes(projection=revised_projection)
    ax.set_extent([155, 250, -70, 5], crs=ccrs.PlateCarree())
    
    gridlines = ax.gridlines(draw_labels=False, xlocs=np.arange(-180, 181, 5), ylocs=np.arange(-90, 91, 5), color='gray', linestyle='--')
    xticks = np.arange(155, 251, 5)
    yticks = np.arange(-70, 6, 5)
    ax.set_xticks(xticks, crs=ccrs.PlateCarree())
    ax.set_yticks(yticks, crs=ccrs.PlateCarree())
    
    ax.xaxis.set_major_formatter(LONGITUDE_FORMATTER)
    ax.yaxis.set_major_formatter(LATITUDE_FORMATTER)
    
    ax.coastlines('10m', edgecolor='black', linewidth=2.5)
    ax.add_feature(cfeature.BORDERS, edgecolor='black', linewidth=2.5)
    
    # Plot the histogram as a heatmap on the geographical map
    cax = ax.pcolormesh(xedges, yedges, hist, vmin=0, vmax=50, cmap='PuRd', shading='auto')
    print(cax)
    
    # Add a colorbar
    cax = fig.add_axes([0.92, 0.25, 0.02, 0.5])
    norm = Normalize(vmin=0, vmax=50)
    colorbar = ColorbarBase(cax, cmap='PuRd', norm=norm, extend='max', ticks=np.arange(0, 51, 2))
    colorbar.set_label('Frequency of tropical cyclone passage', size=25)
    cax.tick_params(labelsize=15)
    
    # Iterate over the tropical cyclones
    for tc in monthly_TCs:
        for sid, locations in tc.items():
            # Unzip the list of lat-lon tuples into two separate lists
            lats, lons = zip(*locations)
            # Plot the track of the tropical cyclone
            ax.plot(lons, lats, color='black', linewidth=0.2, transform=ccrs.PlateCarree())
    
    # Add a title and subtitle
    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

# List of months

months = ['January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September', 'October', 'November', 'December']

# Generate a plot for each month

for month in months:
    hist = plot_monthly_tracks(month)
    # break

# End of program

我尝试检查“hist”变量中二维直方图的内容，认为它是空的。但是，Spyder 中的变量资源管理器表明它不为空。

“hist”变量的内容

我还尝试打印变量“xedges”和“yedges”的内容，但控制台显示它们包含直方图中每个“框”的单元格边界的预期值。

此外，经度标签渲染不正确：反子午线以西（东半球内）的经度线标记为°W，而反子午线以东（西半球内）的经度线标记为° E.这是在我尝试将经度值从 [-180, 180] 表示法转换为 [0, 360] 表示法并将投影置于反子午线上的中心之后（因为不这样做会导致地图显示整个南半球-180 至 180 度）。

结果图

这个问题的原因是什么？我到底该如何解决这个问题？或者有没有更好的解决方法来用Python制作这些类型的图？

非常感谢您对此事的回应和建议。谢谢！

Answer 1

是对付逆经络的问题。许多图书馆假定经度在 -180 度到 180 度之间。就您而言，由于您使用的经度范围为 [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 LONGITUDE_FORMATTER, LATITUDE_FORMATTER

dtypes = {
    'SID': str,
    'ISO_TIME': str,
    'LAT': str,
    'LON': str
}

all_TC_tracks = pd.read_csv(
    r"C:\Users\serge.degossondevare\Downloads\ibtracs.SP.list.v04r00.csv",
    usecols=['SID', 'ISO_TIME', 'LAT', 'LON'],
    dtype=dtypes,
    low_memory=False
)

all_TC_tracks['ISO_TIME'] = pd.to_datetime(all_TC_tracks['ISO_TIME'], errors='coerce', format='%Y-%m-%d %H:%M:%S')

all_TC_tracks = all_TC_tracks.dropna(subset=['ISO_TIME'])

all_TC_tracks['MONTH'] = all_TC_tracks['ISO_TIME'].dt.month

all_TC_tracks['LAT'] = pd.to_numeric(all_TC_tracks['LAT'], errors='coerce')
all_TC_tracks['LON'] = pd.to_numeric(all_TC_tracks['LON'], errors='coerce')

all_TC_tracks = all_TC_tracks.dropna(subset=['LAT', 'LON'])

all_TC_tracks['LON'] = np.where(all_TC_tracks['LON'] < 0, all_TC_tracks['LON'] + 360, all_TC_tracks['LON'])

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)

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))  
    print(hist)
    
    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) 
    yedges = np.arange(-70, 6, 1)  
    
    print(xedges)
    print(yedges)
    
    fig = plt.figure(figsize=(25., 25.), dpi=250)
    revised_projection = ccrs.PlateCarree(central_longitude=180)
    ax = plt.axes(projection=revised_projection)
    ax.set_extent([155, 250, -70, 5], crs=ccrs.PlateCarree())
    
    gridlines = ax.gridlines(draw_labels=False, xlocs=np.arange(-180, 181, 5), ylocs=np.arange(-90, 91, 5), color='gray', linestyle='--')
    xticks = np.arange(155, 251, 5)
    yticks = np.arange(-70, 6, 5)
    ax.set_xticks(xticks, crs=ccrs.PlateCarree())
    ax.set_yticks(yticks, crs=ccrs.PlateCarree())
    
    ax.xaxis.set_major_formatter(LONGITUDE_FORMATTER)
    ax.yaxis.set_major_formatter(LATITUDE_FORMATTER)
    
    ax.coastlines('10m', edgecolor='black', linewidth=2.5)
    ax.add_feature(cfeature.BORDERS, edgecolor='black', linewidth=2.5)
    
    cax = ax.pcolormesh(xedges, yedges, hist, vmin=0, vmax=50, cmap='PuRd', shading='auto', transform=ccrs.PlateCarree())
    print(cax)
    
    cax = fig.add_axes([0.92, 0.25, 0.02, 0.5])
    norm = Normalize(vmin=0, vmax=50)
    colorbar = ColorbarBase(cax, cmap='PuRd', norm=norm, extend='max', ticks=np.arange(0, 51, 2))
    colorbar.set_label('Frequency of tropical cyclone passage', size=25)
    cax.tick_params(labelsize=15)
    
    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)

这给了你（我只拍摄了部分情节，因为 SO 上的图像有 2MB 的限制）。