为什么我的代码计算的变换坐标之间的距离不正确?
Why is my code calculating an incorrect distance between transformed coordinates?
我有一个很大的点文件,我正在尝试找出这些点与另一组点之间的距离。最初我使用 geopandas 的 to_crs
函数来转换 crs,这样当我做 df.distance(point)
时我可以获得以米为单位的准确距离测量。但是,由于文件很大,光是转换文件的crs就花了很长时间。 code 运行 2个小时了,还是没转换完。因此,我改用了这段代码。
inProj = Proj(init='epsg:4326')
outProj = Proj(init='epsg:4808')
for index, row in demand_indo_gdf.iterrows():
o = Point(row['origin_longitude'], row['origin_latitude'])
o_proj = Point(transform(inProj, outProj, o.x, o.y))
for i, r in bus_indo_gdf.iterrows():
stop = r['geometry']
stop_proj = Point(transform(inProj, outProj, stop.x, stop.y))
print ('distance:', o_proj.distance(stop_proj), '\n\n')
我认为单独转换 crs 并执行我的分析可能会更快。对于这组点数:
o = (106.901024 -6.229162)
stop = (106.804 -6.21861)
我将这个 EPSG 4326 坐标转换为局部投影 EPSG 4808,得到了这个:
o_proj = (0.09183386384156803 -6.229330112968891)
stop_proj = (-0.005201753272169649 -6.218776788266844)
这给出了 0.09760780527657992 的距离度量。 Google 地图为我提供了坐标 o
和 stop
的距离测量值,为 10.79 公里。看起来我的代码的距离度量给出的答案比实际距离小 10^-3 倍。为什么?我的代码正确吗?
您使用的坐标是以度为单位,而不是以米为单位。
你用的Points
class可能不关心这个,计算它们之间的笛卡尔距离
请使用Haversine equation, e.g. the one from here:
from math import radians, cos, sin, asin, sqrt
def haversine(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
"""
# convert decimal degrees to radians
lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
c = 2 * asin(sqrt(a))
r = 6371 # Radius of earth in kilometers. Use 3956 for miles
return c * r
那么你的代码returns正确结果:
from pyproj import Proj, transform
inProj = Proj(init='epsg:4326')
outProj = Proj(init='epsg:4808')
o = (106.901024, -6.229162)
stop = (106.804, -6.21861)
o_proj = transform(inProj, outProj, o[0], o[1])
stop_proj = transform(inProj, outProj, stop[0], stop[1])
print(haversine(o_proj[0], o_proj[1], stop_proj[0], stop_proj[1]))
我有一个很大的点文件,我正在尝试找出这些点与另一组点之间的距离。最初我使用 geopandas 的 to_crs
函数来转换 crs,这样当我做 df.distance(point)
时我可以获得以米为单位的准确距离测量。但是,由于文件很大,光是转换文件的crs就花了很长时间。 code 运行 2个小时了,还是没转换完。因此,我改用了这段代码。
inProj = Proj(init='epsg:4326')
outProj = Proj(init='epsg:4808')
for index, row in demand_indo_gdf.iterrows():
o = Point(row['origin_longitude'], row['origin_latitude'])
o_proj = Point(transform(inProj, outProj, o.x, o.y))
for i, r in bus_indo_gdf.iterrows():
stop = r['geometry']
stop_proj = Point(transform(inProj, outProj, stop.x, stop.y))
print ('distance:', o_proj.distance(stop_proj), '\n\n')
我认为单独转换 crs 并执行我的分析可能会更快。对于这组点数:
o = (106.901024 -6.229162)
stop = (106.804 -6.21861)
我将这个 EPSG 4326 坐标转换为局部投影 EPSG 4808,得到了这个:
o_proj = (0.09183386384156803 -6.229330112968891)
stop_proj = (-0.005201753272169649 -6.218776788266844)
这给出了 0.09760780527657992 的距离度量。 Google 地图为我提供了坐标 o
和 stop
的距离测量值,为 10.79 公里。看起来我的代码的距离度量给出的答案比实际距离小 10^-3 倍。为什么?我的代码正确吗?
您使用的坐标是以度为单位,而不是以米为单位。
你用的Points
class可能不关心这个,计算它们之间的笛卡尔距离
请使用Haversine equation, e.g. the one from here:
from math import radians, cos, sin, asin, sqrt
def haversine(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
"""
# convert decimal degrees to radians
lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
c = 2 * asin(sqrt(a))
r = 6371 # Radius of earth in kilometers. Use 3956 for miles
return c * r
那么你的代码returns正确结果:
from pyproj import Proj, transform
inProj = Proj(init='epsg:4326')
outProj = Proj(init='epsg:4808')
o = (106.901024, -6.229162)
stop = (106.804, -6.21861)
o_proj = transform(inProj, outProj, o[0], o[1])
stop_proj = transform(inProj, outProj, stop[0], stop[1])
print(haversine(o_proj[0], o_proj[1], stop_proj[0], stop_proj[1]))