Mesure de la distance dans Mercator sphérique vs UTM zoné

11

J'ai des points en WGS84 lat / long et j'aimerais mesurer de "petites" distances (moins de 5 km) entre eux.

Je peux utiliser la formule haversine de http://www.movable-type.co.uk/scripts/latlong.html et cela fonctionne très bien.

J'aimerais utiliser des bibliothèques Python Shapely, pour que je puisse faire plus d'opérations que la distance, et parce qu'à l'échelle avec laquelle je travaille, une terre plate est une assez bonne approximation. Pour projeter de manière fiable les coordonnées géographiques en coordonnées cartésiennes, j'utilise celles de Python proj4, mais semble générer des erreurs plus importantes que je ne le souhaiterais.

Si j'utilise la zone UTM locale, j'obtiens des différences entre haversine de quelques mètres, ce qui est bien. Mais je ne veux pas avoir à travailler sur la zone UTM (les points pourraient être mondiaux), j'ai donc essayé avec "Mercator sphérique" mais maintenant les différences entre les distances haversine et projetées sont bien supérieures à 100%. Est-ce vraiment bon pour Mercator sphérique? Tout ce que je veux vraiment, c'est une projection cartésienne réalisable pour deux points à moins de 5 km l'un de l'autre n'importe où dans le monde.

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)

À ce stade, la distance haversine entre eux est de 394 m, et en utilisant la zone utm 27, 395 m. Mais si j'utilise une sphère Mercator, la distance cartésienne est de 904 m, ce qui est loin.

python coordinate-system distance Karl P
la source
La zone UTM est facile à "travailler" en fonction des longitudes. Choisissez une longitude lambda typique, -180 <= lambda <180, et utilisez-la pour calculer le numéro de zone comme Int ((180 + lambda) / 6) +1. Utilisez le signe de la latitude pour choisir entre le nord et le sud. Vous n'avez pas besoin d'utiliser les zones polaires spéciales aux hautes latitudes; en fait, très près d'un poteau, vous pouvez utiliser presque toutes les zones UTM.
whuber

Réponses:

17

Oui, vous obtiendrez ce genre d'erreurs avec une projection globale de Mercator: elle est précise à l'équateur et la distorsion augmente de façon exponentielle avec la latitude éloignée de l'équateur. La distorsion de la distance est exactement de 2 (100%) à 60 degrés de latitude. À vos latitudes de test (64,14 degrés), je calcule une distorsion de 2,294, exactement en accord avec le rapport 904/394 = 2,294. (Plus tôt, j'ai calculé 2,301 mais c'était basé sur une sphère, pas sur l'ellipsoïde WGS84. La différence (de 0,3%) nous donne une idée de la précision que vous pourriez gagner à utiliser une projection à base d'ellipsoïde par rapport à la formule Haversine basée sur une sphère. )

Il n'y a rien de tel qu'une projection globale qui donne des distances très précises partout. C'est une des raisons pour lesquelles le système de zone UTM est utilisé!

Une solution consiste à utiliser la géométrie sphérique pour tous vos calculs, mais vous l'avez rejetée (ce qui est raisonnable si vous allez effectuer des opérations complexes, mais la décision mérite d'être réexaminée).

Une autre solution consiste à adapter la projection aux points comparés . Par exemple, vous pouvez utiliser en toute sécurité un Mercator transverse (comme dans le système UTM) avec un méridien situé près du centre de la région d'intérêt. Déplacer le méridien est une chose simple à faire: il suffit de soustraire la longitude du méridien de toutes les longitudes et d'utiliser une seule projection TM centrée sur le méridien principal (avec un facteur d'échelle de 1, plutôt que le 0,9996 du système UTM). Pour votre travail, cela aura tendance à être plusprécis que l'utilisation de l'UTM lui-même. Il donnera des angles corrects (TM est conforme) et sera remarquablement précis pour les points séparés par seulement quelques dizaines de kilomètres: attendez-vous à une précision supérieure à six chiffres. En fait, je serais enclin à attribuer toute petite différence entre ces distances TM adaptées et les distances Haversine à la différence entre l'ellipsoïde (utilisé pour la projection TM) et la sphère (utilisée par Haversine), plutôt qu'à la distorsion dans le projection.

whuber
la source
Cela semble assez parfait, je suppose que je dois alors créer ma propre chaîne d'initialisation pour proj4, plutôt que de pouvoir utiliser l'une des chaînes EPSG existantes?
Karl P
1
+1 adapter la projection aux points. Je préfère la plaque transversale carrée à la Mercator transversale, mais sur des surfaces suffisamment petites ("à grande échelle"), presque toute projection "centrée" près de la région d'intérêt donnera une bonne précision.
David Cary
@David Idée intéressante. Sur la sphère, la plaque transversale Carree (Cassini) sera proche de la formule approximative que j'ai donnée sur gis.stackexchange.com/posts/2964/edit (ce qui pourrait être une solution acceptable ici). Les formules pour TM et TPC sont similaires sur la sphère; sur un ellipsoïde, le TPC est légèrement plus simple. TM est probablement pris en charge par plus de logiciels.
whuber
@Karl Vous pouvez utiliser n'importe laquelle des zones TM si vous le souhaitez. Il suffit de déplacer toutes les longitudes pour qu'un point central de votre région d'intérêt coïncide avec le méridien central de la zone choisie. Multipliez toutes les distances par 1 / 0,9996 (et multipliez toutes les zones par le carré de ce facteur), ne changez aucun angle ou relèvement, et - si vos calculs produisent de nouveaux points - redéfinissez simplement leurs longitudes sur le système de coordonnées d'origine .
whuber