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

10

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?

robintw
quelle

Antworten:

6

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)
Allen
quelle
Danke - das scheint die realistischste Lösung zu sein. Wenn ich mir die Dokumente ansehe, scheint es, dass ich dies nur zwischen einem Referenzpunkt und einer Reihe anderer Punkte tun kann, also müsste ich es in einer Schleife tun, um alle meine Punkte durchzugehen. Kennen Sie einen effizienteren Weg, dies in R zu tun?
Robintw
Es ist keine Schleife erforderlich. Sie können der Funktion zwei Punktmengen zuweisen, und es wird eine Matrix mit Abständen zwischen den einzelnen Punktkombinationen zurückgegeben. Die bearbeitete Antwort enthält Beispielcode.
Allen
2

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
Michalis Avraam
quelle
1

Ä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äß.

Basto
quelle
Danke für den Vorschlag. Ich kann jedoch nicht sehen, wie mir das helfen wird. Gemäß den Dokumenten ( help.arcgis.com/de/arcgisdesktop/10.0/help/index.html#//… ) "liegt der Abstand in der linearen Einheit des Koordinatensystems der Eingabemerkmale." in lat / lon gibt mir sicher ergebnisse in dezimalgrad? (Ich habe hier keine Maschine mit ArcGIS zum Testen)
robintw
In diesem Fall würde ich wahrscheinlich eine "schnelle und schmutzige" Lösung verwenden, indem ich X- und Y-Felder in Ihre Datentabelle einfüge und auf Geometrie berechnen klicke und X und Y in Meter auswähle. Wenn diese Option nicht ausgewählt werden kann, ändern Sie das Koordinatensystem Ihres MXD. Ich habe zuvor an einem Projekt gearbeitet, bei dem mein Kunde Long / Lat-, X / Y- und Gauss-Krueger-R / H-Werte in jeder Shape-Datei haben wollte. Um komplizierte Berechnungen zu vermeiden, war es am einfachsten, Projektionen zu ändern und die Geometrie zu berechnen.
Basto
0

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)

Mike T.
quelle