Рассчитать расстояние в км до ближайших точек (указано в широте / длине) с помощью ArcGIS DEsktop и / или R?

10

У меня есть два набора точечных данных в ArcGIS, оба из которых приведены в координатах WGS84 широта / долгота, и точки распределены по всему миру. Я хотел бы найти ближайшую точку в наборе данных A к каждой точке в наборе данных B и получить расстояние между ними в километрах.

Это похоже на идеальное использование инструмента Ближний, но это дает мне результаты в системе координат входных точек: то есть в десятичных градусах. Я знаю, что могу перепроектировать данные, но я понимаю ( из этого вопроса ), что трудно (если не невозможно) найти проекцию, которая даст точные расстояния по всему миру.

Ответы на этот вопрос предлагают использовать формулу Хаверсайна для расчета расстояний, используя координаты широта-долгота напрямую. Есть ли способ сделать это и получить результат в км, используя ArcGIS? Если нет, то как лучше всего подойти к этому?

robintw
источник

Ответы:

6

Хотя это не решение ArcGIS, ваша проблема может быть решена в R путем экспорта ваших точек из Arc и использования spDists функции из spпакета. Функция находит расстояния между контрольной точкой (точками) и матрицей точек, в километрах, если вы установили longlat=T.

Вот быстрый и грязный пример:

library(sp)
## Sim up two sets of 100 points, we'll call them set a and set b:
a <- SpatialPoints(coords = data.frame(x = rnorm(100, -87.5), y = rnorm(100, 30)), proj4string=CRS("+proj=longlat +datum=WGS84"))
b <- SpatialPoints(coords = data.frame(x = rnorm(100, -88.5), y = rnorm(100, 30.5)), proj4string=CRS("+proj=longlat +datum=WGS84"))

## Find the distance from each point in a to each point in b, store
##    the results in a matrix.
results <- spDists(a, b, longlat=T)
шестигранный
источник
Спасибо - это кажется самым реалистичным решением. Глядя на документы, кажется, что я могу сделать это только между контрольной точкой и набором других точек, поэтому мне придется делать это в цикле, чтобы просмотреть все мои точки. Вы знаете более эффективный способ сделать это в R?
robintw
Зацикливание не требуется, вы можете дать функции два набора точек, и она вернет матрицу с расстояниями между каждой комбинацией точек. Отредактированный ответ для включения примера кода.
Аллен
3

Это не решение ArcGIS, но использование модели данных Round Earth в пространственной базе данных поможет. Вычислить расстояние от земли в базе данных, поддерживающей это, было бы довольно легко. Я могу предложить вам два чтения:

http://postgis.net/workshops/postgis-intro/geography.html

http://blog.safe.com/2012/08/round-earth-data-in-oracle-postgis-and-sql-server/

Джеб
источник
2

Вам нужен расчет расстояния, который работает с широтой / долготой. Vincenty это один я хотел бы использовать ( с точностью до 0,5 мм). Я играл с ним раньше, и это не так уж сложно использовать.

Код немного длинный, но работает. Учитывая две точки в WGS, он вернет расстояние в метрах.

Вы можете использовать его как скрипт Python в ArcGIS или обернуть его вокруг другого скрипта, который просто перебирает два шейп-файла точек и строит для вас матрицу расстояний. Или, возможно, проще представить результаты GENERATE_NEAR_TABLE, найдя 2-3 ближайших объекта (чтобы избежать осложнений искривления Земли).

import math

ellipsoids = {
    #name        major(m)   minor(m)            flattening factor
    'WGS-84':   (6378137,   6356752.3142451793, 298.25722356300003),
    'GRS-80':   (6378137,   6356752.3141403561, 298.25722210100002),
    'GRS-67':   (6378160,   6356774.5160907144, 298.24716742700002),

}

