如何在PyGMT中将全球XYZ数据集转为网格并绘制指定z值？

一、PyGMT核心问题解决：将一维XYZ数据转为二维网格 #

你的xarray结构是一维（index维度），而PyGMT的grdimage需要二维网格数据（以经度、纬度为维度），因此必须先完成格式转换，以下两种方法按需选择：

方法1：针对规则网格数据直接重塑（高效无插值） #

从数据维度（16200个点）和经纬度范围判断，数据是规则间隔网格（经度间隔2°，共180个点；纬度间隔2°，共90个点），可直接重塑为二维数组：

import pandas as pd
import xarray as xr
import pygmt

# 读取XYZ数据（根据实际文件格式调整分隔符）
data = pd.read_csv("your_data.xyz", sep="\s+", header=0)

# 指定要展示的z字段（比如LAB）
target_z = "LAB"

# 提取唯一经纬度值，确保顺序正确
lons = data["LONG"].unique()
lats = data["LAT"].unique()

# 将z值重塑为二维数组（维度：lat × lon）
z_array = data[target_z].values.reshape(len(lats), len(lons))

# 创建PyGMT可识别的xarray网格
grid = xr.DataArray(
    data=z_array,
    dims=["lat", "lon"],
    coords={"lat": lats, "lon": lons}
)

# 绘制图像
fig = pygmt.Figure()
fig.basemap(region="d", projection="N12c", frame=True)
fig.grdimage(grid=grid)
fig.coast(shorelines="0.5p,black")
fig.show()

方法2：用PyGMT工具转换（支持规则/不规则网格） #

如果数据是不规则网格，或不确定规则性，使用xyz2grd（规则网格）或sphinterpolate（球面插值，适合全球数据）：

import pandas as pd
import pygmt

data = pd.read_csv("your_data.xyz", sep="\s+", header=0)
target_z = "LAB"

# 提取LONG, LAT, 目标z字段，按顺序排列
xyz_data = data[["LONG", "LAT", target_z]]

# 方法A：规则网格用xyz2grd
grid = pygmt.xyz2grd(
    data=xyz_data,
    region="d",  # 全球范围
    spacing=2,   # 数据间隔，根据实际调整
    output="temp_grid.nc"  # 可选：保存网格到文件
)

# 方法B：不规则网格用球面插值sphinterpolate
# grid = pygmt.sphinterpolate(
#     data=xyz_data,
#     region="d",
#     spacing=2,
#     method="cubic",  # 插值方法：线性/立方等
#     radius=1000      # 搜索半径（单位：km）
# )

# 绘制图像
fig = pygmt.Figure()
fig.basemap(region="d", projection="N12c", frame=True)
fig.grdimage(grid=grid)
fig.coast(shorelines="0.5p,black")
fig.show()

二、Matplotlib Basemap警告解决 #

Basemap的contourf不支持tri参数，若要处理散点/非网格数据，改用plt.tricontourf：

import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import numpy as np
import pandas as pd

data = pd.read_csv("your_data.xyz", sep="\s+", header=0)
target_z = "LAB"

map = Basemap(projection='ortho', lat_0=15, lon_0=160, resolution='l')
map.drawcoastlines(linewidth=0.25)

# 将经纬度转换为地图投影坐标
x, y = map(data["LONG"].values, data["LAT"].values)

# 生成合理的等值线区间
clevs = np.linspace(data[target_z].min(), data[target_z].max(), 70)

# 使用tricontourf绘制散点数据的填色图
plt.tricontourf(x, y, data[target_z].values, clevs)
plt.colorbar(label=target_z)
plt.title(f"{target_z} 全球分布（西太平洋中心）")
plt.show()

内容的提问来源于stack exchange，提问作者Geo_toto