使用gbuffer在R中缓存（地理）空间点

Question

我正在尝试使用100公里的半径缓冲我的数据集中的点。我正在使用rgeos软件包中的gBuffer函数。以下是我目前的代码：

head( sampledf )
#  postalcode      lat       lon       city province
#1     A0A0A0 47.05564 -53.20198     Gander       NL
#4     A0A1C0 47.31741 -52.81218 St. John's       NL

coordinates( sampledf ) <- c( "lon", "lat" )
proj4string( sampledf ) <- CRS( "+proj=longlat +datum=WGS84" )
distInMeters <- 1000
pc100km <- gBuffer( sampledf, width=100*distInMeters, byid=TRUE )

我收到以下警告:

在gBuffer（sampledf，width =100 * distInMeters，byid = TRUE）中： 空间对象未投影; GEOS期望平面坐标

从我的理解和阅读中，我需要更改数据集的坐标参考系统（CRS），特别是将其从“地理”投影为“投影”。 我不确定如何更改。我要补充的是，这些都是加拿大地址。所以NAD83似乎是我选择的自然投影，但我可能错了。

非常感谢任何/所有帮助。

- vathymut

---

这虽然是一件小事，但是NAD83并不是一个投影。它是一个大地水准面（North American Datum 1983），代表了地球形状的模型（椭球体与球体）。 - Foxhound013

Answer 1

经过更深入的挖掘，发现使用“投影”坐标参考系统就像这样简单：

# To get Statscan CRS, see here:
# http://spatialreference.org/ref/epsg/3347/
pc <- spTransform( sampledf, CRS( "+init=epsg:3347" ) ) 

EPSG3347，被STATSCAN使用（适用于加拿大地址），使用兰伯特等积圆锥投影。请注意，NAD83不合适：它是一个“地理”的CRS，而不是一个“投影”的CRS。为了缓冲这些点

pc100km <- gBuffer( pc, width=100*distm, byid=TRUE )
# Add data, and write to shapefile
pc100km <- SpatialPolygonsDataFrame( pc100km, data=pc100km@data )
writeOGR( pc100km, "pc100km", "pc100km", driver="ESRI Shapefile" ) 

- vathymut

---

重要提示：除非您的对象在加拿大，否则请勿使用 epsg:3347。请注意保持原意并使翻译更加通俗易懂，但不要添加解释或返回其他内容。 - MichaelChirico

如果返回的是 sf "data.frame"，则应使用 st_transform 而不是 spTransform。 - Mox

Answer 2

正如@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
#>         lon      lat
#> 1 -53.20198 47.05564
#> 2 -52.81218 47.31741

make_GeodesicBuffer <- function(pts, width) {

  # A) Construct buffers as points at given distance and bearing ---------------

  dg <- seq(from = 0, to = 360, by = 5)

  # Construct equidistant points defining circle shapes (the "buffer points")
  buff.XY <- geosphere::destPoint(p = pts, 
                                  b = rep(dg, each = length(pts)), 
                                  d = width)

  # B) Make SpatialPolygons -------------------------------------------------

  # Group (split) "buffer points" by id
  buff.XY <- as.data.frame(buff.XY)
  id  <- rep(1:dim(pts)[1], times = length(dg))
  lst <- split(buff.XY, id)

  # Make SpatialPolygons out of the list of coordinates
  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"))
  # Disaggregate (split in unique polygons)
  spolys <- sp::disaggregate(spolys)
  return(spolys)
}

pts_buf_100km <- make_GeodesicBuffer(as.matrix(pts), width = 100*10^3)

# Make a kml file and check the results on Google Earth
library(plotKML)
#> plotKML version 0.5-9 (2019-01-04)
#> URL: http://plotkml.r-forge.r-project.org/
kml(pts_buf_100km, file.name = "pts_buf_100km.kml")
#> KML file opened for writing...
#> Writing to KML...
#> Closing  pts_buf_100km.kml

^{这段内容是由reprex package (v0.2.1)在2019年2月11日创建的。}

---

为了玩一下，我将该函数封装在一个包中 - geobuffer

以下是一个示例：

# install.packages("devtools") # if you do not have devtools, then install it
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创建的。}

其他人可能会提出更好的解决方案，但目前为止，这对我的问题很有效，希望也能解决其他人的问题。

- Valentin_Ștefan

---

我希望你能用R包绘制这个图。我试着复现它，但在最后才意识到你已经将其导出并在其他地方绘制了。我想知道是哪个神奇的包绘制了那种地图！尽管如此，这仍然是一个很好的例子！ - geneorama

嗨@geneorama，是的，为了举例说明，我只是用plotKML::kml制作了kml文件，然后在Google Earth中打开并截屏。或者，您可以尝试mapview，它基于leaflet提供交互式地图。不确定这是否符合您的需求。 - Valentin_Ștefan