如何使用Python绘制带坐标的2D围栏图 地球物理反演数据可视化咨询

前置依赖准备 #

你可以优先用Python生态下成熟的地学绘图工具栈完成需求，需要提前安装的库如下：

• 地理投影与地图底图绘制用cartopy，替代已停止维护的basemap库
• 数值计算与文件读取用numpy，文本格式的dat文件直接用numpy.loadtxt读取，二进制dat文件用numpy.fromfile或者struct库解析
• 图像文件读取用PIL.Image或者matplotlib.pyplot.imread
• 绘图渲染用matplotlib

二维平面地图叠加反演结果操作流程 #

针对dat格式数值反演结果 #

• 先确认你手头的两个端点坐标（通常为研究区的左上角和右下角经纬度/平面投影坐标），以及反演结果矩阵的行列数，矩阵行数对应y轴（纬度/北向）采样数，列数对应x轴（经度/东向）采样数
• 用numpy.linspace生成对应范围的坐标序列，再通过numpy.meshgrid生成和反演矩阵维度匹配的坐标网格
• 选择匹配你数据的地图投影：小范围研究区直接用等经纬度投影cartopy.crs.PlateCarree()即可，UTM投影数据对应选择cartopy.crs.UTM()并传入对应的带号参数
• 用pcolormesh或者contourf函数将反演矩阵和坐标网格绑定，指定transform参数和你选择的投影一致，即可准确叠加到地图上

针对图像格式反演结果 #

• 读取图像后先确认坐标方向：图像默认原点为左上角，和地理数据y轴北向增大的方向相反，可用numpy.flipud()翻转图像矩阵匹配坐标方向
• 用matplotlib.offsetbox.OffsetImage包裹图像对象，再通过AnnotationBbox指定图像两个对角的坐标，直接锚定到地图的对应位置即可
• 可以调整alpha参数设置图像透明度，方便和底图边界对照校准偏移

沿z轴分层展示操作说明 #

如果需要按深度沿z轴堆叠多层反演结果，可选择matplotlib的3D轴模式：

• 每层的x/y轴对应地理坐标，z轴设置为对应层的深度值
• 逐层调用pcolormesh并指定z轴偏移量即可实现多层堆叠效果，和你预期的实现逻辑完全匹配

最简可运行示例代码 #

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np

# 替换为你自己的研究区端点坐标，示例为左上角(100°E, 30°N)、右下角(105°E, 25°N)
lon_min, lat_max = 100, 30
lon_max, lat_min = 105, 25

# 读取你的反演dat数据，此处用随机数据模拟
# inv_data = np.loadtxt("你的反演结果文件.dat", skiprows=表头行数)
inv_data = np.random.rand(50, 50) 
ny, nx = inv_data.shape

# 生成匹配维度的坐标网格
lon_arr = np.linspace(lon_min, lon_max, nx)
lat_arr = np.linspace(lat_min, lat_max, ny)
lon_grid, lat_grid = np.meshgrid(lon_arr, lat_arr)

# 绘图
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
# 添加底图特征
ax.coastlines(resolution='10m')
ax.gridlines(draw_labels=True, linestyle='--')

# 叠加反演数据
pcm = ax.pcolormesh(lon_grid, lat_grid, inv_data, cmap='jet', alpha=0.6, transform=ccrs.PlateCarree())
plt.colorbar(pcm, label='反演参数值', shrink=0.8)
plt.show()

调试校准技巧 #

• 叠加数据前先在地图上绘制你已知的两个端点的标记点，确认坐标位置和实际研究区匹配后再叠加反演结果
• 初始测试时将反演数据透明度设置为0.5左右，和底图的河流、行政边界对照，确认无坐标偏移
• 无效值或者背景值可以提前设置为np.nan，叠加时会自动透明不遮挡底图

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