У меня есть два кадра данных длинных точек.DatasetA - это набор местоположений.Набор данных B представляет собой набор местоположений в 50 милях (в продольном направлении) от набора данных A.
Я могу построить два набора данных без инцидентов, как показано в коде ниже.Однако, когда я вычисляю расстояние между точками DatasetA и DatasetB (ближайшие точки), я получаю разные расстояния:
> points$Distance2radius
[1] 44.13807 41.92467 39.39219 36.55992 33.44940
У меня возникают проблемы с пониманием, почему они будут разными.Я бы предположил, что использование формулы Haversine в distm будет учитывать любые сферические воздействия на расчет расстояния.
Любая помощь будет оценена.
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