Algoritmo de intersección que maneja correctamente el meridiano de 180 ° y los polos?

¿Existe un algoritmo de intersección bien conocido que maneje correctamente 180 meridianos y polos?

Por ejemplo, supongamos que tenemos una lista de valores de latitud y longitud que representan la Antártida. Digamos que también tenemos un polígono simple que representa un avión. Queremos saber si el avión está sobre la Antártida.

La intersección genérica del polígono 2D no funciona en este caso si simplemente usa la latitud para y y la longitud para x, porque el sistema de coordenadas planas tendrá bordes en -180 y 180 de longitud, y -90 y 90 de latitud. La Antártida se saldrá de la página en tres de los bordes.

Los SIG (como campo) no han funcionado demasiado cuando se trata de lidiar realmente con la superficie del globo.

Por ejemplo, su problema no está completamente definido. A diferencia de 2D, donde sabemos que los bordes de un polígono están compuestos de líneas rectas, ¿qué son en el globo? Los arcos de grandes círculos, que minimizan la distancia entre vértices, son una buena opción, pero no la única. Las líneas rectas (que viajan por debajo de la superficie del globo) son otra opción, por ejemplo. (Las líneas de rumbo también son una opción, pero potencialmente una tonta, especialmente cerca de los polos ...) Aquí hay un excelente enlace Problemas de interpretación de borde WRT: http://blog.opengeo.org/2010/08/10/shape-of -a-polígono /

Otro gran problema con el enfoque de proyecto a 2d y luego intersección, además de las singularidades que otros han mencionado, es que cualquier punto de intersección recién creado (donde se cruzan los bordes del polígono) estará fuera de lugar y dependerá de la proyección utilizada. (Creo que de aquí proviene la recomendación para la densificación: al agregar una tonelada de vértices intermedios, se obtiene un error reducido de la proyección del centro de los bordes del polígono).

Suponiendo que no desea hacer todos los compromisos y soluciones alternativas de proyección a 2D y está buscando pensar y codificar algo usted mismo, he hecho un poco de código prototipo (¡y solo prototipo!) Para un trato con un cliente con este.

Aquí hay un bosquejo del enfoque. Necesitará saber qué es un vector , así como el significado de los productos de punto y cruz . (Advertencia: Wikipedia es conveniente para enlaces rápidos, pero suficiente si es su primera introducción a los temas. Un buen tutorial de gráficos en 3D lo ayudará).

• Representar un punto en la esfera por una unidad vectorial cartesiana en 3D. El punto en la superficie de la tierra es donde, si extendieras el vector en un rayo infinitamente largo, se cruzaría con la superficie de la tierra.
• Representar grandes círculos por un plano a través del origen. (En 3D, un solo vector unitario es suficiente para definir un plano a través del origen; es lo normal al plano). El gran círculo es la intersección de todo el plano con la superficie de la tierra.
• Puede encontrar los puntos de intersección de dos grandes círculos intersectando sus dos planos.
• Defina un arco de un gran círculo por dos puntos. El vector que define el gran círculo normal es el producto cruzado de los vectores de los puntos inicial y final.
• Para determinar si un punto que sabemos que está en el gran círculo se encuentra dentro del arco, cree dos planos de manera que: sean perpendiculares al plano del gran círculo, uno contiene el punto de inicio y el otro contiene el punto final, y están orientados uno frente al otro. Entonces el punto está en el arco si se encuentra en el "interior" de ambos planos. (Para ayudar a visualizar: usted ha creado la sonrisa de pac-man mientras muerde el arco. Si el punto de prueba está entre sus mandíbulas, entonces yace en el arco, como ya sabemos que está en el gran círculo).
• Para determinar si dos arcos se cruzan: encuentre los dos puntos de intersección de sus grandes círculos correspondientes, luego pruebe cada punto para ver si se encuentra dentro de ambos arcos.
• Una definición menos ambigua: un polígono es una colección de puntos, cada uno de los cuales está conectado por aristas que consisten en arcos de grandes círculos. Los puntos están ordenados de tal manera que si caminas a lo largo de la superficie de la tierra a lo largo de los bordes del polígono, el "interior" del polígono estará a tu izquierda. Dejemos los polígonos complejos (islas, agujeros y auto-intersecciones) por ahora.
• Puede saber si un punto está en el lado derecho o izquierdo de un plano mediante el signo del producto escalar de sus vectores correspondientes. (Esto es equivalente a, a medida que viaja alrededor del gran círculo, si el punto en la superficie de la esfera se encuentra a su izquierda o derecha).
• Una prueba precisa para determinar si un punto está dentro de un polígono: ¿se encuentra en el lado izquierdo de todos los bordes?
• Ahora tenemos la capacidad de probar si un punto está dentro de un polígono y determinar los cruces de bordes: ¡los ingredientes para las intersecciones polígono-polígono! El margen de este comentario es demasiado pequeño para escribir un algoritmo completo, pero los pasos básicos son: (a) encontrar todas las intersecciones, luego (b) caminar bordes, alternando en qué polígono estás caminando cuando te encuentras con puntos de intersección.
• Una vez que todo lo anterior esté funcionando, comience a pensar en estrategias de indexación para hacerlo más rápido, ya que el contorno del punto en el polígono I es O (n) en el número de bordes, y la intersección del polígono O (m * n) en número de bordes

