基于PyVista的地震震源3D绘图未呈现预期效果问题排查

我使用Python的pyvista包绘制地震震源的3D分布图，参考gemgis.readthedocs.io上的教程，代码可正常运行，但生成的绘图效果与参考示例不符：

我的绘图：
参考绘图：

已整理好经度、纬度、深度及震级数据，代码如下：

import geopandas as gpd
import matplotlib.pyplot as plt
import pyvista as pv
import numpy as np
import gemgis as gg
df = gpd.read_file("data gempa 2010-2021 talaud.csv")
gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.X, df.Y))
gdf = gg.vector.extract_xy(gdf=gdf)
test = pv.Sphere(radius=1000, center=gdf.loc[0][['X', "Y", "Z"]].tolist())


spheres = pv.MultiBlock([pv.Sphere(radius=float(gdf.loc[i]["Magnitude"])*200, center=gdf.loc[i][['X', 'Y', 'Z']].tolist()) for i in range(len(gdf))])
for i in range(len(spheres.keys())):
    spheres[spheres.keys()[i]]['Magnitude'] = np.zeros(len(spheres[spheres.keys()[i]].points)) + float(gdf.loc[i]['Magnitude'])


sargs = dict(fmt="%.1f", color='black')


p = pv.Plotter()
p.add_mesh(spheres,scalars='Magnitude', cmap='Reds', clim=[0,6], scalar_bar_args=sargs)
p.set_background('white')
p.show_grid(color='black')
p.enable_stereo_render()
p.show()

核心问题分析 #

• 坐标尺度不匹配：
  你的数据直接使用经纬度（X/Y为度）作为3D坐标，但深度Z的单位通常是米或千米，两者尺度差异极大——1度经度在赤道附近约等于111千米，这会导致Z方向的深度变化在画面中被完全压缩，所有点看起来都挤在一个平面上，和参考示例的立体效果差距明显。参考教程应该是先将经纬度转换为了平面直角坐标系（如UTM投影），让X/Y/Z的单位统一（比如都为米）。
• 球体半径尺度错误：
  当前设置的Magnitude*200半径，相对于经纬度的尺度来说过小，所以球体呈现为小点。而参考示例中坐标是投影后的米级单位，半径设置才能对应合适的视觉大小。
• 未做坐标投影转换：
  gg.vector.extract_xy仅提取几何点的X/Y到属性列，没有进行坐标投影转换，这是关键遗漏步骤。

修正方案 #

1. 转换到平面投影坐标系：
   先查看数据原始坐标系（执行print(gdf.crs)），通常是WGS84（EPSG:4326），然后转换为区域对应的UTM投影（比如印尼Talaud群岛附近可用EPSG:32751）：

   # 转换坐标系统
   gdf = gdf.to_crs(epsg=32751)
   # 更新X/Y为投影后的坐标值
   gdf['X'] = gdf.geometry.x
   gdf['Y'] = gdf.geometry.y
   

2. 统一深度单位：
   如果深度Z的单位是千米，转换为米以匹配投影坐标单位：

   gdf['Z'] = gdf['Z'] * 1000  # 若深度单位为千米则执行此步
   

3. 调整球体半径：
   根据投影后的米级单位调整半径，比如设置为Magnitude*1000，可根据实际视觉效果微调：

   spheres = pv.MultiBlock([pv.Sphere(radius=float(gdf.loc[i]["Magnitude"])*1000, center=gdf.loc[i][['X', 'Y', 'Z']].tolist()) for i in range(len(gdf))])
   

4. 可选：优化MultiBlock创建逻辑：
   用循环遍历的方式更直观，避免loc索引的频繁调用：

   spheres = pv.MultiBlock()
   for idx, row in gdf.iterrows():
       sphere = pv.Sphere(radius=row['Magnitude']*1000, center=[row['X'], row['Y'], row['Z']])
       sphere['Magnitude'] = np.zeros(len(sphere.points)) + row['Magnitude']
       spheres.append(sphere)
   

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