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).
Answer
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
No comments:
Post a Comment