Calcular a distância em km até os pontos mais próximos (dados em lat / long) usando ArcGIS DEsktop e / ou R?

10

Eu tenho dois conjuntos de dados de pontos no ArcGIS, ambos dados nas coordenadas latinas / longínquas WGS84 e os pontos estão espalhados por todo o mundo. Gostaria de encontrar o ponto mais próximo no conjunto de dados A para cada ponto no conjunto de dados B e obter a distância entre eles em quilômetros.

Parece um uso perfeito da ferramenta Near, mas isso me dá resultados no sistema de coordenadas dos pontos de entrada: ou seja, graus decimais. Sei que poderia re-projetar os dados, mas deduzo (a partir dessa pergunta ) que é difícil (se não impossível) encontrar uma projeção que ofereça distâncias precisas em todo o mundo.

As respostas a essa pergunta sugerem o uso da fórmula de Haversine para calcular distâncias usando as coordenadas de latitude-longitude diretamente. Existe uma maneira de fazer isso e obter um resultado em km usando o ArcGIS? Caso contrário, qual é a melhor maneira de abordar isso?

robintw
fonte

Respostas:

6

Embora essa não seja uma solução do ArcGIS, seu problema pode ser resolvido no R exportando seus pontos do Arc e usando a spDists função do sppacote. A função localiza as distâncias entre um (s) ponto (s) de referência e uma matriz de pontos, em quilômetros, se você definir longlat=T.

Aqui está um exemplo rápido e sujo:

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
fonte
Obrigado - esta parece ser a solução mais realista. Olhando para os documentos, parece que eu só posso fazer isso entre um ponto de referência e um conjunto de outros pontos, então eu teria que fazer isso em um loop para passar por todos os meus pontos. Você conhece uma maneira mais eficiente de fazer isso no R?
robintw
Não é necessário loop, você pode atribuir à função dois conjuntos de pontos e ele retornará uma matriz com distâncias entre cada combinação de pontos. Resposta editada para incluir código de exemplo.
Allen
2

Você precisa de um cálculo de distância que funcione com Lat / Long. Vincenty é o que eu usaria (precisão de 0,5 mm). Já brinquei com ele antes e não é muito difícil de usar.

O código é um pouco longo, mas funciona. Dados dois pontos no WGS, ele retornará uma distância em metros.

Você pode usar isso como um script Python no ArcGIS, ou envolvê-lo em outro script que simplesmente itere sobre os dois Shapefiles de ponto e constrói uma matriz de distância para você. Ou, provavelmente, é mais fácil alimentar os resultados de GENERATE_NEAR_TABLE com a localização das 2-3 características mais próximas (para evitar complicações da curvatura da terra).

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

Fiz experiências semelhantes com pequenos conjuntos de dados usando a ferramenta Point Distance. Fazendo isso, você não pode encontrar automaticamente os pontos mais próximos no seu conjunto de dados A, mas pelo menos obter uma saída de tabela com resultados úteis de km ou m. Em uma próxima etapa, você pode selecionar a menor distância para cada ponto do conjunto de dados B da tabela.

Mas essa abordagem dependeria da quantidade de pontos em seus conjuntos de dados. Pode não funcionar corretamente com grandes conjuntos de dados.

basto
fonte
Obrigado pela sugestão. No entanto, não consigo ver como isso vai me ajudar. De acordo com os documentos ( help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//… ) "a distância está na unidade linear do sistema de coordenadas dos recursos de entrada.", Que como meus recursos de entrada são em lat / lon certamente me dará resultados em graus decimais? (Eu não tenho uma máquina com ArcGIS sobre ele aqui para teste)
robintw
Nesse caso, eu provavelmente usaria uma solução "rápida e suja" adicionando campos X e Y em sua tabela de dados e clique em Calcular Geometria escolhendo X e Y em metros. Se não for possível escolher essa opção, altere o sistema de coordenadas do seu MXD. Eu estava trabalhando em um projeto antes, onde meu cliente desejava valores long / lat, X / Y e Gauss-Krueger R / H em todos os arquivos Shape. Para evitar cálculos complicados, alterar as projeções e calcular a geometria foi a maneira mais fácil de resolver.
basto
0

Se você precisar de medições geodésicas robustas e de alta precisão, use GeographicLib , que é originalmente escrito em várias linguagens de programação, incluindo C ++, Java, MATLAB, Python, etc.

Veja CFF Karney (2013) "Algoritmos para geodésica" para uma referência literária. Observe que esses algoritmos são mais robustos e precisos que o algoritmo de Vincenty, por exemplo, perto de antípodas.

Para calcular a distância em metros entre dois pontos, obtenha o s12atributo de distância da solução geodésica inversa . Por exemplo, com o pacote geographiclib para 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 crie uma função de conveniência, que também converte metros para quilômetros:

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

Agora, encontre o ponto mais próximo entre as listas Ae B, cada um com 100 pontos:

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

Mike T
fonte