r - 在点 p1 周围创建几个点，到该点具有指定的距离。算法错误

Question

我尝试在点p1周围创建几个点，该点的距离为 km_distance。不知何故， km_distance有奇怪的行为，不能按预期工作。我在某个地方有错误。

  p1<-c(11.829831, 48.110075)
  km_distance<-13
  LatDec <- p1[2]
  LonDec <- p1[1]
  
  b <- geosphere::bearing(p1, p2,a = 6378137, f = 1/298.257223563)
  Km <- km_distance
  
  ER <- 6378 #Mean Earth radius in kilometers.
  
  if(b < 0) {
    AngDeg <- seq(b + span_degree,b - span_degree)
  }else
  {
    AngDeg <- seq(b-span_degree,b+span_degree)
  }
  #AngDeg <- seq(b+100,b-100) #angles in degrees 
  Lat1Rad <- LatDec*pi/180#Latitude of the center of the circle in radians
  Lon1Rad <- LonDec*pi/180#Longitude of the center of the circle in radians
  AngRad <- AngDeg*pi/180#angles in radians
  Lat2Rad <- asin(sin(Lat1Rad)*cos(Km/ER) + cos(Lat1Rad)*sin(Km/ER)*cos(AngRad)) #Latitude of each point of the circle rearding to angle in radians
  Lon2Rad <- Lon1Rad + atan2(sin(AngRad)*sin(Km/ER)*cos(Lat1Rad),cos(Km/ER) - sin(Lat1Rad)*sin(Lat2Rad))#Longitude of each point of the circle rearding to angle in radians
  Lat2Deg <- Lat2Rad*180/pi#Latitude of each point of the circle rearding to angle in degrees (conversion of radians to degrees deg = rad*(180/pi) )
  Lon2Deg <- Lon2Rad*180/pi#Longitude of each point of the circle rearding to angle in degrees (conversion of radians to degrees deg = rad*(180/pi) )
  
  point_on_circle <- cbind(Lon2Deg,Lat2Deg)
  point_on_circle<-as.data.frame(point_on_circle)
  point_on_circle <- point_on_circle[seq(1, nrow(point_on_circle), 10), ]

如果我跑

apply(point_on_circle, 1,function(x)   geodist::geodist_vec(p1[1], p2[2], as.numeric(x[1]), as.numeric(x[2]), paired = TRUE, measure = "geodesic"))

我预计到点的距离约为 13 公里，但我的距离要大得多。

标签： rgeolocationdata.tabletidyversegeospatial

Answer 1

像这样的方法可以工作..

p1 <- c(11.829831, 48.110075)
km_distance <- 13 * 1000  # <--!! distance in metres!!
n = 10 #how many points to plot?

#set centre point
start <- matrix( data = p1, ncol = 2, dimnames = list( "", c("lon", "lat")  ) )
#pre-allocate
L <- vector( mode = "list", length = 1+n )
#and initialise centre
L[[1]] <- start

library( geosphere )
set.seed(123) #for reproduction only, remove in production!
for ( i in 2:(n+1) ) {
  L[[i]] <- geosphere::destPoint( p = p1, #start = centre point p1
                                  b = runif(1, 0, 360), #random bearing bewteen 0 and 360 degrees
                                  d = km_distance #distance from p1
  )
}
#transform to a data.table
library( data.table )
ans <- data.table::rbindlist( lapply( L, as.data.table ), use.names = TRUE, fill = TRUE)
#         lon      lat
# 1: 11.82983 48.11008
# 2: 11.96358 48.08844
# 3: 11.91312 48.02324
# 4: 11.82692 48.11503
# 5: 11.80252 48.00736
# 6: 11.80454 48.05945
# 7: 11.80861 48.16114
# 8: 11.74009 48.08061
# 9: 11.92464 48.19400
#10: 11.83717 48.11020
#11: 11.97674 48.05750

#visual confirmation
library(leaflet)
leaflet( data = ans) %>% addTiles() %>%
  addCircleMarkers( lng = ~lon, lat = ~lat )