Pheew

Este enfoque tiene algunas grandes ventajas: todas las operaciones anteriores se reducen a solo multiplicaciones y sumas. (Después de convertir los datos a esta representación: p. Ej., Coordinado lat / long requiere algo de trigonometría para obtener un vector XYZ cartesiano). No hay singularidades en los polos o en el sistema de coordenadas, y no hay muchos casos especiales de los que preocuparse (ejemplo, caso especial : bordes superpuestos paralelos).

Al echar un vistazo al código, parece que el paquete Spheres que alguien más ha vinculado sigue algunos de estos enfoques, aunque también parece un poco cocido.

PostGIS también puede usar un enfoque similar sobre su tipo de datos de geografía , pero no he mirado debajo del capó. Sé que para la indexación espacial, al menos, usan un árbol R sobre cartesiano 3D.

(Nota: esta respuesta se hizo lo suficientemente larga como para que probablemente la edite en una publicación de blog ... ¡Comentarios / comentarios muy bienvenidos!)

De hecho, estoy buscando una solución general que funcione para muchos polígonos grandes.

La mayoría de las respuestas que he leído hasta ahora se centran naturalmente en encontrar una proyección adecuada para sus características. Esto puede ser contrario a la intuición, pero aquí debe considerar dónde una proyección no funciona, no dónde funciona bien.

Algunas proyecciones no funcionan solo en un punto: estas se encuentran entre las llamadas proyecciones "globales". Un buen ejemplo es el equidistante azimutal: el único punto que no puede proyectar es el diametralmente opuesto a su punto de origen. (En realidad, no es estrictamente cierto que no se puede proyectar ese último punto; estoy pasando por alto un problema técnico. Deje que sea suficiente que la proyección experimente problemas graves en el polo opuesto).

Por lo tanto, si puede encontrar un único punto P fuera de la unión de todos los polígonos con los que planea trabajar, entonces todo lo que necesita hacer es trabajar dentro de una proyección global que falla solo en este punto único. ( Por ejemplo , use un acimutal equidistante oblicuo centrado en el punto diametralmente opuesto a P ).

(Puede encontrar dicho punto automáticamente. Por ejemplo, cree un polígono que cubra la tierra, amortigüe ligeramente la unión de todas sus características, elimine este búfer del polígono de tierra y seleccione cualquier punto arbitrario dentro de lo que queda. El motivo del almacenamiento en búfer es que es mejor mantenerse alejado de los puntos singulares donde falla una proyección, cuanto más mejor, este es un problema de precisión computacional).

Hay muchas proyecciones globales con propiedades deseables, que incluyen equidistante, conforme ( p . Ej. , Estereográfica) e igual área. Por lo tanto, elegir una proyección global para su análisis no es una limitación severa.

¡La principal limitación es que el software más popular que existe, ArcGIS, no es compatible con muchas versiones oblicuas de proyecciones globales! :-(

¡Puedes resolver esto con proyecciones!

El hecho de que los "extremos" de las escalas estén en -180 y +180, no significa que tenga que considerar el mapa con esos extremos.

Utilice EPSG: 32761 ( http://spatialreference.org/ref/epsg/32761/ )

Convierta sus valores de latitud y longitud en coordenadas en el espacio UPS, y luego puede usar cálculos de geometría regulares para determinar si su avión está sobre la Antártida.

Aquí encontrará más información sobre el sistema de coordenadas estereográficas polares universales: http://en.wikipedia.org/wiki/Universal_Polar_Stereographic_coordinate_system

Si solo está tratando de determinar un algoritmo de intersección para polígonos gigantes, entonces trate la tierra como un esferoide (o incluso una esfera) y use geometría esférica para analizar sus polígonos.

Encontré una biblioteca llamada Spheres que tiene GPL, por lo que podría revisar su código para ver con precisión qué algoritmos usan para la intersección de arco y los "puntos dentro de un polígono esférico".

Sin embargo, cuando desee DIBUJAR estos polígonos en una superficie plana (como el papel o la pantalla), deberá elegir una proyección , o al menos un estilo de proyección.

Si bien hay problemas con el uso de coordenadas cartesianas, no veo el problema en el ejemplo que ha proporcionado. Esto es lo que me parece, con un polígono que cubre la Antártida delimitado por los bordes del espacio de coordenadas:

_{(fuente: geohack.net )}

La intersección entre el plano y el continente se puede calcular de forma segura sin realizar transformaciones complejas.

Una solución a este problema es utilizar algún tipo de sistema de coordenadas tridimensional para la Tierra en lugar de un sistema de coordenadas bidimensional. Si a los polígonos 2D se les da una gran cantidad virtual de profundidad, como la mitad del centro de la tierra, entonces los polígonos no necesitan ser modelados como curvas en la superficie de la tierra. Debería ser suficiente para que solo los vértices estén en la superficie.

Una posible solución es cortar el polígono en subpolígonos en todos los bordes (180 meridianos y polos). El software podría incluso mantener algún tipo de referencia entre el polígono original y los subpolígonos.

Tome el centroide del avión y úselo como punto de tangencia para una proyección azimutal . También necesitaría densificar las geometrías antes de proyectarlas para casos como el área de la tierra que actualmente ve la luz del día.