如何基于Matplotlib Basemap绘制震源深度剖面图？

嘿，我明白你想要实现类似GMT的震源深度剖面图的需求——把特定方位角和宽度范围内的地震震源投影到X-Z平面上对吧？其实不用切换到GMT，咱们可以在你现有的Matplotlib + Basemap代码基础上，通过坐标转换和数据筛选来实现这个功能。下面是具体的实现步骤和代码示例：

核心思路 #

1. 定义你想要的剖面线：比如指定剖面的起点和终点经纬度（对应你说的方位角），同时设置剖面的宽度阈值。
2. 把所有地震的经纬度转换为平面投影坐标（比如用Basemap的transform方法转成地图投影坐标）。
3. 计算每个地震到剖面线的垂直距离，筛选出距离在宽度范围内的地震。
4. 计算这些筛选后的地震沿剖面线的投影距离（作为X轴），深度作为Z轴（注意深度是向下的，所以要反转Y轴）。

完整代码示例 #

import requests
from csv import DictReader
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

# 下载地震数据（和你原代码一致）
DATA_URL='http://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/4.5_month.csv'
print("Downloading", DATA_URL)
resp = requests.get(DATA_URL)
quakes = list(DictReader(resp.text.splitlines()))

# 提取并转换数据
lngs = np.array([float(q['longitude']) for q in quakes])
lats = np.array([float(q['latitude']) for q in quakes])
depths = np.array([float(q['depth']) for q in quakes])
mags = np.array([2 ** float(q['mag']) for q in quakes])

# --------------------------
# 1. 定义剖面参数
# --------------------------
# 示例：设置剖面的起点和终点经纬度（对应你想要的方位角）
# 这里随便选了一段跨太平洋的剖面，你可以根据需求修改
profile_start = (140, 35)  # (经度, 纬度)
profile_end = (-120, 35)
profile_width = 200  # 剖面宽度，单位：公里（可以根据需求调整）

# --------------------------
# 2. 初始化Basemap并转换坐标
# --------------------------
# 用Mercator投影来做平面转换，你也可以换成UTM投影更精准
m = Basemap(projection='merc', llcrnrlat=20, urcrnrlat=50,
            llcrnrlon=130, urcrnrlon=-110, resolution='l')

# 把经纬度转成地图投影的x,y坐标（单位：米）
x, y = m(lngs, lats)
# 剖面起点和终点的投影坐标
x0, y0 = m(profile_start[0], profile_start[1])
x1, y1 = m(profile_end[0], profile_end[1])

# --------------------------
# 3. 筛选剖面范围内的地震
# --------------------------
# 计算剖面线的向量
line_vec = np.array([x1 - x0, y1 - y0])
# 计算每个地震点到起点的向量
point_vec = np.array([x - x0, y - y0])

# 计算点到直线的垂直距离（单位：米）
cross_product = np.cross(line_vec, point_vec)
line_length = np.linalg.norm(line_vec)
distance_to_profile = np.abs(cross_product) / line_length

# 转换为公里，筛选出在剖面宽度内的地震
mask = distance_to_profile <= profile_width * 1000
x_selected = x[mask]
y_selected = y[mask]
depths_selected = depths[mask]
mags_selected = mags[mask]

# --------------------------
# 4. 计算沿剖面的距离（作为X轴）
# --------------------------
# 计算每个点在剖面线上的投影长度（单位：公里）
proj_length = np.dot(point_vec[mask], line_vec) / line_length
proj_length_km = proj_length / 1000

# --------------------------
# 5. 绘制震源深度剖面图
# --------------------------
plt.figure(figsize=(12, 6))
# 深度向下为正，所以反转Y轴
plt.scatter(proj_length_km, depths_selected, s=mags_selected/10, c='darkred', alpha=0.6)
plt.gca().invert_yaxis()
plt.xlabel('Distance along profile (km)')
plt.ylabel('Depth (km)')
plt.title('Earthquake Hypocenter Profile')
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

# 你也可以同时画出2D地图和剖面的位置，方便对照
plt.figure(figsize=(14, 8))
m.bluemarble(alpha=0.42)
m.drawcoastlines(color='#555566', linewidth=1)
# 画出剖面线
m.plot([x0, x1], [y0, y1], color='yellow', linewidth=3, zorder=11)
# 画出筛选出的地震点
m.scatter(x_selected, y_selected, mags_selected/10, c='red', alpha=0.5, zorder=10)
plt.show()

关键细节说明 #

• 剖面参数调整：你可以修改profile_start、profile_end来设置不同的方位角，profile_width控制剖面的横向范围。
• 投影选择：如果需要更精准的距离计算，建议换成UTM投影（Basemap支持utm投影类型，需要指定带号），Mercator投影在高纬度会有距离变形。
• 深度处理：因为地震深度是向下的，所以用invert_yaxis()让Y轴从上到下数值增大，符合常规的剖面展示习惯。
• 震级缩放：原代码里的2**mag会让震级大的点尺寸过大，我在剖面里除以了10，你可以根据视觉效果调整缩放比例。

这样就能实现类似GMT的震源深度剖面图啦，完全基于你现有的Python工具链，不用切换到其他软件~

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