如何用R的sf包替代弃用工具实现指定国家外世界遮罩并在Leaflet展示

之前使用sp、rgdal、rgeos包实现了Leaflet地图中「除指定国家（如德国）外，其余区域遮罩变白」的效果，但这些包已被弃用。现需要用sf（或结合terra）复现该功能，且确保适配leaflet::addPolygons()。

旧实现的遮罩生成函数：

mask_world_except_specified_region <- function(path, shp_filename) {
  specified_area <- rgdal::readOGR(path, shp_filename)
  world <- raster::rasterToPolygons(raster::raster(ncol = 1, nrow = 1, crs = sp::proj4string(specified_area)))
  specified_area_mask <- rgeos::gDifference(world, specified_area)
  tmp_id_df <- data.frame(ID = "1")
  rownames(tmp_id_df ) <- names(specified_area_mask)
  sp::SpatialPolygonsDataFrame(specified_area_mask, tmp_id_df)
}

基于sf的替代实现 #

步骤说明： #

• 使用sf::st_read()读取指定区域的矢量数据，替代rgdal::readOGR()
• 创建覆盖全球的多边形（用sf原生构造，替代raster::rasterToPolygons()）
• 使用sf::st_difference()计算全球多边形与指定区域的差集，替代rgeos::gDifference()
• 为结果添加属性列，确保适配Leaflet的addPolygons()

新遮罩生成函数 #

library(sf)

mask_world_except_specified_region_sf <- function(path, shp_filename) {
  # 读取指定区域的sf对象，静默输出
  specified_area <- st_read(path, shp_filename, quiet = TRUE)
  
  # 创建与指定区域同CRS的全球边界多边形
  world_bbox <- st_bbox(c(xmin = -180, ymin = -90, xmax = 180, ymax = 90), crs = st_crs(specified_area))
  world_poly <- st_as_sfc(world_bbox)
  
  # 计算差集：全球范围减去指定区域，得到遮罩区域
  mask_area <- st_difference(world_poly, specified_area)
  
  # 添加属性列，满足leaflet对带属性矢量对象的要求
  mask_sf <- st_sf(ID = "1", geometry = mask_area)
  
  return(mask_sf)
}

完整地图绘制代码 #

library(leaflet)

create_leaflet_map_germany <- function(df, mask_sf) {
  leaflet(data = df, options = leafletOptions(minZoom = 6, dragging = TRUE)) %>%
    addProviderTiles(provider = "CartoDB.PositronNoLabels", group = "Carto") %>%
    addProviderTiles(provider = "Esri.WorldStreetMap", group = "Esri") %>%
    addTiles(group = "OSM") %>%
    # 传入sf格式的遮罩层
    addPolygons(data = mask_sf, fillColor = "white", fillOpacity = 1.0, color = "black", weight = 1) %>%
    setView(lng = 10.3608414, lat = 51.1593825, zoom = 6) %>%
    addCircleMarkers(~lon, ~lat, color = '#195365', radius = 7, fillColor = '#195365', popup = ~Popup) %>%
    addResetMapButton() %>%
    addLayersControl(baseGroups = c("Carto", "Esri", "OSM"), position = "topright")
}

# 示例点位数据
df <- data.frame(
  lon = c(10, 11.34),
  lat = c(53.34, 48.08),
  Popup = c("Moin", "Servus")
)

# 生成德国区域的遮罩sf对象
germany_mask_sf <- mask_world_except_specified_region_sf("my/path", "gadm_36_DEU_1")

# 绘制地图
create_leaflet_map_germany(df, germany_mask_sf)

注意事项 #

• 若指定区域的矢量数据CRS不是WGS84（EPSG:4326），需用st_transform(specified_area, 4326)转换后再计算差集，确保与Leaflet兼容
• 处理复杂多边形时，若出现拓扑错误，可先用st_make_valid(specified_area)修复几何对象

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