def distanceVincenty(lat1, long1, lat2, long2, ellipsoid='WGS-84'):
    """Computes the Vicenty distance (in meters) between two points
    on the earth. Coordinates need to be in decimal degrees.
    """
    # Check if we got numbers
    # Removed to save space
    # Check if we know about the ellipsoid
    # Removed to save space
    major, minor, ffactor = ellipsoids[ellipsoid]
    # Convert degrees to radians
    x1 = math.radians(lat1)
    y1 = math.radians(long1)
    x2 = math.radians(lat2)
    y2 = math.radians(long2)
    # Define our flattening f
    f = 1 / ffactor
    # Find delta X
    deltaX = y2 - y1
    # Calculate U1 and U2
    U1 = math.atan((1 - f) * math.tan(x1))
    U2 = math.atan((1 - f) * math.tan(x2))
    # Calculate the sin and cos of U1 and U2
    sinU1 = math.sin(U1)
    cosU1 = math.cos(U1)
    sinU2 = math.sin(U2)
    cosU2 = math.cos(U2)
    # Set initial value of L
    L = deltaX
    # Set Lambda equal to L
    lmbda = L
    # Iteration limit - when to stop if no convergence
    iterLimit = 100
    while abs(lmbda) > 10e-12 and iterLimit >= 0:
        # Calculate sine and cosine of lmbda
        sin_lmbda = math.sin(lmbda)
        cos_lmbda = math.cos(lmbda)
        # Calculate the sine of sigma
        sin_sigma = math.sqrt(
                (cosU2 * sin_lmbda) ** 2 + 
                (cosU1 * sinU2 - 
                 sinU1 * cosU2 * cos_lmbda) ** 2
        )
        if sin_sigma == 0.0:
            # Concident points - distance is 0
            return 0.0
        # Calculate the cosine of sigma
        cos_sigma = (
                    sinU1 * sinU2 + 
                    cosU1 * cosU2 * cos_lmbda
        )
        # Calculate sigma
        sigma = math.atan2(sin_sigma, cos_sigma)
        # Calculate the sine of alpha
        sin_alpha = (cosU1 * cosU2 * math.sin(lmbda)) / (sin_sigma)
        # Calculate the square cosine of alpha
        cos_alpha_sq = 1 - sin_alpha ** 2
        # Calculate the cosine of 2 sigma
        cos_2sigma = cos_sigma - ((2 * sinU1 * sinU2) / cos_alpha_sq)
        # Identify C
        C = (f / 16.0) * cos_alpha_sq * (4.0 + f * (4.0 - 3 * cos_alpha_sq))
        # Recalculate lmbda now
        lmbda = L + ((1.0 - C) * f * sin_alpha * (sigma + C * sin_sigma * (cos_2sigma + C * cos_sigma * (-1.0 + 2 * cos_2sigma ** 2)))) 
        # If lambda is greater than pi, there is no solution
        if (abs(lmbda) > math.pi):
            raise ValueError("No solution can be found.")
        iterLimit -= 1
    if iterLimit == 0 and lmbda > 10e-12:
        raise ValueError("Solution could not converge.")
    # Since we converged, now we can calculate distance
    # Calculate u squared
    u_sq = cos_alpha_sq * ((major ** 2 - minor ** 2) / (minor ** 2))
    # Calculate A
    A = 1 + (u_sq / 16384.0) * (4096.0 + u_sq * (-768.0 + u_sq * (320.0 - 175.0 * u_sq)))
    # Calculate B
    B = (u_sq / 1024.0) * (256.0 + u_sq * (-128.0 + u_sq * (74.0 - 47.0 * u_sq)))
    # Calculate delta sigma
    deltaSigma = B * sin_sigma * (cos_2sigma + 0.25 * B * (cos_sigma * (-1.0 + 2.0 * cos_2sigma ** 2) - 1.0/6.0 * B * cos_2sigma * (-3.0 + 4.0 * sin_sigma ** 2) * (-3.0 + 4.0 * cos_2sigma ** 2)))
    # Calculate s, the distance
    s = minor * A * (sigma - deltaSigma)
    # Return the distance
    return s
Михалис Авраам
источник
1

Я сделал похожий опыт с маленькими наборами данных, используя инструмент Point Distance. При этом вы не можете автоматически найти ближайшие точки в наборе данных A, но, по крайней мере, получите вывод таблицы с полезными результатами в км или m. На следующем шаге вы можете выбрать кратчайшее расстояние до каждой точки набора данных B вне таблицы.

Но этот подход будет зависеть от количества точек в ваших наборах данных. Это может не работать должным образом с большими наборами данных.

Басто
источник
Спасибо за предложение. Однако я не вижу, как это мне поможет. Согласно документации ( help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//… ) «расстояние находится в линейной единице системы координат входных объектов», которая, как и мои входные объекты, в широте / долготе наверняка даст мне результаты в десятичных градусах? (У меня нет машины с ArcGIS для тестирования)
robintw
В этом случае я бы, вероятно, использовал бы «быстрое и грязное» решение, добавив поля X и Y в таблицу данных и щелкнув «Рассчитать геометрию», выбрав «X» и «Y» в метрах. Если невозможно выбрать эту опцию, измените систему координат вашего MXD. Раньше я работал над проектом, где мой клиент хотел получить значения long / lat, X / Y и Gauss-Krueger R / H в каждом файле Shape. Чтобы избежать сложных вычислений, проще всего было просто изменить проекции и рассчитать геометрию.
Басто
0

Если вам нужны высокоточные и надежные геодезические измерения, используйте GeographicLib , которая изначально написана на нескольких языках программирования, включая C ++, Java, MATLAB, Python и т. Д.

См. CFF Karney (2013) «Алгоритмы для геодезических» для литературной ссылки. Обратите внимание, что эти алгоритмы являются более надежными и точными, чем алгоритм Винсенти, например, рядом с антиподами.

Чтобы вычислить расстояние в метрах между двумя точками, получите s12атрибут расстояния из обратного геодезического решения . Например, с пакетом geographiclib для Python

from geographiclib.geodesic import Geodesic
g = Geodesic.WGS84.Inverse(-41.32, 174.81, 40.96, -5.50)
print(g)  # shows:
{'a12': 179.6197069334283,
 'azi1': 161.06766998615873,
 'azi2': 18.825195123248484,
 'lat1': -41.32,
 'lat2': 40.96,
 'lon1': 174.81,
 'lon2': -5.5,
 's12': 19959679.26735382}

Или сделайте удобную функцию, которая также конвертирует из метров в километры:

dist_km = lambda a, b: Geodesic.WGS84.Inverse(a[0], a[1], b[0], b[1])['s12'] / 1000.0
a = (-41.32, 174.81)
b = (40.96, -5.50)
print(dist_km(a, b))  # 19959.6792674 km

Теперь, чтобы найти ближайшую точку между списками Aи B, каждый из которых 100 баллов:

from random import uniform
from itertools import product
A = [(uniform(-90, 90), uniform(-180, 180)) for x in range(100)]
B = [(uniform(-90, 90), uniform(-180, 180)) for x in range(100)]
a_min = b_min = min_dist = None
for a, b in product(A, B):
    d = dist_km(a, b)
    if min_dist is None or d < min_dist:
        min_dist = d
        a_min = a
        b_min = b

print('%.3f km between %s and %s' % (min_dist, a_min, b_min))

22,481 км между (84,57916462672875, 158,67545706102192) и (84,70326937581333, 156,9784597422855)

Майк Т
источник