正如@MichaelChirico所指出的那样,使用rgeos::gBuffer()投影数据应谨慎处理。我不是大地测量学专家,但从这篇ESRI文章中所理解的(
Understanding Geodesic Buffering),投影然后应用gBuffer实际上意味着产生欧几里得缓冲区,而不是大圆缓冲区。欧几里得缓冲区受到投影坐标系统引入的扭曲的影响。如果您的分析涉及跨越大面积的广泛缓冲区,特别是在更广泛的纬度范围内(我认为加拿大是一个很好的候选国家),这些扭曲可能是值得担忧的事情。
我曾经遇到过同样的问题,并将我的问题针对gis.stackexchange -
Euclidean and Geodesic Buffering in R。我认为当时我提出的R代码以及给出的答案对于这个问题也是相关的。
主要思路是利用
geosphere::destPoint()
。更多细节和更快的替代方法,请参见上面提到的gis.stackexchange链接。这是我之前在你的两个点上尝试过的方法:
library(geosphere)
library(sp)
pts <- data.frame(lon = c(-53.20198, -52.81218),
lat = c(47.05564, 47.31741))
pts
make_GeodesicBuffer <- function(pts, width) {
dg <- seq(from = 0, to = 360, by = 5)
buff.XY <- geosphere::destPoint(p = pts,
b = rep(dg, each = length(pts)),
d = width)
buff.XY <- as.data.frame(buff.XY)
id <- rep(1:dim(pts)[1], times = length(dg))
lst <- split(buff.XY, id)
poly <- lapply(lst, sp::Polygon, hole = FALSE)
polys <- lapply(list(poly), sp::Polygons, ID = NA)
spolys <- sp::SpatialPolygons(Srl = polys,
proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))
spolys <- sp::disaggregate(spolys)
return(spolys)
}
pts_buf_100km <- make_GeodesicBuffer(as.matrix(pts), width = 100*10^3)
library(plotKML)
kml(pts_buf_100km, file.name = "pts_buf_100km.kml")
这段内容是由reprex package (v0.2.1)在2019年2月11日创建的。
![enter image description here](https://istack.dev59.com/X1MQZ.webp)
为了玩一下,我将该函数封装在一个包中 - geobuffer
以下是一个示例:
devtools::install_github("valentinitnelav/geobuffer")
library(geobuffer)
pts <- data.frame(lon = c(-53.20198, -52.81218),
lat = c(47.05564, 47.31741))
pts_buf_100km <- geobuffer_pts(xy = pts, dist_m = 100*10^3)
这篇文章是由reprex package (v0.2.1)于2019-02-11创建的。
其他人可能会提出更好的解决方案,但目前为止,这对我的问题很有效,希望也能解决其他人的问题。