Почему мой код вычисляет неверное расстояние между преобразованными координатами? - PullRequest
2 голосов
/ 19 июня 2019

У меня большой массив точек, и я пытаюсь найти расстояние между этими точками и другим набором точек.Изначально я использовал функцию геопандий to_crs для преобразования crs, чтобы я мог получить точное измерение расстояния в метрах, когда я делаю df.distance(point).Тем не менее, поскольку файл очень большой, преобразование crs файла заняло слишком много времени.Код работал в течение 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 раза меньше, чем фактическое расстояние.Почему это так?Мой код правильный?

1 Ответ

1 голос
/ 19 июня 2019

Используемые вами координаты даны в градусах, а не в метрах.

Используемый вами класс Points, вероятно, не заботится об этом и вычисляет декартово расстояние между ними.

Пожалуйста, используйте уравнение Haversine , например один из здесь :

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

Тогда ваш код возвращает правильный результат:

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]))
...