基于经纬度计算船舶到海岸最短距离的代码问题排查

核心问题1：坐标顺序颠倒 #

Shapely的Point构造参数是**(经度, 纬度)，但你传入的(36.20972222, 125.7061111)是(纬度, 经度)**的顺序。这会导致目标点被错误定位，甚至刚好落在海岸线数据的平面投影范围内，从而计算出最近点与目标点重合、距离为0的错误结果。

核心问题2：球面坐标直接用平面算法计算 #

Shapely的nearest_points是基于笛卡尔平面的计算逻辑，而经纬度属于球面地理坐标，直接用平面算法计算最近点会产生巨大误差，尤其是在高纬度地区。正确的做法是先将地理坐标转换为适合的投影坐标系（如UTM，该坐标系能将球面坐标转换为平面坐标，且局部区域误差极小），再进行最近点计算。

修正后的代码 #

from shapely.geometry import Point
from shapely.ops import nearest_points
from geopy.distance import geodesic
import geopandas as gpd

# 读取海岸线数据，并确认坐标系（默认是WGS84，EPSG:4326）
coastline_data = gpd.read_file('./ne_10m_coastline.shp')
# 转换为UTM投影（自动匹配目标点对应的UTM带）
target_lat, target_lon = 36.20972222, 125.7061111
utm_crs = gpd.GeoSeries([Point(target_lon, target_lat)]).estimate_utm_crs()
coastline_proj = coastline_data.to_crs(utm_crs)
coastline_union = coastline_proj.geometry.unary_union

# 创建正确的目标点（经度在前，纬度在后），并转换为UTM投影
target_point = Point(target_lon, target_lat)
target_proj = gpd.GeoSeries([target_point], crs="EPSG:4326").to_crs(utm_crs).iloc[0]

# 在投影坐标系下计算最近点
nearest_proj = nearest_points(target_proj, coastline_union)[0]
# 将最近点转回WGS84地理坐标系
nearest_geo = gpd.GeoSeries([nearest_proj], crs=utm_crs).to_crs("EPSG:4326").iloc[0]

# 用geodesic计算球面距离（参数为(纬度, 经度)）
distance = geodesic(
    (target_point.y, target_point.x),
    (nearest_geo.y, nearest_geo.x)
).kilometers

print(f"海岸最近点坐标：{nearest_geo.y:.6f}, {nearest_geo.x:.6f}")
print(f"到海岸的最短距离：{distance:.2f} 公里")

额外说明 #

• 如果你坚持用地理坐标直接计算，可以使用geopandas的distance方法（该方法会自动处理球面距离计算），但效率远低于投影坐标系下的平面计算：

  # 直接用GeoSeries计算球面距离
  target_gs = gpd.GeoSeries([target_point], crs="EPSG:4326")
  min_distance = coastline_data.distance(target_gs.iloc[0]).min()
  # 这里的min_distance是度数，需要转换为公里（1度≈111公里，仅近似，不如geodesic准确）
  distance_km = min_distance * 111
  

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