Berechnen Sie die Entfernung in km zu den nächstgelegenen Punkten (angegeben in lat / long) mit ArcGIS DEsktop und / oder R?

Ich habe zwei Punktdatensätze in ArcGIS, die beide in WGS84-Lat / Lon-Koordinaten angegeben sind und die Punkte auf der ganzen Welt verteilt sind. Ich möchte den nächstgelegenen Punkt in Datensatz A zu jedem Punkt in Datensatz B finden und die Entfernung zwischen ihnen in Kilometern ermitteln.

Dies scheint eine perfekte Verwendung des Near-Werkzeugs zu sein, aber das gibt mir Ergebnisse im Koordinatensystem der Eingabepunkte: dh Dezimalgrade. Ich weiß, dass ich die Daten neu projizieren könnte, aber ich stelle ( aus dieser Frage ) fest, dass es schwierig (wenn nicht unmöglich) ist, eine Projektion zu finden, die genaue Entfernungen auf der ganzen Welt liefert.

Die Antworten auf diese Frage legen nahe, die Haversine-Formel zu verwenden, um Entfernungen direkt anhand der Längen- und Breitengradkoordinaten zu berechnen. Gibt es eine Möglichkeit, mit ArcGIS ein Ergebnis in km zu erzielen? Wenn nicht, wie geht man das am besten an?

Obwohl dies keine ArcGIS-Lösung ist, kann Ihr Problem in R gelöst werden, indem Sie Ihre Punkte aus Arc exportieren und die spDists Funktion aus dem spPaket verwenden. Die Funktion ermittelt die Abstände zwischen einem oder mehreren Referenzpunkten und einer Punktmatrix in Kilometern, wenn Sie diese festlegen longlat=T.

Hier ist ein schnelles und schmutziges Beispiel:

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)

Es ist keine ArcGIS-Lösung, aber die Verwendung eines Round Earth-Datenmodells in einer räumlichen Datenbank würde den Trick tun. Die Berechnung der Erdentfernung in einer Datenbank, die dies unterstützt, wäre ziemlich einfach. Ich kann Ihnen zwei Lesungen vorschlagen:

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

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

Sie benötigen eine Entfernungsberechnung, die mit Lat / Long funktioniert. Vincenty ist derjenige, den ich verwenden würde (0,5 mm Genauigkeit). Ich habe schon einmal damit gespielt und es ist nicht zu schwer zu benutzen.

Der Code ist etwas lang, funktioniert aber. Bei zwei Punkten in WGS wird eine Entfernung in Metern zurückgegeben.

Sie können dies als Python-Skript in ArcGIS verwenden oder es um ein anderes Skript wickeln, das einfach die beiden Punktformdateien durchläuft und eine Distanzmatrix für Sie erstellt. Oder es ist wahrscheinlich einfacher, die Ergebnisse von GENERATE_NEAR_TABLE mit der Suche nach den 2-3 nächstgelegenen Merkmalen zu versorgen (um Komplikationen der Erdkrümmung zu vermeiden).

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

Ähnliche Erfahrungen habe ich mit kleinen Datensätzen mit dem Punktabstandstool gemacht. Auf diese Weise können Sie nicht automatisch die nächstgelegenen Punkte in Ihrem Datensatz A finden, sondern erhalten zumindest eine Tabellenausgabe mit nützlichen km- oder m-Ergebnissen. In einem nächsten Schritt können Sie den kürzesten Abstand zu jedem Punkt von Datensatz B aus der Tabelle auswählen.

Dieser Ansatz hängt jedoch von der Anzahl der Punkte in Ihren Datensätzen ab. Bei großen Datenmengen funktioniert dies möglicherweise nicht ordnungsgemäß.

Wenn Sie hochpräzise und robuste geodätische Messungen benötigen, verwenden Sie GeographicLib , das nativ in mehreren Programmiersprachen geschrieben ist, darunter C ++, Java, MATLAB, Python usw.

Siehe CFF Karney (2013) "Algorithmen für Geodäten" für eine literarische Referenz. Beachten Sie, dass diese Algorithmen robuster und genauer sind als der Algorithmus von Vincenty, beispielsweise in der Nähe von Antipoden.

Um die Entfernung in Metern zwischen zwei Punkten zu berechnen, rufen Sie das s12Entfernungsattribut aus der inversen geodätischen Lösung ab . Zum Beispiel mit dem geographiclib- Paket für 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}

Oder machen Sie eine Komfortfunktion, die auch von Metern in Kilometer umrechnet:

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

Finden Sie nun den nächstgelegenen Punkt zwischen Listen Aund Bmit jeweils 100 Punkten:

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 zwischen (84.57916462672875, 158.67545706102192) und (84.70326937581333, 156.9784597422855)