如何在R语言sf包缓冲区合并操作中消除交叉线
解决缓冲区合并后绘图出现交叉线的问题
在基于WorldClim的bio1栅格生成随机点、创建大范围缓冲区并合并后,转回经纬度坐标系绘图时出现了跨地图的交叉线,这是跨日界线的多边形在经纬度投影下的典型显示问题。
原问题代码
library(raster) library(sf) bioc1 <- getData('worldclim', var='bio', res=10) # 加载WorldClim栅格 bio1 <- bioc1[[1]] # 选择bio1图层 plot(bio1) # 查看栅格 ebio <- extent(c(-208.341377620,194.932040093,-50.879788320,41.18930971)) # 点的范围 set.seed(120) df <- sampleRandom(x=bio1, size=10000, na.rm=TRUE, ext=ebio, xy=TRUE) # 生成随机点 df <- data.frame(df) df.sf4236 <- sf::st_as_sf(df, coords=c("x", "y"), crs=raster::crs(bio1)) # 转为sf对象 plot(bio1) plot(df.sf4236, add=TRUE) # 投影到Eckert IV坐标系(米单位)以创建缓冲区 eckertIV <- "+proj=eck4 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" df.sf54012 <- sf::st_transform(df.sf4236, crs = eckertIV) # 创建2000000米缓冲区并合并 df.buf <- sf::st_buffer(df.sf54012, dist = 2000000) %>% sf::st_union() %>% sf::st_sf() %>% sf::st_transform(crs = raster::crs(bio1)) # 绘图出现交叉线 plot(bio1) plot(df.sf4236, add=TRUE) plot(df.buf, border = "red", lwd = 3, add = TRUE)
问题原因
交叉线是因为跨180°日界线的多边形在经纬度坐标系(WGS84)下显示时,球面转平面的投影特性会导致多边形被错误地拉成横跨整个地图的线条。Eckert IV是伪圆柱投影,全球范围的缓冲区在转回经纬度坐标系时,跨日界线的部分无法正确适配平面的经纬度范围。
解决方法
通过两步处理修复:
- 在投影坐标系下裁剪缓冲区,只保留目标范围内的部分,减少跨日界线的无效区域;
- 使用
st_wrap_dateline()将跨日界线的多边形分割为多个不跨界线的子多边形,使其在经纬度坐标系下正常显示。
修改后的完整代码
library(raster) library(sf) # 加载WorldClim数据 bioc1 <- getData('worldclim', var='bio', res=10) bio1 <- bioc1[[1]] plot(bio1) # 定义点的范围 ebio <- extent(c(-208.341377620,194.932040093,-50.879788320,41.18930971)) # 生成随机点 set.seed(120) df <- sampleRandom(x=bio1, size=10000, na.rm=TRUE, ext=ebio, xy=TRUE) df <- data.frame(df) # 转为sf对象(WGS84坐标系) df.sf4236 <- sf::st_as_sf(df, coords=c("x", "y"), crs=raster::crs(bio1)) plot(bio1) plot(df.sf4236, add=TRUE) # 投影到Eckert IV坐标系(米单位) eckertIV <- "+proj=eck4 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" df.sf54012 <- sf::st_transform(df.sf4236, crs = eckertIV) # 创建缓冲区并合并 df.buf <- sf::st_buffer(df.sf54012, dist = 2000000) %>% sf::st_union() %>% sf::st_sf() # 1. 在Eckert IV坐标系下,用目标范围裁剪缓冲区 ebio_sf <- st_as_sfc(ebio, crs = raster::crs(bio1)) %>% st_transform(crs = eckertIV) df.buf_clipped <- st_intersection(df.buf, ebio_sf) # 2. 转回WGS84并处理跨日界线问题 df.buf_fixed <- df.buf_clipped %>% st_transform(crs = raster::crs(bio1)) %>% st_wrap_dateline(options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180")) # 绘图验证 plot(bio1) plot(df.sf4236, add=TRUE) plot(df.buf_fixed, border = "red", lwd = 3, add = TRUE)
关键修改说明
- 新增
ebio_sf:将原范围转为sf对象并投影到Eckert IV,用它裁剪缓冲区,去除超出目标区域的部分; st_intersection():裁剪后只保留需要的区域,避免不必要的跨日界线多边形;st_wrap_dateline():自动分割跨日界线的多边形,让其在经纬度地图上正确显示,不再出现交叉线。
内容的提问来源于stack exchange,提问作者Ridwan Shittu
相关产品推荐
相关产品推荐

