如何在R语言中从GPX文件生成独立的地形与路径STL文件？

你的问题核心在于：调用save_3dprint时会保存当前整个3D场景，而你在渲染路线前已经绘制了地形，所以导出的soRibbon.stl包含了地形+路线的组合，而非单独的路线。另外，render_path的参数传递写法错误（用了赋值语句lat <- ...），这会导致参数传递异常。

解决步骤 #

1. 分离地形与路线的生成流程：先生成并保存地形STL，然后清空当前3D场景，单独构建路线的3D模型再保存。
2. 修正render_path的参数传递：直接传入向量，而非赋值语句。
3. 确保路线是地形上方的带状结构：通过offset参数控制路线相对于地形的高度，同时注意zscale的一致性（和地形生成时的zscale保持一致，否则高度比例会错）。

修正后的完整代码 #

library(rayshader)
library(sf)
library(xml2)
library(elevatr)
library(raster)
library(tidyverse)

gpx_text <- '<?xml version="1.0" encoding="utf-8"?>
<gpx xmlns="http://www.topografix.com/GPX/1/1" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.1" xsi:schemaLocation="http://www.topografix.com/GPX/1/1 http://www.topografix.com/GPX/1/1/gpx.xsd" creator="https://mtbproject.com">
<metadata><name>Sope Creek Start</name><author><name/><link href="https://www.mtbproject.com"/></author></metadata><trk><name>Sope Creek Start</name><link href="https://www.mtbproject.com"/></trk><trkseg><trkpt lon="-84.445526" lat="33.90424"><ele>241.35</ele><time>2022-01-08T08:08:21-07:00</time></trkpt><trkpt lon="-84.446101" lat="33.905142"><ele>241.63</ele><time>2022-01-08T08:09:21-07:00</time></trkpt><trkpt lon="-84.447467" lat="33.906283"><ele>241.85</ele><time>2022-01-08T08:10:21-07:00</time></trkpt><trkpt lon="-84.448536" lat="33.907252"><ele>241.9</ele><time>2022-01-08T08:11:21-07:00</time></trkpt><trkpt lon="-84.448841" lat="33.907923"><ele>241.84</ele><time>2022-01-08T08:12:21-07:00</time></trkpt><trkpt lon="-84.448877" lat="33.908989"><ele>241.81</ele><time>2022-01-08T08:13:21-07:00</time></trkpt><trkpt lon="-84.44859" lat="33.910092"><ele>241.84</ele><time>2022-01-08T08:14:21-07:00</time></trkpt><trkpt lon="-84.448697" lat="33.91048"><ele>241.99</ele><time>2022-01-08T08:15:21-07:00</time></trkpt><trkpt lon="-84.449713" lat="33.911665"><ele>242.3</ele><time>2022-01-08T08:16:21-07:00</time></trkpt><trkpt lon="-84.450539" lat="33.911755"><ele>242.87</ele><time>2022-01-08T08:17:21-07:00</time></trkpt><trkpt lon="-84.450701" lat="33.911867"><ele>243.58</ele><time>2022-01-08T08:18:21-07:00</time></trkpt><trkpt lon="-84.450135" lat="33.913231"><ele>243.99</ele><time>2022-01-08T08:19:21-07:00</time></trkpt><trkpt lon="-84.449722" lat="33.913954"><ele>244.05</ele><time>2022-01-08T08:20:21-07:00</time></trkpt><trkpt lon="-84.448868" lat="33.914767"><ele>244.53</ele><time>2022-01-08T08:21:21-07:00</time></trkpt><trkpt lon="-84.44753" lat="33.915557"><ele>246.48</ele><time>2022-01-08T08:22:21-07:00</time></trkpt><trkpt lon="-84.446838" lat="33.915892"><ele>250.05</ele><time>2022-01-08T08:23:21-07:00</time></trkpt><trkpt lon="-84.446344" lat="33.916526"><ele>254.4</ele><time>2022-01-08T08:24:21-07:00</time></trkpt><trkpt lon="-84.446613" lat="33.917927"><ele>258.69</ele><time>2022-01-08T08:25:21-07:00</time></trkpt><trkpt lon="-84.446604" lat="33.918315"><ele>262.3</ele><time>2022-01-08T08:26:21-07:00</time></trkpt><trkpt lon="-84.446919" lat="33.919001"><ele>265.04</ele><time>2022-01-08T08:27:21-07:00</time></trkpt></trkseg></trk></gpx>
'

gpx  <- read_xml(gpx_text, encoding = "UTF-8")
ns <- xml_ns(gpx)

trkpts <- xml_find_all(gpx, ".//d1:trkpt", ns)
trail_df <- tibble(
  lon = as.numeric(xml_attr(trkpts, "lon")),
  lat = as.numeric(xml_attr(trkpts, "lat")),
  ele = as.numeric(xml_text(xml_find_all(trkpts, "d1:ele", ns)))
)

# 处理坐标转换（原代码这部分没问题）
mean_coords <- colMeans(trail_df[, c("lon", "lat")])
lon <- mean_coords["lon"]
lat <- mean_coords["lat"]
zone <- floor((lon + 180) / 6) + 1
epsg_code <- ifelse(lat >= 0, 32600 + zone, 32700 + zone)
trail_sp <- st_as_sf(trail_df %>% select(lon, lat), coords = c("lon", "lat"), crs = 4326)

# --------------------------
# 第一步：生成并保存地形STL
# --------------------------
terrain_raster <- get_elev_raster(locations = trail_sp, z = 13, clip = "bbox")
terrain_matrix <- raster_to_matrix(terrain_raster)

# 渲染地形3D场景
tex <- sphere_shade(terrain_matrix, texture = "imhof1")
ray <- ray_shade(terrain_matrix, sunaltitude = 40, sunangle = 315)
amb <- ambient_shade(terrain_matrix)

shaded <- tex %>%
  add_shadow(ray, 0.5) %>%
  add_shadow(amb, 0.3) %>%
  plot_3d(terrain_matrix, 
          zscale = 10,
          fov = 0,
          theta = 45, phi = 45,
          windowsize = c(800, 800),
          plot_new = TRUE)

# 保存地形STL
save_3dprint(
  "soTerrain.stl",
  maxwidth = 125,
  unit = "mm",
  rotate = TRUE
)

# 清空当前3D场景，为生成路线做准备
rgl::clear3d()

# --------------------------
# 第二步：生成并保存路线STL
# --------------------------
ribbon_height <- 1  # 路线相对于地形的高度（单位：米，因为zscale=10，所以实际3D中高度是ribbon_height*10）
ribbon_thickness <- 3  # 路线宽度

# 修正参数传递：直接传向量，不用赋值语句
render_path(
  lat = trail_df$lat,
  long = trail_df$lon,
  altitude = trail_df$ele,
  extent = terrain_raster,  # 用地形的extent，确保坐标匹配
  zscale = 10,  # 和地形生成时的zscale一致，保证高度比例正确
  offset = ribbon_height,
  linewidth = ribbon_thickness,
  color = "red"
)

# 保存路线STL
save_3dprint("soRibbon.stl", maxwidth = 125, unit = "mm", rotate = TRUE)

关键修改点说明 #

• 清空场景：用rgl::clear3d()在生成路线前清空地形的3D场景，确保路线STL只包含路线模型。
• 修正参数传递：render_path的参数用=而非<-，避免赋值导致的参数错误。
• 匹配extent：路线的extent用terrain_raster而非trail_sp，确保路线坐标和地形的坐标系统完全对齐。
• zscale一致性：地形和路线的zscale保持相同，保证路线的高度偏移和地形的垂直比例一致。

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