基于.shp与.tif文件绘制地图视图及高程剖面的实现方案

已备好坐标系匹配的route.shp路径文件和dem.tif地形栅格文件，以下是简洁实现代码：

一、依赖包安装与加载 #

替换老旧的rgdal为更现代的sf（空间数据处理），配合ggplot2（静态可视化）、plotly（交互式可视化）完成需求：

# 首次运行需安装依赖包
install.packages(c("sf", "ggplot2", "plotly", "raster", "patchwork", "viridis", "htmlwidgets"))

# 加载包
library(sf)
library(ggplot2)
library(plotly)
library(raster)
library(patchwork)
library(viridis)
library(htmlwidgets)

二、数据读取 #

用sf的st_read替代readOGR，处理空间数据更高效稳定：

# 读取路径矢量文件
route <- st_read(".", "route")

# 读取DEM栅格文件
dem <- raster("dem.tif")

三、静态可视化（ggplot2风格，1列2行） #

生成地形路径叠加图+路径高程剖面图的组合图，保存为plot.png：

# 1. 地形+路径叠加图
p1 <- ggplot() +
  geom_raster(data = as.data.frame(dem, xy = TRUE), aes(x = x, y = y, fill = layer)) +
  scale_fill_viridis_c(option = "magma", name = "高程(m)") +
  geom_sf(data = route, color = "#ff4d4d", linewidth = 1) +
  theme_minimal() +
  labs(title = "地形与路径分布", x = "经度", y = "纬度")

# 2. 路径高程剖面图
# 提取路径节点的高程值
route_elev <- extract(dem, route, along = TRUE)[[1]]
# 计算路径节点的累计距离
route_coords <- st_coordinates(route)[,1:2]
route_dist <- c(0, cumsum(sqrt(diff(route_coords[,1])^2 + diff(route_coords[,2])^2)))
# 绘制剖面图
p2 <- ggplot(data.frame(dist = route_dist, elev = route_elev), aes(x = dist, y = elev)) +
  geom_line(color = "#2660a4", linewidth = 1) +
  theme_minimal() +
  labs(title = "路径高程剖面", x = "路径距离(m)", y = "高程(m)")

# 组合图并保存
(p1 / p2) + plot_layout(heights = c(2, 1))
ggsave("plot.png", width = 8, height = 10, dpi = 300)

四、交互式可视化（plotly同步滑块） #

生成带同步滑块的双图交互页面，点击地形图上的路径节点，剖面图会自动定位到对应位置，保存为plot.html：

# 转换数据格式适配plotly
dem_df <- as.data.frame(dem, xy = TRUE)
route_df <- as.data.frame(st_coordinates(route))
names(route_df) <- c("x", "y")
route_df$elev <- route_elev
route_df$dist <- route_dist

# 地形路径交互图
plot_map <- plot_ly() %>%
  add_raster(data = dem_df, x = ~x, y = ~y, z = ~layer, colors = viridis(20), name = "高程") %>%
  add_lines(data = route_df, x = ~x, y = ~y, color = I("#ff4d4d"), line = list(width = 3), name = "路径") %>%
  layout(title = "地形与路径分布", xaxis = list(title = "经度"), yaxis = list(title = "纬度"))

# 高程剖面交互图
plot_profile <- plot_ly(data = route_df, x = ~dist, y = ~elev) %>%
  add_lines(color = I("#2660a4"), line = list(width = 3), name = "高程剖面") %>%
  layout(title = "路径高程剖面", xaxis = list(title = "路径距离(m)"), yaxis = list(title = "高程(m)"))

# 绑定双图交互事件，实现滑块同步
subplot(plot_map, plot_profile, nrows = 2, heights = c(0.7, 0.3)) %>%
  layout(showlegend = FALSE) %>%
  onRender("function(el, x) {
    var mapPlot = document.querySelectorAll('.plotly')[0];
    var profilePlot = document.querySelectorAll('.plotly')[1];
    mapPlot.on('plotly_click', function(event) {
      if(event.points.length > 0) {
        var idx = event.points[0].pointNumber;
        var targetDist = route_df.dist[idx];
        profilePlot.plotly.relayout({'xaxis.range': [targetDist - 100, targetDist + 100]});
      }
    });
  }")

# 保存为独立HTML文件
saveWidget(as_widget(last_plot()), "plot.html", selfcontained = TRUE)

注意事项 #

• 若不需要地形底图的色彩映射，可调整scale_fill_viridis_c的参数或替换配色
• 同步滑块的范围（±100m）可根据路径实际长度修改
• 确保route为单条线要素，若为多条需先通过route <- route[1,]筛选单条

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