1

I'm calculating distances between points in arcpy or ArcGIS API for Python, and am probably doing something wrong. Hopefully, someone can show me what. Take this as an interactive session with ArcGIS Pro's Python interpreter.

I'm starting with:

earthRadiusInKm = 6378.137

sr = arcgis.geometry.SpatialReference(wkid=4326)
type(sr)    # arcgis.geometry._types.SpatialReference
sr          # '{ "wkid": 4326 }'

# Looking good so far 
p1 = arcgis.geometry.Point (x=0, y=0, spatialReference=sr)
p1    # { "x": 0, "y": 0, "spatialReference": { "wkid": 4326 } }
p2 = arcgis.geometry.Point (x=1, y=1, spatialReference=sr)

p1.distance_to(p2)           # 1.4142135623...
math.radians(p1.distance_to(p2)) * earthRadiusInKm # 157.4295
# which, the cos(theta) being ~= theta for small angles, looks correct

Going to skip a few specifics (since I can't cut and paste from where I did my testing)

The long and short of it is, though, that, no matter what coordinates I put in (including, say, just off from the poles, or near the Date Line), Point.distance_to() is returning the Pythagorean distance between those two points. So, it's very close on small angular differences, except for the edge cases above, and gets further off the larger the degree difference is.

I assume that's because the SpatialReference is somehow wrongly-formed, but I have no idea how it is (the documentation on SpatialReference is lacking some specifics. I have no clue what sort of iterable it would take as an unlabeled argument, for instance). I also tried playing around with different spatial reference arguments (just passing an integer of 4326 in, say), but always got the same results.

Can anyone point me in the right direction?

I have tagged this with ArcPy, even though my examples are written in ArcGIS API for Python, since it appears that ArcGIS API for Python uses ArcPy internally (I've always gotten 0 from distance_to() when ArcPy was not already imported).

1

I have a vague memory of reading somewhere that the distance to only does euclidean distances, but I can't find a reference to it.

Also you might actually try using PointGeometry, which accepts the spatial reference when you create the object.

https://pro.arcgis.com/en/pro-app/arcpy/classes/point.htm

https://pro.arcgis.com/en/pro-app/arcpy/classes/pointgeometry.htm

Also, the angleAndDistanceTo function allows you to pick GEODESIC as the method which would give you the distance calcualtion for a spherical (geographic) coordinate system.

In arcpy, this would be like the example

sr = arcpy.SpatialReference(4326)
point = arcpy.Point(0,0)
ptGeometry = arcpy.PointGeometry(point, sr)
point2 = arcpy.Point(1,1)
ptGeometry2 = arcpy.PointGeometry(point2, sr)
ptGeometry.angleAndDistanceTo(ptGeometry2,"GEODESIC")

Edit: Here is the reference for the ArcGIS API for Python and the angle_distance_to function: https://developers.arcgis.com/python/api-reference/arcgis.geometry.html#geometry

|improve this answer|||||
  • Thanks. ArcGIS API has Geometry.angle_distance_to (other_geometry, method="GEODESIC"), which, I expect, calls PointGeometry.angleAndDistanceTo() under the covers. The documentation is vague almost to the point of uselessness, but as near as I can tell from experimenting, it returns a tuple with bearing from point1 to point2, and geodesic distance in meters. – DaveWork Oct 14 '19 at 16:56

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.