I have 4 coordinates on map.

lat1 = radians(47.434756)
lon1 = radians(42.971468)
lat2 = radians(47.542419)
lon2 = radians(42.946061)

I want calculate distance between this.

from math import sin, cos, sqrt, atan2, radians
R = 6373.0

lat1 = radians(47.434756)
lon1 = radians(42.971468)
lat2 = radians(47.542419)
lon2 = radians(42.946061)

dlon = lon2 - lon1
dlat = lat2 - lat1

a = sin(dlat / 2)**2 + cos(lat1) * cos(lat2) * sin(dlon / 2)**2
c = 2 * atan2(sqrt(a), sqrt(1 - a))

distance = R * c/100

print("Result:", distance)

My code say that distance is 0.1212. ArcGIS says 0,11062 (shape length).

But when I measured with a ruler - its shows 9.23km.

How do I correctly calculate distance with ArcPy?

enter image description here


There's no need to use spherical haversine computation when ArcPy already has the full Inverse (aka Reverse) Problem of Geodesy implemented via the Esri Projection Engine.

The key here is to use the arcpy.PointGeometry type (arcpy.Point is a simple helper class to store double precision values), with an explicit SpatialReference parameter, set to your preferred geographic coordinate system. This example uses WGS84:

import arcpy

sr = arcpy.SpatialReference(4326)               # WGS84
sr.setFalseOriginAndUnits(-400,-400,10000000)   # ten millionth of a degree ~= 1.1112 cm resolution

p1 = arcpy.PointGeometry(arcpy.Point(42.971468,47.434756),sr)  # always use spatial_reference!
p2 = arcpy.PointGeometry(arcpy.Point(42.946061,47.542419),sr)
a1,d1 = p1.angleAndDistanceTo(p2,'GEODESIC')
a2,d2 = p2.angleAndDistanceTo(p1,'GEODESIC')

print("Distance = {:8.3f} meters, Bearing1->2 = {:9.5f} degrees".format(d1,a1))
print("Distance = {:8.3f} meters, Bearing2->1 = {:9.5f} degrees".format(d2,a2))

The results are in degrees and meters (unless the SpatialReference is projected, in which case it's in the units of the projection):

Distance = 12122.175 meters, Bearing1->2 =  -9.07843 degrees
Distance = 12122.175 meters, Bearing2->1 = 170.90285 degrees

Note that, even for such close points, the angle of 2->1 is not a full complement of the angle from 1->2, since the path follows the curvature of the Earth.

If the lat and lon are flipped from the values you provided (Point takes x,y parameters, which is lon,lat)...

p3 = arcpy.PointGeometry(arcpy.Point(47.434756,42.971468),sr)
p4 = arcpy.PointGeometry(arcpy.Point(47.542419,42.946061),sr)
a3,d3 = p3.angleAndDistanceTo(p4,'GEODESIC')
a4,d4 = p4.angleAndDistanceTo(p3,'GEODESIC')

print("Distance = {:8.3f} meters, Bearing3->4 = {:9.5f} degrees".format(d3,a3))
print("Distance = {:8.3f} meters, Bearing4->3 = {:9.5f} degrees".format(d4,a4))

then the distance is 9.23 kilometers.

Distance = 9227.107 meters, Bearing3->4 = 107.77524 degrees
Distance = 9227.107 meters, Bearing4->3 = -72.15139 degrees

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.