我来帮你搞定这个问题！要沿着起点到已知点的轨迹方向，算出距离起点250海里的新坐标，我们可以用两种方式实现——要么用R的地理空间包快速搞定，要么手动写球面三角公式来理解底层逻辑。下面分别说明：

geosphere是R里专门处理球面地理计算的包，能帮我们省去复杂的三角计算。步骤很简单：

1. 先安装并加载包（如果还没装的话）
2. 定义起点和已知点的坐标（注意：这个包的函数默认是先经度，后纬度的顺序，别搞反！）
3. 计算起点到已知点的方位角
4. 从起点出发，沿着这个方位角走250海里，算出目标点坐标

代码示例：

# 安装并加载geosphere包（如果未安装）
if (!require(geosphere)) {
  install.packages("geosphere")
  library(geosphere)
}

# 定义起点和已知点坐标（lon在前，lat在后）
start_point <- c(lon_0 = 7.017196, lat_0 = 4.842816)
known_point <- c(lon_1 = 8.099835, lat_1 = 4.108957)

# 计算起点到已知点的方位角（单位：度）
bearing_angle <- bearing(start_point, known_point)

# 计算目标点：250海里转换为米（1海里=1852米）
target_distance_m <- 250 * 1852
target_point <- destPoint(start_point, bearing_angle, target_distance_m)

# 格式化输出结果
colnames(target_point) <- c("经度(lon)", "纬度(lat)")
print(target_point)

如果你想搞清楚底层逻辑，可以用球面三角学的公式自己实现。核心思路是：先算出起点到已知点的方位角，再用这个方位角和目标距离，计算新的经纬度。

核心公式（弧度制）： #

• 目标纬度：φ2 = arcsin( sinφ1 * cos(d/R) + cosφ1 * sin(d/R) * cosθ )
• 目标经度：λ2 = λ1 + arctan2( sinθ * sin(d/R) * cosφ1, cos(d/R) - sinφ1 * sinφ2 )

其中：

• φ1, λ1：起点的纬度、经度（转成弧度）
• d：目标距离（250海里转成米）
• R：地球半径（取赤道半径6378137米）
• θ：起点到已知点的方位角（弧度）

代码示例：

# 定义常量
earth_radius_m <- 6378137  # 地球赤道半径（米）
nm_to_m <- 1852            # 1海里=1852米

# 起点和已知点的原始坐标（lat在前，lon在后）
lat0 <- 4.842816
lon0 <- 7.017196
lat1 <- 4.108957
lon1 <- 8.099835

# 把经纬度转换为弧度（三角函数需要弧度输入）
lat0_rad <- lat0 * pi / 180
lon0_rad <- lon0 * pi / 180
lat1_rad <- lat1 * pi / 180
lon1_rad <- lon1 * pi / 180

# 计算起点到已知点的方位角（弧度）
delta_lon <- lon1_rad - lon0_rad
x <- sin(delta_lon) * cos(lat1_rad)
y <- cos(lat0_rad) * sin(lat1_rad) - sin(lat0_rad) * cos(lat1_rad) * cos(delta_lon)
bearing_rad <- atan2(x, y)

# 计算目标距离（250海里转成米）
target_d_m <- 250 * nm_to_m
d_over_r <- target_d_m / earth_radius_m  # 距离与地球半径的比值（弧度）

# 计算目标点的纬度（弧度）
lat2_rad <- asin( sin(lat0_rad) * cos(d_over_r) + cos(lat0_rad) * sin(d_over_r) * cos(bearing_rad) )

# 计算目标点的经度（弧度）
lon2_rad <- lon0_rad + atan2( sin(bearing_rad) * sin(d_over_r) * cos(lat0_rad), cos(d_over_r) - sin(lat0_rad) * sin(lat2_rad) )

# 把弧度转回度，保留6位小数
lat2 <- round(lat2_rad * 180 / pi, 6)
lon2 <- round(lon2_rad * 180 / pi, 6)

# 输出结果
cat("目标点坐标：\n纬度(lat) =", lat2, "\n经度(lon) =", lon2, "\n")

两种方法算出的结果应该几乎一致，geosphere包的实现更健壮，能处理接近极点、跨日界线等边缘情况，所以日常开发更推荐用第一种方法。

内容的提问来源于stack exchange，提问作者Alex Germain