Medición de distancia en Mercator esférico vs UTM zonificado

11

Tengo puntos en WGS84 lat / long y me gustaría medir distancias "pequeñas" (menos de 5 km) entre ellos.

Puedo usar la fórmula de Haversine de http://www.movable-type.co.uk/scripts/latlong.html y funciona muy bien.

Sin embargo, me gustaría usar las bibliotecas Python Shapely, para poder hacer más operaciones que solo la distancia, y porque en la escala con la que estoy trabajando, una tierra plana es una aproximación lo suficientemente buena. Para proyectar de manera confiable las coordenadas geográficas a un coord cartesiano, estoy usando Python proj4, pero parece tener errores más grandes de lo que me gustaría.

Si uso la zona UTM local, obtengo diferencias entre la haversina de un par de metros, lo cual está bien. Pero no quiero tener que resolver la zona UTM (los puntos podrían ser mundiales), así que probé con "Mercator esférico", pero ahora las diferencias entre la distancia de Haversine y la proyectada están por encima del 100%. ¿Es esto realmente correcto para Mercator esférico? Todo lo que realmente quiero es una proyección cartesiana viable para dos puntos a menos de 5 km de distancia en cualquier parte del mundo.

from shapely.geometry import Point
from pyproj import Proj

proj = Proj(proj='utm',zone=27,ellps='WGS84')
#proj = Proj(init="epsg:3785")  # spherical mercator, should work anywhere...

point1_geo = (-21.9309694, 64.1455718)
point2_geo = (-21.9372481, 64.1478206)
point1 = proj(point1_geo[0], point1_geo[1])
point2 = proj(point2_geo[0], point2_geo[1])

point1_cart = Point(point1)
point2_cart = Point(point2)

print "p1-p2 (haversine)", hdistance(point1_geo, point2_geo)
print "p1-p2 (cartesian)", point1_cart.distance(point2_cart)

En este punto, la distancia haversina entre ellos es de 394m, y utilizando la zona utm 27, 395m. Pero si uso Mercator esférico, la distancia cartesiana es de 904 m, que está muy lejos.

python coordinate-system distance Karl P
fuente
La zona UTM es fácil de "resolver" en función de las longitudes. Elija una longitud típica lambda, -180 <= lambda <180, y úsela para calcular el número de zona como Int ((180 + lambda) / 6) +1. Usa el signo de la latitud para decidir entre el norte y el sur. No necesita usar las zonas polares especiales en latitudes altas; de hecho, realmente cerca de un poste puedes usar casi cualquier zona UTM.
whuber

Respuestas:

17

Sí, obtendrá este tipo de errores con una proyección global de Mercator: es precisa en el ecuador y la distorsión aumenta exponencialmente con la latitud del ecuador. La distorsión de la distancia es exactamente 2 (100%) a 60 grados de latitud. En sus latitudes de prueba (64.14 grados) calculo una distorsión de 2.294, coincidiendo exactamente con la relación 904/394 = 2.294. (Anteriormente calculé 2.301 pero eso se basó en una esfera, no en el elipsoide WGS84. La diferencia (del 0.3%) nos da una idea de la precisión que puede obtener al usar una proyección basada en elipsoide versus la fórmula de Haversine basada en la esfera. )

No existe una proyección global que produzca distancias muy precisas en todas partes. ¡Esa es una razón por la que se utiliza el sistema de zona UTM!

Una solución es usar geometría esférica para todos sus cálculos, pero lo ha rechazado (lo cual es razonable si va a realizar operaciones complejas, pero la decisión podría valer la pena).

Otra solución es adaptar la proyección a los puntos que se comparan . Por ejemplo, podría usar de forma segura un Mercator transversal (como en el sistema UTM) con un meridiano ubicado cerca del centro de la región de interés. Mover el meridiano es algo simple: solo resta la longitud del meridiano de todas las longitudes y usa una sola proyección TM centrada en el meridiano principal (con un factor de escala de 1, en lugar del 0.9996 del sistema UTM). Para su trabajo, esto tenderá a ser máspreciso que usar UTM en sí. Dará ángulos correctos (TM es conforme) y será notablemente preciso para puntos separados por solo unas pocas decenas de kilómetros: espere una precisión mejor que seis dígitos. De hecho, me inclinaría a atribuir cualquier pequeña diferencia entre estas distancias TM adaptadas y las distancias de Haversine a la diferencia entre el elipsoide (utilizado para la proyección TM) y la esfera (utilizada por Haversine), en lugar de la distorsión en el proyección.

whuber
fuente
Eso suena bastante perfecto, supongo que necesito hacer mi propia cadena de inicio para proj4 en lugar de poder usar cualquiera de las cadenas EPSG existentes.
Karl P
1
+1 adapta la proyección a los puntos. Prefiero la placa transversal carrée sobre la transversal de Mercator, pero en áreas lo suficientemente pequeñas ("gran escala"), casi cualquier proyección "centrada" cerca de la región de interés dará una buena precisión.
David Cary
@David Interesante idea. En la esfera, la placa transversal Carree (Cassini) estará cerca de la fórmula aproximada que proporcioné en gis.stackexchange.com/posts/2964/edit (que podría ser una solución aceptable aquí). Las fórmulas para TM y TPC son similares en la esfera; en un elipsoide, TPC es un poco más simple. TM es probablemente compatible con más software.
whuber
@Karl Puedes usar cualquiera de las zonas TM si lo deseas. Simplemente cambie todas las longitudes para que un punto central en su región de interés coincida con el meridiano central de la zona elegida. Multiplique todas las distancias por 1 / 0.9996 (y multiplique todas las áreas por el cuadrado de este factor), no cambie ningún ángulo ni orientación y, si sus cálculos producen nuevos puntos, simplemente cambie sus longitudes nuevamente al sistema de coordenadas original. .
whuber