如何用ggplot2的geom_contour复现R基础绘图的等高线图

我来帮你解决这个ggplot2复现等高线的问题！你的核心问题在于ggplot2的geom_contour需要规则网格的x/y/z数据，而你当前构建simDF的方式没有生成正确的网格对应关系，导致等高线无法正常渲染。下面是完整的修正方案：

问题根源分析 #

你用expand.grid(fSim, fSim)生成了copula的均匀边际u和v，但后续错误地把NEP直接关联到分位数，没有保留网格的一一对应关系。基础绘图的contour函数接受矩阵格式的z值，而ggplot需要长格式的、每个(x,y)坐标对应唯一z值的数据集，这是两者的核心差异。

完整修正代码 #

先保留你原来的拟合和模拟逻辑（这部分是正确的），然后重点修正网格数据的构建和ggplot绘图部分：

library(evd) 
library(copula) 
library(dplyr)
library(ggplot2)

# 生成模拟数据（设置种子保证结果可重复）
set.seed(123)
dur <- rgev(500, 2.854659, 2.170122, -0.007829) 
mag <- rgev(500, 0.02482, 0.01996, 0.04603) 

# 拟合GEV边际分布
fDurGev <- fgev(dur) 
fMagGev <- fgev(mag) 

# 将原始数据转换为copula需要的均匀边际
durVec <- pgev(dur, fDurGev[[1]][1], fDurGev[[1]][2], fDurGev[[1]][3]) 
magVec <- pgev(mag, fMagGev[[1]][1], fMagGev[[1]][2], fMagGev[[1]][3]) 
durMagMat <- as.matrix(cbind(duration = durVec, magnitude = magVec)) 

# 拟合Clayton Copula
theta <- coef(fitCopula(claytonCopula(dim = 2), durMagMat, method = "itau")) 
clayCop <- claytonCopula(theta, dim = 2) 

# 计算观测数据的copula概率指标
fCopDurMag <- pCopula(durMagMat, clayCop) 
copPts <- data.frame(duration = dur, magnitude = mag, 
                     copNEP = fCopDurMag, copEP = (1 - fCopDurMag), copRI = (1 / fCopDurMag)) 

# 生成模拟的均匀边际序列（减少长度加快计算，原1000也可使用）
fSim <- seq(0.05, 0.99998, length.out = 100) 
quaDur <- qgev(fSim, fDurGev[[1]][1], fDurGev[[1]][2], fDurGev[[1]][3]) 
quaMag <- qgev(fSim, fMagGev[[1]][1], fMagGev[[1]][2], fMagGev[[1]][3]) 

# 构建均匀边际的网格，并计算对应的copula CDF值（NEP）
expDurMagMat <- expand.grid(u = fSim, v = fSim) 
simPred <- pCopula(as.matrix(expDurMagMat), clayCop) 

# -------------------------- 重点修正：构建ggplot需要的网格数据 --------------------------
simDF <- expDurMagMat %>%
  mutate(
    quaDur = qgev(u, fDurGev[[1]][1], fDurGev[[1]][2], fDurGev[[1]][3]),
    quaMag = qgev(v, fMagGev[[1]][1], fMagGev[[1]][2], fMagGev[[1]][3]),
    NEP = simPred
  )

# 生成重现期对应的NEP阈值
RI <- c(1.25, 2 ,5, 10, 20, 50, 100, 200, 500)
NEP_levels <- 1 - (1 / RI)

# 生成copula模拟样本（用于背景点）
rndPred <- data.frame(rCopula(5000, clayCop))
rndPred$rndDur <- qgev(rndPred[,1], fDurGev[[1]][1], fDurGev[[1]][2], fDurGev[[1]][3])
rndPred$rndMag <- qgev(rndPred[,2], fMagGev[[1]][1], fMagGev[[1]][2], fMagGev[[1]][3])

ggplot绘图代码 #

ggplot() +
  # 绘制copula模拟的背景点
  geom_point(data = rndPred, aes(x = rndDur, y = rndMag), 
             color = "light grey", alpha = 0.5, size = 0.8) +
  # 绘制原始观测点
  geom_point(data = copPts, aes(x = duration, y = magnitude), 
             color = "red", size = 0.8) +
  # 绘制重现期等高线
  geom_contour(data = simDF, aes(x = quaDur, y = quaMag, z = NEP),
               breaks = NEP_levels, color = "black", linewidth = 1) +
  # 给等高线添加重现期标签（直接显示RI而非NEP）
  geom_contour_label(data = simDF, aes(x = quaDur, y = quaMag, z = NEP),
                     breaks = NEP_levels, 
                     label = sapply(NEP_levels, function(x) round(1/(1-x), 1)),
                     color = "black", size = 3) +
  labs(x = "Duration (time)", y = "Magnitude (x)") +
  theme_classic() +
  # 匹配基础绘图的坐标轴范围
  coord_cartesian(xlim = range(quaDur), ylim = range(quaMag))

关键修改说明 #

1. 保留规则网格结构：通过expand.grid(u, v)确保每个均匀边际对(u,v)对应唯一的(quaDur, quaMag)和NEP值，形成ggplot能识别的规则网格
2. 直接显示重现期：用geom_contour_label将NEP转换为直观的重现期RI标签，无需手动转换
3. 优化视觉效果：调整点的大小、透明度，让图层层次更清晰
4. 统一坐标轴范围：用coord_cartesian确保和基础绘图的坐标轴范围一致，保持图表一致性

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