把地震目录按 0.1 度格网做成热力图，卡在颜色分级上

手上有一份 1970 年至今、约 12 万条事件的地震目录，字段只有时间、经纬度、震级、深度，CSV 大约 8MB，在笔记本上跑。想做的其实很简单：把震中按 0.1°×0.1° 的格子计数，出一张能一眼看出哪些区域事件密集的热力图。第一版直接上 matplotlib 的 hist2d，bins=360×180，出来的图中间一团黑，边缘全白，完全没法看。

问题出在计数分布：大部分格子是 0 或 1，少数俯冲带上的格子能到四位数。线性色标把所有低值都压没了。我先试了对数色标，用 LogNorm：

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm

lon = df["longitude"].to_numpy()
lat = df["latitude"].to_numpy()

fig, ax = plt.subplots(figsize=(12, 6))
h = ax.hist2d(lon, lat, bins=[np.arange(-180, 180.1, 0.1),
                            np.arange(-90, 90.1, 0.1)],
              norm=LogNorm(vmin=1, vmax=None), cmap="magma")
fig.colorbar(h[3], ax=ax, label="事件数（对数）")

对数色标出来后，密度结构清楚了，但出现两个新问题。一是 0 值格子会被 LogNorm 当成无效值处理，显示成白色，跟「没有数据」混在一起；二是 1 和 2 的格子被拉得太开，看起来像热点，其实只是偶然落进去的一两条。

我最后的处理是先把计数矩阵算出来，再手动把 0 改成 np.nan，色标从 1 开始，并且把 vmax 卡在 99 分位而不是最大值：

counts, xedges, yedges = np.histogram2d(lon, lat, bins=[360, 180],
                                        range=[[-180, 180], [-90, 90]])
counts[counts == 0] = np.nan
vmax = np.nanpercentile(counts, 99)
ax.pcolormesh(xedges, yedges, counts.T, norm=LogNorm(vmin=1, vmax=vmax))

为什么这么做：地震目录本身是被检测能力塑造过的，越早的年代、越偏远的海域，漏记越多。热力图上的「空白」不一定真的没有地震，可能只是没有台站。所以色标卡在 99 分位，是为了不让少数几个俯冲带格子吃掉整个动态范围，代价是超过 vmax 的格子会被饱和成同一个颜色，看上去会低估最密集区域的差异。

失效的地方也在这里。0.1° 在赤道附近约 11km，到高纬度经度方向会明显缩水，全球等经纬格网在高纬会变形，如果只关心某个区域，最好换成本地投影再分格。另外，同一地点不同年代的事件数放在一起时，早期目录稀疏会让「热点」偏向近几十年，这不是地质信号，是观测历史。我踩过一次坑：直接用 hist2d 的自动 bins，x 和 y 的格子宽度不一致，图上的密度对比完全失真，后来固定成 360×180 才稳定。

现在这张图还只是「事件计数」，没有按震级加权。如果改成按能量累加，色标的分布会完全不同，1 个 M6 可能顶掉上千个 M2，你们在做类似目录可视化时，是直接计数还是按能量加权？