Por ejemplo, tengo coordenadas para tres puntos base en una costa y necesito encontrar las coordenadas del punto frente a la costa que es equidistante de los tres. Es un ejercicio simple de geometría, pero todas las mediciones deben tener en cuenta la geodesia.
Si me acercara a esto de una manera euclidiana, podría medir los caminos geodésicos que conectan los puntos base, encontrar los puntos medios de los lados del triángulo resultante y crear ortodromos perpendiculares a cada uno de esos caminos. Los tres loxodromos probablemente convergerían en el punto equidistante. Si este es el método correcto, debe haber una forma más fácil de hacerlo en Arc.
Respuestas:
Esta respuesta se divide en múltiples secciones:
Análisis y reducción del problema , que muestra cómo encontrar el punto deseado con rutinas "enlatadas".
Ilustración: un prototipo de trabajo , dando un código de trabajo.
Ejemplo , que muestra ejemplos de las soluciones.
Dificultades , discutir posibles problemas y cómo hacer frente a ellos.
Implementación de ArcGIS , comentarios sobre la creación de una herramienta ArcGIS personalizada y dónde obtener las rutinas necesarias.
Análisis y reducción del problema
Comencemos por observar que en el modelo esférico (perfectamente redonda) habrá siempre ser una solución --en realidad, exactamente dos soluciones. Dados los puntos base A, B y C, cada par determina su "bisectriz perpendicular", que es el conjunto de puntos equidistantes de los dos puntos dados. Esta bisectriz es una geodésica (gran círculo). La geometría esférica es elíptica : cualquiera de las dos geodésicas se cruzan (en dos puntos únicos). Por lo tanto, los puntos de intersección de la bisectriz de AB y la bisectriz de BC son, por definición, equidistantes de A, B y C, resolviendo así el problema. (Vea la primera figura a continuación).
Las cosas se ven más complicadas en un elipsoide, pero debido a que es una pequeña perturbación de la esfera, podemos esperar un comportamiento similar. (El análisis de esto nos llevaría demasiado lejos). Sin embargo, las complicadas fórmulas utilizadas (internamente dentro de un SIG) para calcular distancias precisas en un elipsoide no son una complicación conceptual: el problema es básicamente el mismo. Para ver cuán simple es realmente el problema, vamos a exponerlo de manera algo abstracta. En esta declaración, "d (U, V)" se refiere a la distancia verdadera y totalmente precisa entre los puntos U y V.
Estas tres distancias dependen de la X desconocida . Por lo tanto, las diferencias en distancias u (X) = d (X, A) - d (X, B) yv (X) = d (X, B) - d (X, C) son funciones de valor real de X. De nuevo, de manera algo abstracta, podemos agrupar estas diferencias en un par ordenado. También usaremos (lat, lon) como coordenadas para X, lo que nos permite considerarlo también como un par ordenado, digamos X = (phi, lambda). En esta configuración, la función
es una función de una porción de un espacio bidimensional que toma valores en un espacio bidimensional y nuestro problema se reduce a
Aquí es donde vale la pena la abstracción: existe un gran software excelente para resolver este problema (puramente numérico, multidimensional, búsqueda de raíces). La forma en que funciona es que escribe una rutina para calcular F , luego la pasa al software junto con cualquier información sobre restricciones en su entrada ( phi debe estar entre -90 y 90 grados y lambda debe estar entre -180 y 180 grados). Se arranca por una fracción de segundo y devuelve (típicamente) solo un valor de ( phi , lambda ), si puede encontrar uno.
Hay detalles que manejar, porque hay un arte en esto: hay varios métodos de solución para elegir, dependiendo de cómo F "se comporta"; ayuda a "dirigir" el software dándole un punto de partida razonable para su búsqueda (esta es una forma en que podemos obtener la solución más cercana , en lugar de cualquier otra); y generalmente necesita especificar qué tan precisa le gustaría que fuera la solución (para que sepa cuándo detener la búsqueda). (Para obtener más información sobre lo que los analistas de SIG necesitan saber sobre estos detalles, que surgen mucho en los problemas de SIG, visite Recomendar temas para que se incluyan en un curso de Ciencias de la Computación para Tecnologías Geoespaciales y busque en la sección "Miscelánea" al final. )
Ilustración: un prototipo de trabajo
El análisis muestra que necesitamos programar dos cosas: una estimación inicial cruda de la solución y el cálculo de F en sí.
La estimación inicial se puede hacer mediante un "promedio esférico" de los tres puntos base. Esto se obtiene al representarlos en coordenadas cartesianas geocéntricas (x, y, z), promediando esas coordenadas, y proyectando ese promedio de regreso a la esfera y reexpresándolo en latitud y longitud. El tamaño de la esfera es irrelevante y, por lo tanto, los cálculos se hacen directos: dado que este es solo un punto de partida, no necesitamos cálculos elipsoidales.
Para este prototipo de trabajo usé Mathematica 8.
(La
Ifcondición final comprueba si el promedio puede fallar claramente para indicar una longitud; si es así, vuelve a una media aritmética recta de las latitudes y longitudes de su entrada; tal vez no sea una gran opción, pero al menos una válida. Para aquellos que usan este código como guía de implementación, tenga en cuenta que los argumentos de MathematicaArcTanse invierten en comparación con la mayoría de las otras implementaciones: su primer argumento es la coordenada x, su segundo es la coordenada y, y devuelve el ángulo formado por el vector ( x, y).)En cuanto a la segunda parte, debido a que Mathematica , como ArcGIS y casi todos los demás SIG, contiene código para calcular distancias precisas en el elipsoide, no hay casi nada que escribir. Simplemente llamamos a la rutina de búsqueda de raíz:
El aspecto más notable de esta implementación es cómo esquiva la necesidad de restringir la latitud (
f) y la longitud (q) calculando siempre el módulo 180 y 360 grados, respectivamente. Esto evita tener que restringir el problema (que a menudo crea complicaciones). Los parámetros de control,MaxIterationsetc., se modifican para que este código proporcione la mayor precisión posible.Para verlo en acción, apliquémoslo a los tres puntos básicos dados en una pregunta relacionada :
Las distancias calculadas entre esta solución y los tres puntos son
(estos son metros). Están de acuerdo con el undécimo dígito significativo (que es demasiado preciso, en realidad, ya que las distancias rara vez son precisas mejor que un milímetro más o menos). Aquí hay una imagen de estos tres puntos (negro), sus tres bisectrices mutuas y la solución (rojo):
Ejemplo
Para probar esta implementación y obtener una mejor comprensión de cómo se comporta el problema, aquí hay una gráfica de contorno de la discrepancia cuadrática media en las distancias para tres puntos base ampliamente espaciados. (La discrepancia RMS se obtiene calculando las tres diferencias d (X, A) -d (X, B), d (X, B) -d (X, C) yd (X, C) -d (X , A), promediando sus cuadrados y sacando la raíz cuadrada. Es igual a cero cuando X resuelve el problema y, de lo contrario, aumenta a medida que X se aleja de una solución, y así mide cuán "cercanos" estamos para ser una solución en cualquier ubicación. )
Los puntos base (60, -120), (10, -40) y (45,10) se muestran en rojo en esta proyección de Plate Carree; la solución (49.2644488, -49.9052992), que requirió 0.03 segundos para calcular, está en amarillo. Su discrepancia RMS es inferior a tres nanómetros , a pesar de que todas las distancias relevantes son miles de kilómetros. Las áreas oscuras muestran valores pequeños del RMS y las áreas claras muestran valores altos.
Este mapa muestra claramente que hay otra solución cerca (-49.2018206, 130.0297177) (calculada a un RMS de dos nanómetros al establecer el valor de búsqueda inicial diametralmente opuesto a la primera solución).
Trampas
Inestabilidad numérica
Cuando los puntos base están casi colineales y muy juntos, todas las soluciones estarán a casi medio mundo de distancia y serán extremadamente difíciles de precisar con precisión. La razón es que pequeños cambios en una ubicación en todo el mundo, moviéndolo hacia o desde los puntos base, inducirán solo cambios increíblemente pequeños en las diferencias de distancias. Simplemente no hay suficiente precisión y precisión incorporada en el cálculo habitual de distancias geodésicas para precisar el resultado.
Por ejemplo, comenzando con los puntos base en (45.001, 0), (45, 0) y (44.999,0), que están separados a lo largo del meridiano principal por solo 111 metros entre cada par,
triobtiene la solución (11.8213, 77.745 ) Las distancias desde él a los puntos base son 8,127,964.998 77; 8.127.964.998 41; y 8,127,964.998 65 metros, respectivamente. ¡Están de acuerdo con el milímetro más cercano! No estoy seguro de cuán preciso puede ser este resultado, pero no me sorprendería en lo más mínimo si otras implementaciones devolvieran ubicaciones lejos de esta, mostrando casi la misma igualdad de las tres distancias.Tiempo de cómputo
Estos cálculos, debido a que implican una búsqueda considerable utilizando cálculos de distancia complicados, no son rápidos, por lo general requieren una fracción notable de un segundo. Las aplicaciones en tiempo real deben ser conscientes de esto.
Implementación de ArcGIS
Python es el entorno de secuencias de comandos preferido para ArcGIS (a partir de la versión 9). El paquete scipy.optimize tiene un buscador de raíz multivariado
rootque debe hacer lo queFindRoothace en el código de Mathematica . Por supuesto, ArcGIS ofrece cálculos precisos de distancia elipsoidal. El resto, entonces, son todos los detalles de implementación: decida cómo se obtendrán las coordenadas del punto base (de una capa, escrita por el usuario, de un archivo de texto, del mouse) y cómo se presentará la salida (como coordenadas visualizado en la pantalla? como un punto gráfico? como un nuevo objeto de punto en una capa?), escriba esa interfaz, conecte el código de Mathematica que se muestra aquí (sencillo), y estará listo.fuente
Como observa, este problema surge al determinar los límites marítimos; a menudo se le conoce como el problema de "tres puntos" y puede buscarlo en Google y encontrar varios documentos que lo aborden. Uno de estos documentos es mío (!) Y ofrezco una solución precisa y rápidamente convergente. Consulte la Sección 14 de http://arxiv.org/abs/1102.1215
El método consta de los siguientes pasos:
La fórmula necesaria para la solución de tres puntos en la proyección se da en el documento. Mientras utilice una proyección equidistante azimutal precisa, la respuesta será exacta. La convergencia es cuadrática, lo que significa que solo se necesitan unas pocas iteraciones. Esto seguramente superará a los métodos generales de búsqueda de raíces sugeridos por @whuber.
No puedo ayudarte directamente con ArcGIS. Puede obtener mi paquete de Python para hacer cálculos geodésicos desde https://pypi.python.org/pypi/geographiclib y codificar la proyección en función de esto es simple.
Editar
Cayley consideró el problema de encontrar el punto triple en el caso degenerado de @ whuber (45 + eps, 0) (45,0) (45-eps, 0) en En las líneas geodésicas de un esferoide achatado , Phil. revista (1870), http://books.google.com/books?id=4XGIOoCMYYAC&pg=PA15
En este caso, el punto triple se obtiene siguiendo una geodésica de (45,0) con azimut 90 y encontrando el punto en el que se desvanece la escala geodésica. Para el elipsoide WGS84, este punto es (-0.10690908732248, 89.89291072793167). La distancia desde este punto a cada uno de (45.001,0), (45,0), (44.999) es 10010287.665788943 m (dentro de un nanómetro más o menos). Esto es aproximadamente 1882 km más que la estimación de Whuber (lo que solo muestra cuán inestable es este caso). Para una tierra esférica, el punto triple sería (0,90) o (0, -90), por supuesto.
ADENDA: Aquí hay una implementación del método equidistante azimutal usando Matlab
Probar esto usando Octave me sale
como el punto triple para Nueva York, Tokio y Wellington.
Este método es inexacto para los puntos colineales vecinos, por ejemplo, [45.001,0], [45,0], [44.999,0]. En ese caso, debe resolver para M 12 = 0 en una geodésica que emana de [45,0] en el acimut 90. La siguiente función hace el truco (usando el método de Newton):
Por ejemplo, esto da:
fuente
Tenía curiosidad por ver qué tan rápido el enfoque de @ cffk converge en una solución, así que escribí una prueba usando arcobjects, que produjo esta salida. Las distancias son en metros:
Aquí está el código fuente. (Editar) Cambió FindCircleCenter para manejar intersecciones (puntos centrales) que caen del borde de la proyección azimutal:
También hay un enfoque alternativo en la edición de junio de 2013 de MSDN Magazine, Amoeba Method Optimization usando C # .
Editar
El código publicado anteriormente convergió a la antípoda en algunos casos. He alterado el código para que produzca esta salida para los puntos de prueba de @ cffk.
Aquí está la salida que ahora produce:
Aquí está el código revisado:
Editar
Aquí están los resultados que obtengo con esriSRProjCS_WGS1984N_PoleAziEqui
fuente