R中的缓冲区(geo)空间点与gbuffer

Question

我正在尝试缓冲半径为100km的数据集中的点.我正在使用gBuffer包中的功能rgeos.这是我到目前为止所拥有的:

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在我看来是一个自然的选择,但我可能错了.

任何/所有帮助将不胜感激.

Answer 1

通过更多的挖掘,事实证明使用"投影"坐标参考系统非常简单

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

STATSCAN使用的EPSG3347(适用于加拿大地址)使用lambert conformal conic投影.请注意,NAD83是不合适的:它是"地理",而不是"预计"的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" ) 

Answer 2

rgeos::gBuffer()正如@MichaelChirico 指出的那样，应谨慎应用预测要使用的数据。我不是大地测量学方面的专家，但据我从这篇 ESRI 文章（了解测地线缓冲）中了解到，投影然后应用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 包(v0.2.1)于 2019-02-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 包(v0.2.1)于 2019-02-11 创建}

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