在点 p1 周围创建多个与该点具有指定距离的点算法错误

如何解决在点 p1 周围创建多个与该点具有指定距离的点算法错误

我尝试在 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 公里的距离，但我得到的距离要大得多。

解决方法

这样的方法可行..

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,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 )