使用 geosphere::distm 计算的两点之间的距离与使用小册子绘制的距离不同答案

【问题标题】：Distances between two points calculated with geosphere::distm are different than distances plotted with leaflet使用 geosphere::distm 计算的两点之间的距离与使用小册子绘制的距离不同
【发布时间】：2019-06-11 12:27:44
【问题描述】：

我有两个经纬度点的数据框。 DatasetA 是一组位置。数据集 B 是距离数据集 A 50 英里（纵向）的一组位置。

我可以毫无意外地绘制两个数据集，如下面的代码所示。但是，当我计算 DatasetA 和 DatasetB 的点（最近点）之间的距离时，我得到不同的距离：

> points$Distance2radius
[1] 44.13807 41.92467 39.39219 36.55992 33.44940

我很难理解为什么这些会有所不同。我假设在 distm 中使用 Haversine 公式可以解释距离计算中的任何球面影响。

任何帮助将不胜感激。

library(leaflet)
library(geosphere)

### Make a dataframe of some test points ###

## Center of the US
dc.lat <- c(38.0000)
dc.long <- c(-97.0000)

## Make the data frame
lats <- c(dc.lat - 10, dc.lat - 5, dc.lat, 5 + dc.lat, 10 + dc.lat)
points <- data.frame(cbind(dc.long, lats))
names(points) <- c("long" , "lat")
coordinates(points) <- ~ long + lat

## The radius we are interested in, in miles
radius <- 50

## Add points that are the radius miles away from the center

points.at.radius <- data.frame(points)
#points$lat <- points$lat + radius/110.54
points.at.radius$long <- points$long + radius / 69.2
coordinates(points.at.radius) <- ~ long + lat

## Get distances with distm
distances <- distm (points, points.at.radius,
                    fun = distHaversine) / 1609

# Find the closest pint and add this distance to the data set
points$Distance2radius <-
  apply(distances , 1, min)

# Plot these points and the points that are 50 miles away
m <- leaflet() %>%
  addTiles() %>%  # Add default OpenStreetMap map tiles
  addCircleMarkers(data = points,
                   points$long,
                   points$lat,
                   color = "red") %>%
  addPolygons(
    data = gBuffer(
      points,
      width =  radius / 69.2,
      joinStyle = "ROUND",
      byid = FALSE
    ),
    color = "gray70",
    group = "IWER area"
  ) %>%
  addMarkers(data = points.at.radius,
             points.at.radius$long,
             points.at.radius$lat) %>%  addPopups(
               points$long,
               points$lat,
               47.597131,
               points$Distance2radius,
               options = popupOptions(closeButton = FALSE)
             ) 

  m  # Print the map

【问题讨论】：

您忘记投影数据或计算未投影的半径距离。

标签： r shiny leaflet

【解决方案1】：

geobuffer 非常简单。

代码

# Load libraries.
library(leaflet)
library(geosphere)
library(geobuffer)

# Set coordinates for center points.
coordinates <- data.frame(lon = -97,
                          lat = seq(28, 48, 5))

# Create SpatialPoints in the WGS84 spatial reference system,
# a.k.a. unprojected lon/lat.
points <- SpatialPoints(coords = coordinates,
                        proj4string = CRS("+init=epsg:4326"))

# Set radius as 50 miles.
radius <- 50 * 1609

# Created SpatialPolygons in the shape of an octagon. 
buffered_points <- geobuffer::geobuffer_pts(xy = points,
                                            dist_m = radius,
                                            step_dg = 360 / 8)

测试

> distm(points@coords[1,], 
      coordinates(buffered_points@polygons[[1]]@Polygons[[1]])) / 1609

     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,]   50   50   50   50   50   50   50   50   50

您会发现不仅 SpatialPolygon #1 如此，其他三个对象也是如此。

结果

leaflet() %>%
    addTiles() %>%
    addCircleMarkers(data = points,
                     color = "red") %>%
    addPolygons(
        data = buffered_points,
        color = "gray70"
    )

【讨论】：