R语言中如何绘制含NA值的多站点长时序时间序列图

一直在R语言环境中查找缺失值可视化的相关实现方案，尽管现有可选的优质方案数量众多，但始终未找到能完全匹配自身需求的可运行代码。

示例测试数据 #

测试用数据框df的构造代码如下：

data.frame(
  stringsAsFactors = FALSE,
              Date = c("01/01/2000","02/01/2000",
                       "03/01/2000","04/01/2000","05/01/2000","06/01/2000",
                       "07/01/2000","08/01/2000","09/01/2000"),
            Site.1 = c(NA,0.952101337,0.066766616,
                       0.77279551,0.715427011,NA,NA,NA,0.925705179),
            Site.2 = c(0.85847963,0.663818831,NA,NA,
                       0.568488712,0.002833073,0.349365844,0.652482654,
                       0.334879886),
            Site.3 = c(0.139854891,0.057024999,
                       0.297705256,0.914754178,NA,0.14108163,0.282896932,
                       0.823245136,0.153609705),
            Site.4 = c(0.758317946,0.284147119,
                       0.756356853,NA,NA,0.313465424,NA,0.013689324,0.654615632)
) -> df

预期呈现效果 #

绘制结果需与参考样图保持一致：

业务适配要求 #

实际业务数据规模较大，共包含51个监测站点、约9000个日期节点的长时序观测记录，绘图逻辑需要适配该数据规模，保证运行效率。

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