¿Calcular la distancia en km a puntos más cercanos (en lat/long) utilizando ArcGIS DEsktop o R?

Question

Tengo dos puntos de los conjuntos de datos en ArcGIS, ambos de los cuales son dados en WGS84 (lat/lon coordenadas y los puntos están distribuidos en todo el mundo. Me gustaría encontrar el punto más cercano en Un conjunto de datos para cada punto del conjunto de datos B, y obtener la distancia entre ellos en kilómetros.

Esto parece como un perfecto uso de la herramienta de proximidad, pero que me da resultados en el sistema de coordenadas de los puntos de entrada: es decir, grados decimales. Sé que podría re-proyectar los datos, pero tengo entendido (a partir de esta pregunta) que es difícil (si no imposible) encontrar una proyección que dará precisa de las distancias en todo el mundo.

Las respuestas a esta pregunta que sugieren el uso de la fórmula de Haversine para el cálculo de distancias a partir de la latitud y longitud en coordenadas directamente. Hay una manera de hacer esto y conseguir un resultado en km utilizando ArcGIS? Si no, ¿cuál es la mejor manera de abordar esto?

Answer 1

Aunque no es una solución de ArcGIS, tu problema puede ser resuelto en R por exportar sus puntos de arco y usando el spDists de la función de la sp paquete. La función encuentra en kilómetros las distancias entre un punto de referencia y matriz de puntos, si se establece longlat=T .

Aquí está un ejemplo rápido y sucio:

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)

Answer 2

No es una solución ArcGIS, pero el uso de una Tierra Redonda modelo de datos en una base de datos espaciales haría el truco. Cálculo de la distancia de la tierra en la base de datos que soporte este iba a ser bastante fácil. Puedo sugerir dos lecturas:

http://workshops.opengeo.org/postgis-intro/geography.html

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

Answer 3

Usted necesita un cálculo de la distancia que funciona con Lat/Long. Vincenty es el que yo uso (0,5 mm de precisión). He jugado con él antes, y no es demasiado difícil de usar.

El código es un poco largo, pero funciona. Dados dos puntos en WGS, devolverá una distancia en metros.

Usted puede utilizar esto como una secuencia de comandos de Python en ArcGIS, o se envuelve alrededor de otro script que simplemente se itera sobre el Punto dos de los Shapefiles y construye una matriz de distancias para usted. O, es probablemente más fácil para alimentar a los resultados de GENERATE_NEAR_TABLE con la búsqueda de los 2-3 más cercano características (para evitar las complicaciones de la tierra de la curvatura).

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

Answer 4

He hecho experiencias similares con pequeños conjuntos de datos mediante la herramienta de distancia de los puntos. Al hacerlo, usted no puede encontrar automáticamente los puntos más cercanos en el Dataset A, pero al menos Haz una tabla de salida con útil km o m resultados. En un próximo paso puede seleccionar la distancia más corta a cada punto del conjunto de datos B de la tabla.

Pero este enfoque dependerá de la cantidad de puntos en sus bases de datos. No podría funcionar correctamente con conjuntos de datos grandes.

Answer 5

Si se necesita una alta precisión y robusto geodésica mediciones, uso GeographicLib, que es originalmente escritos en varios lenguajes de programación, incluyendo C++, Java, MATLAB, Python, etc.

Ver C. F. F. Karney (2013) "Algoritmos para geodesics" para una referencia literaria. Tenga en cuenta que estos algoritmos son más robustos y precisos que Vincenty del algoritmo, por ejemplo cerca de las antípodas.

Para calcular la distancia en metros entre dos puntos, obtener el s12 distancia atributo de la inversa de la geodésica solución. E. g., con el geographiclib paquete de 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}

O hacer una función de conveniencia, que también convierte de metros a kilómetros de distancia:

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

Ahora para encontrar el punto más cercano entre las listas de A y B, cada una con 100 puntos:

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 km entre (84.57916462672875, 158.67545706102192) y (84.70326937581333, 156.9784597422855)