Hitung jarak dalam km ke titik terdekat (diberikan dalam lat / panjang) menggunakan ArcGIS DEsktop dan / atau R?

10

Saya memiliki dua set data titik di ArcGIS, keduanya diberikan dalam koordinat WGS84 lat / lon dan poin-poinnya tersebar di seluruh dunia. Saya ingin mencari titik terdekat di Dataset A ke setiap titik di Dataset B, dan mendapatkan jarak di antara mereka dalam kilometer.

Ini tampaknya seperti penggunaan sempurna alat Dekat, tapi itu memberi saya hasil dalam sistem koordinat titik input: yaitu, derajat desimal. Saya tahu saya bisa memproyeksikan ulang data, tetapi saya mengumpulkan ( dari pertanyaan ini ) bahwa sulit (jika bukan tidak mungkin) untuk menemukan proyeksi yang akan memberikan jarak akurat di seluruh dunia.

Jawaban atas pertanyaan itu menyarankan menggunakan rumus Haversine untuk menghitung jarak menggunakan koordinat garis lintang-bujur secara langsung. Apakah ada cara untuk melakukan ini dan mendapatkan hasil dalam km menggunakan ArcGIS? Jika tidak, apa cara terbaik untuk mendekati ini?

Robintw
sumber

Jawaban:

6

Meskipun ini bukan solusi ArcGIS, masalah Anda dapat diselesaikan dalam R dengan mengekspor poin Anda dari Arc dan menggunakan spDists fungsi dari sppaket. Fungsi menemukan jarak antara titik referensi dan matriks titik, dalam kilometer jika Anda atur longlat=T.

Ini contoh cepat dan kotor:

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
sumber
Terima kasih - ini sepertinya solusi paling realistis. Melihat dokumen sepertinya saya hanya bisa melakukan ini antara titik referensi dan satu set poin lain, jadi saya harus melakukannya dalam satu lingkaran untuk melewati semua poin saya. Apakah Anda tahu cara yang lebih efisien untuk melakukan ini di R?
robintw
Tidak diperlukan perulangan, Anda bisa memberi fungsi dua set titik dan itu akan mengembalikan matriks dengan jarak antara setiap kombinasi titik. Jawaban yang diedit untuk memasukkan kode contoh.
allen
2

Anda membutuhkan perhitungan jarak yang berfungsi dengan Lat / Long. Vincenty adalah yang akan saya gunakan (akurasi 0,5mm). Saya telah bermain dengannya sebelumnya, dan tidak terlalu sulit untuk digunakan.

Kode ini agak panjang, tetapi berhasil. Diberi dua poin dalam WGS, ia akan mengembalikan jarak dalam meter.

Anda bisa menggunakan ini sebagai skrip Python di ArcGIS, atau membungkusnya dengan skrip lain yang hanya diulangi pada dua Point Shapefile dan membuat matriks jarak untuk Anda. Atau, mungkin lebih mudah untuk memberi makan hasil GENERATE_NEAR_TABLE dengan menemukan 2-3 fitur terdekat (untuk menghindari komplikasi kelengkungan bumi).

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
sumber
1

Saya membuat pengalaman serupa dengan dataset kecil menggunakan alat Point Distance. Dengan melakukan itu, Anda tidak dapat secara otomatis menemukan titik terdekat di Dataset A Anda, tetapi setidaknya dapatkan output tabel dengan hasil km atau m yang bermanfaat. Pada langkah berikutnya Anda bisa memilih jarak terpendek ke setiap titik Dataset B di luar tabel.

Tetapi pendekatan ini akan tergantung pada jumlah poin dalam dataset Anda. Mungkin tidak berfungsi dengan baik dengan kumpulan data besar.

basto
sumber
Terima kasih untuk sarannya. Namun, saya tidak bisa melihat bagaimana itu akan membantu saya. Menurut dokumen ( help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//… ) "jaraknya ada dalam unit linear dari sistem koordinat fitur input.", Yang sebagai fitur input saya adalah di lat / lon pasti akan memberi saya hasil dalam derajat desimal? (Saya belum punya mesin dengan ArcGIS di sini untuk menguji)
robintw
Dalam hal ini, saya mungkin akan menggunakan solusi "cepat dan kotor" dengan menambahkan bidang X dan Y di dataat Anda dan klik Calculate Geometry memilih X dan Y dalam meter. Jika tidak memungkinkan untuk memilih opsi ini, ubah sistem koordinat MXD Anda. Saya sedang mengerjakan sebuah proyek sebelumnya, di mana klien saya menginginkan nilai panjang / lat, X / Y dan Gauss-Krueger R / H semua dalam setiap file Shape. Untuk menghindari perhitungan yang rumit, mengubah proyeksi dan menghitung geometri adalah cara yang paling mudah.
basto
0

Jika Anda membutuhkan pengukuran geodesik presisi tinggi dan kuat, gunakan GeographicLib , yang aslinya ditulis dalam beberapa bahasa pemrograman, termasuk C ++, Java, MATLAB, Python, dll.

Lihat CFF Karney (2013) "Algoritma untuk geodesik" untuk referensi sastra. Perhatikan bahwa algoritma ini lebih kuat dan akurat daripada algoritma Vincenty, misalnya di dekat antipode.

Untuk menghitung jarak dalam meter antara dua titik, dapatkan s12atribut jarak dari solusi geodesik terbalik . Misalnya, dengan paket geographiclib untuk 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}

Atau buat fungsi kenyamanan, yang juga mengubah dari meter ke kilometer:

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

Sekarang untuk menemukan titik terdekat antara daftar Adan B, masing-masing dengan 100 poin:

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 antara (84.57916462672875, 158.67545706102192) dan (84.70326937581333, 156.9784597422855)

Mike T
sumber