如何在瑞士边界多边形Shapefile的ggplot中添加GEOtiff地形数据

Hey there! 我来帮你把地形 relief 数据（GeoTIFF格式）整合到你的瑞士ggplot地图里。下面是完整的实现方案，包含代码和关键注意事项：

首先要注意的核心问题：所有空间数据的坐标参考系统（CRS）必须完全一致，否则叠加后地图会错位。下面是一步步的代码和说明：

步骤1：加载所需包 #

你原来的包已经覆盖了大部分需求，raster包可以直接处理GeoTIFF文件：

# 加载依赖包
library(rgdal)
library(readxl)
library(sp)
library(ggplot2)
library(maptools)
library(plyr)
library(raster)

步骤2：导入并预处理基础数据 #

把边界和站点数据转换成ggplot能识别的格式：

# 导入瑞士边界Shapefile并转成数据框
gb <- readOGR("swissBOUNDARIES3D_1_3_TLM_KANTONSGEBIET.shp")
gb_df <- fortify(gb) # 将sp对象转换为ggplot兼容的数据框

# 导入气象站点坐标
coord <- read_excel("SMN-Stationen_20151222.xlsx")
# 可以用head(coord)确认你的经度/纬度列名，比如假设是"lon"和"lat"

步骤3：导入并处理地形Relief数据 #

这是核心新增部分，我们需要读取GeoTIFF并统一CRS：

# 读取地形GeoTIFF文件
relief_raster <- raster("你的地形文件路径.tif") # 替换成你的实际文件路径

# 检查并统一CRS（如果地形数据和边界的CRS不一致，自动转换）
if (proj4string(relief_raster) != proj4string(gb)) {
  relief_raster <- projectRaster(relief_raster, crs = proj4string(gb))
}

# 将栅格数据转成ggplot可用的数据框（保留xy坐标）
relief_df <- as.data.frame(relief_raster, xy = TRUE)
colnames(relief_df)[3] <- "elevation" # 重命名海拔列，方便后续调用

步骤4：绘制完整地图 #

注意图层顺序：地形在最底层，然后是边界，最后是站点点，这样站点会显示在最上方：

ggplot() +
  # 底层：地形图层
  geom_raster(data = relief_df, aes(x = x, y = y, fill = elevation)) +
  # 设置海拔颜色渐变（viridis调色板色盲友好，也可以自定义）
  scale_fill_viridis_c(option = "viridis", name = "海拔(m)") +
  # 中层：瑞士边界
  geom_polygon(data = gb_df, aes(x = long, y = lat, group = group), 
               fill = NA, color = "black", linewidth = 0.5) +
  # 顶层：气象站点
  geom_point(data = coord, aes(x = 你的经度列名, y = 你的纬度列名), 
             color = "red", size = 2, alpha = 0.8) +
  # 保持地图比例正确，避免拉伸变形
  coord_fixed() +
  # 美化图表
  labs(title = "瑞士地形与气象站点分布", x = "经度", y = "纬度") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

关键注意事项 #

• CRS一致性：这是地图叠加成功的核心！如果不确定CRS，可以用proj4string(gb)和proj4string(relief_raster)分别查看，确保两者一致。
• 图层顺序：ggplot会按照代码中图层的顺序绘制，所以底层的地形要放在最前面，站点放在最后。
• 颜色自定义：如果你不喜欢viridis调色板，可以换成scale_fill_gradient(low = "lightblue", high = "darkblue")这类自定义渐变。
• 站点列名：记得把代码中的你的经度列名和你的纬度列名替换成你实际的列名（比如从head(coord)中查看的列名）。

内容的提问来源于stack exchange，提问作者T. Grüter