Python计算EPSG:4326坐标系下LineString米制长度与R结果不符的问题排查

你遇到的问题核心是坐标系处理逻辑的差异：

为什么R的结果是对的？ #

R的sf包非常智能：当你处理的是EPSG:4326这种地理坐标系（基于椭球体的经纬度，不是平面投影）时，st_length()会自动切换到测地线计算模式——它会基于WGS84椭球计算两点间的真实地表距离，所以得到的311米是准确的直线距离。

为什么Python转成EPSG:3857会得到错误结果？ #

EPSG:3857（Web墨卡托）是为网页地图设计的投影，它能保持方向不变，但代价是高纬度地区的距离会被严重拉伸。你这里的坐标位于北纬~60°，墨卡托投影的距离缩放因子大约是1/cos(59.94°)≈2，所以计算出来的平面长度刚好是真实距离的两倍，也就是621米——这是投影变形导致的，不是你的坐标顺序写错了。

正确的Python实现方式 #

要得到和R一致的准确测地线距离，不需要转成投影坐标系，直接计算椭球上的距离即可：

方法1：用pyproj的Geod类（兼容所有shapely版本） #

from pyproj import Geod
from shapely.geometry import LineString

# 你的原始坐标，转换为(lon, lat)的顺序（对应EPSG:4326的x,y）
latlon = [[59.94661244955431, 10.72052329576658], [59.948009419166226, 10.725362000881619]]
xy = [(lon, lat) for lat, lon in latlon]
line = LineString(xy)

# 使用WGS84椭球初始化Geod对象（和EPSG:4326的基准一致）
geod = Geod(ellps="WGS84")
# 计算LineString的测地线长度，单位是米
length_meters = geod.geometry_length(line)
print(length_meters)  # 输出约311.017米，和R结果完全匹配

方法2：用shapely 2.0+的length_crs方法（更简洁） #

如果你用的是shapely 2.0及以上版本，可以直接调用几何对象的length_crs方法，指定地理坐标系即可：

from shapely.geometry import LineString
from pyproj import CRS

latlon = [[59.94661244955431, 10.72052329576658], [59.948009419166226, 10.725362000881619]]
xy = [(lon, lat) for lat, lon in latlon]
line = LineString(xy)

# 指定地理坐标系
wgs84 = CRS("EPSG:4326")
# 直接计算测地线长度
length_meters = line.length_crs(wgs84)
print(length_meters)  # 同样得到约311米

内容的提问来源于stack exchange，提问作者Kjetil Haukås