Calculer la distance en km jusqu'aux points les plus proches (donnée en lat / long) en utilisant ArcGIS DEsktop et / ou R?

J'ai deux jeux de données ponctuelles dans ArcGIS, qui sont tous deux donnés en coordonnées lat / lon WGS84 et les points sont répartis dans le monde entier. Je voudrais trouver le point le plus proche dans l'ensemble de données A de chaque point dans l'ensemble de données B, et obtenir la distance entre eux en kilomètres.

Cela semble être une utilisation parfaite de l'outil Near, mais cela me donne des résultats dans le système de coordonnées des points d'entrée: c'est-à-dire des degrés décimaux. Je sais que je pourrais re-projeter les données, mais je suppose (à partir de cette question ) qu'il est difficile (voire impossible) de trouver une projection qui donnera des distances précises partout dans le monde.

Les réponses à cette question suggèrent d'utiliser la formule Haversine pour calculer les distances en utilisant directement les coordonnées latitude-longitude. Existe-t-il un moyen de le faire et d'obtenir un résultat en km en utilisant ArcGIS? Sinon, quelle est la meilleure façon d'aborder cela?

Bien que ce ne soit pas une solution ArcGIS, votre problème peut être résolu dans R en exportant vos points depuis Arc et en utilisant la spDists fonction du sppackage. La fonction recherche les distances entre un ou des points de référence et une matrice de points, en kilomètres si vous définissez longlat=T.

Voici un exemple rapide et sale:

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)

Ce n'est pas une solution ArcGIS, mais l'utilisation d'un modèle de données Round Earth dans une base de données spatiale ferait l'affaire. Calculer la distance de la terre dans une base de données prenant en charge ce serait assez facile. Je peux vous proposer deux lectures:

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

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

Vous avez besoin d'un calcul de distance qui fonctionne avec Lat / Long. Vincenty est celui que j'utiliserais (précision 0,5 mm). J'ai déjà joué avec, et ce n'est pas trop difficile à utiliser.

Le code est un peu long, mais ça marche. Compte tenu de deux points en WGS, il retournera une distance en mètres.

Vous pouvez l'utiliser comme un script Python dans ArcGIS, ou l'enrouler autour d'un autre script qui itère simplement sur les deux fichiers de formes de points et crée une matrice de distance pour vous. Ou, il est probablement plus facile d'alimenter les résultats de GENERATE_NEAR_TABLE en trouvant les 2-3 entités les plus proches (pour éviter les complications de la courbure de la Terre).

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

J'ai fait des expériences similaires avec de petits ensembles de données en utilisant l'outil Point Distance. Ce faisant, vous ne pouvez pas trouver automatiquement les points les plus proches dans votre ensemble de données A, mais obtenez au moins une sortie de tableau avec des résultats utiles en km ou en m. Dans une étape suivante, vous pouvez sélectionner la distance la plus courte jusqu'à chaque point du jeu de données B hors de la table.

Mais cette approche dépendrait de la quantité de points dans vos jeux de données. Il peut ne pas fonctionner correctement avec de grands ensembles de données.

Si vous avez besoin de mesures géodésiques de haute précision et robustes, utilisez GeographicLib , qui est écrit en natif dans plusieurs langages de programmation, y compris C ++, Java, MATLAB, Python, etc.

Voir CFF Karney (2013) «Algorithms for geodesics» pour une référence littéraire. Notez que ces algorithmes sont plus robustes et précis que l'algorithme de Vincenty, par exemple près des antipodes.

Pour calculer la distance en mètres entre deux points, obtenez l' s12attribut de distance à partir de la solution géodésique inverse . Par exemple, avec le paquet GeographicLib pour 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}

Ou créez une fonction de confort, qui convertit également des mètres en kilomètres:

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

Maintenant, pour trouver le point le plus proche entre les listes Aet B, chacun avec 100 points:

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) et (84.70326937581333, 156.9784597422855)