如何使用R+Leaflet高亮显示线要素上两点之间的指定线段

要实现线要素上两点之间原有路径的高亮，核心是截取两点之间的线要素子段，再将子段用差异化样式渲染即可，不需要依赖第三方地图服务，完全基于sf生态即可实现。

实现步骤 #

1. 安装加载所需依赖包，除已使用的包外，需要用到lwgeom包提供的线要素截取功能
2. 计算两个点在线要素上的归一化位置（距离起点长度占线总长度的比例）
3. 根据两个位置比例截取对应子线段
4. 将子线段用高亮样式添加到Leaflet地图中

完整可运行代码 #

library(tidyverse)
library(sf) 
library(leaflet)
# 首次运行需安装依赖包
# install.packages("lwgeom")
library(lwgeom)

# 示例线要素（铁路线）
tr <- structure(
  list(
    SHAPE_Length = 22210.9047689, 
    SHAPE = structure(list(
      structure(c(-3.447431, -3.426467, -3.3717497, -3.3504529, 
                  -3.34384, -3.3063404, -3.2875194, -3.2388496, -3.1727098, 
                  -3.136474, -3.1031447, 58.51414, 58.5140392, 58.5066418, 
                  58.5013098, 58.4950488, 58.4901998, 58.4837109, 58.4745371, 
                  58.4552422, 58.4512823, 58.4421872), .Dim = c(11L, 2L),
                class = c("XY", "LINESTRING", "sfg"))), 
      n_empty = 0L, crs = structure(list(
        input = "+init=epsg:4326", 
        wkt = "GEOGCRS[\"WGS 84\",\n    
        DATUM[\"World Geodetic System 1984\",\n        
        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            
        LENGTHUNIT[\"metre\",1]],\n        
        ID[\"EPSG\",6326]],\n    
        PRIMEM[\"Greenwich\",0,\n        
        ANGLEUNIT[\"degree\",0.0174532925199433],\n        
        ID[\"EPSG\",8901]],\n    
        CS[ellipsoidal,2],\n        
        AXIS[\"longitude\",east,\n            
        ORDER[1],\n            
        ANGLEUNIT[\"degree\",0.0174532925199433,\n                
        ID[\"EPSG\",9122]]],\n        
        AXIS[\"latitude\",north,\n            
        ORDER[2],\n            
        ANGLEUNIT[\"degree\",0.0174532925199433,\n                
        ID[\"EPSG\",9122]]],\n    
        USAGE[\n        
        SCOPE[\"unknown\"],\n        
        AREA[\"World.\"],\n        
        BBOX[-90,-180,90,180]]]", 
        class = "crs"), 
      class = c("sfc_LINESTRING", "sfc"), 
      precision = 0, 
      bbox = structure(c(xmin = -3.447431, ymin = 58.4421872, 
                         xmax = -3.1031447, ymax = 58.51414), class = "bbox"))), 
  row.names = 1L, class = c("sf", "data.frame"), 
  sf_column = "SHAPE", 
  agr = structure(c(SHAPE_Length = NA_integer_), 
                  .Label = c("constant", "aggregate", "identity"), class = "factor"))

# 示例线上点位
sg <- structure(list(
  geometry = structure(list(structure(c(-3.12213008366773, 58.4467460240336), 
                                      class = c("XY", "POINT", "sfg")), 
                            structure(c(-3.28773441769833, 58.483750111291), 
                                      class = c("XY", "POINT", "sfg"))), 
                       class = c("sfc_POINT", "sfc"), 
                       precision = 0, 
                       bbox = structure(c(xmin = -3.28773441769833, ymin = 58.4467460240336, 
                                          xmax = -3.12213008366773, ymax = 58.483750111291), 
                                        class = "bbox"), 
                       crs = structure(list(input = "+init=epsg:4326",
                                            wkt = "GEOGCRS[\"WGS 84\",\n    
                                            DATUM[\"World Geodetic System 1984\",\n        
                                            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            
                                            LENGTHUNIT[\"metre\",1]],\n        
                                            ID[\"EPSG\",6326]],\n    
                                            PRIMEM[\"Greenwich\",0,\n        
                                            ANGLEUNIT[\"degree\",0.0174532925199433],\n        
                                            ID[\"EPSG\",8901]],\n    
                                            CS[ellipsoidal,2],\n        
                                            AXIS[\"longitude\",east,\n            
                                            ORDER[1],\n            
                                            ANGLEUNIT[\"degree\",0.0174532925199433,\n                
                                            ID[\"EPSG\",9122]]],\n        
                                            AXIS[\"latitude\",north,\n            
                                            ORDER[2],\n            
                                            ANGLEUNIT[\"degree\",0.0174532925199433,\n                
                                            ID[\"EPSG\",9122]]],\n    
                                            USAGE[\n        
                                            SCOPE[\"unknown\"],\n        
                                            AREA[\"World.\"],\n        
                                            BBOX[-90,-180,90,180]]]", 
                                       class = "crs"), n_empty = 0L)), 
  row.names = 1:2, class = c("sf","data.frame"), 
  sf_column = "geometry", 
  agr = structure(integer(0), .Label = c("constant","aggregate", "identity"), 
                  class = "factor", .Names = character(0)))

# --------------- 核心截取逻辑 ---------------
# 提取线要素
tr_line <- tr$SHAPE[[1]]
# 计算两个点在线上的归一化距离（0为线起点，1为线终点）
dist_vec <- st_distance(sg$geometry, tr_line, along = TRUE) / as.numeric(st_length(tr_line))
start_ratio <- min(dist_vec)
end_ratio <- max(dist_vec)
# 截取两点之间的子线段
highlight_segment <- st_linesubstring(tr_line, start_ratio, end_ratio) %>% 
  st_sfc(crs = st_crs(tr))

# 绘制地图
leaflet() %>%
  addTiles() %>%
  # 原有全线要素，蓝色细线条
  addPolylines(data = tr, color = "#2c7fb8", weight = 3) %>%
  # 高亮子线段，红色粗线条
  addPolylines(data = highlight_segment, color = "#de2d26", weight = 6, opacity = 0.8) %>%
  addMarkers(data = sg)

注意事项 #

如果你的点位存在微小偏移、没有完全落在线上，可以先使用sg_snap <- st_snap(sg, tr, tolerance = 0.0001)（tolerance单位和坐标系一致，此处为度，可根据实际点位偏移情况调整）将点吸附到线上后再进行后续计算，避免截取位置偏差。

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