Calculating a destination point based on an origin, distance and bearing?

1073
9
Jump to solution
04-25-2022 01:36 PM
EricEagle
Occasional Contributor III

I'm porting a tool from R to be run in a managed ArcGIS Pro instance (no access to R).  In the script, destination points are calculated from a known origin based on bearing and distance.

The script relies on the geosphere destPoint() method to do it, and is invoked like this:

x = destPoint(origin, bearing, distance)

x is a location object with .latitude and .longitude properties that return as floats for use in further processing.

This seems like a super straightforward geometric problem but I am struggling to find how to use ArcPy to accomplish it.  For initial testing purposes I gave up and just pip installed GeoPy, using geopy.distance.distance(unit=distance).destination((origin.Y, origin.X), bearing) and it works beautifully but I do not like delivering stuff with janky 3rd party dependencies.

 

How do I do this in ArcPy without re-writing half of GeoPy?

0 Kudos
1 Solution

Accepted Solutions
DanPatterson
MVP Esteemed Contributor
import arcpy
pnt = arcpy.Point(300000, 5000000)
SR = "NAD 1983 CSRS MTM  9"
pnt = arcpy.Point(300000, 5000000)
p_geom = arcpy.PointGeometry(pnt, spatial_reference=SR)

# -- make a new p_geom at 45 degrees from N at 1000m
new_p = p_geom.pointFromAngleAndDistance(angle=45.0, distance=1000, method='GEODESIC')
new_p[0]
 <Point (300707.56992333446, 5000706.502158922, #, #)>

# --- lots more
dir(new_p)
['JSON', 'WKB', 'WKT', '__add__', '__class__', '__cmp__', '__delattr__', 
  '__dict__', '__dir__', '__doc__', '__eq__', '__format__',
...snip
 'angleAndDistanceTo', 'area', 'boundary', 'buffer', 'centroid', 'clip',
 'contains', 'convexHull', 'crosses', 'cut', 'densify', 'difference',
 'disjoint', 'distanceTo', 'equals', 'extent', 'firstPoint', 'generalize',
 'getArea', 'getGeohash', 'getLength', 'getPart', 'hasCurves',
 'hullRectangle', 'intersect', 'isMultipart', 'labelPoint', 'lastPoint',
 'length', 'length3D', 'measureOnLine', 'overlaps', 'partCount',
 'pointCount', 'pointFromAngleAndDistance', 'positionAlongLine', 'projectAs',
 'queryPointAndDistance', 'segmentAlongLine', 'snapToLine',
 'spatialReference', 'symmetricDifference', 'toCoordString', 'touches',
 'trueCentroid', 'type', 'union', 'within']

... sort of retired...

View solution in original post

9 Replies
by Anonymous User
Not applicable

You might try the Bearing Distance to Line tool.

How do you know GeoPy is janky? Did you test it against your R tool and get different results?

Another great option is to use geographiclib.

0 Kudos
EricEagle
Occasional Contributor III

Sorry!  I don’t mean GeoPy is janky only the dependency management inside a Python toolbox for ArcGIS Pro is janky with any third party dependency.

I’d like to avoid the Bearing Distance to Line tool because it is IO bound (needs output to disk) and I don’t need lines and all the overhead, just points derived within memory space.

0 Kudos
by Anonymous User
Not applicable

No apologies needed - It's always worth testing and reporting problems or discontinuity between different solutions. Totally understand the need for tool and environment portability.

It looks like Dan's solution is what you're looking for. For posterity it's worth nothing you should be able to write the Bearing Distance to Line output to memory.

https://pro.arcgis.com/en/pro-app/latest/help/analysis/geoprocessing/basics/the-in-memory-workspace.... 

 

0 Kudos
EricEagle
Occasional Contributor III

Thanks!  I’ll investigate that more.

0 Kudos
JamesHood
Occasional Contributor

From a point feature class as input: 

import arcpy

# samplept is a point feature class with xcoord, ycoord fields as double, calculated geometry to wgs84,
#also added dist field, and bearing fields both as double, both field calculated with the distance and bearing values.
 
input_table = r"C:\Data\Points\points_data.gdb\samplept"
output_line = r"C:\Data\Points\points_data.gdb\bearingLine"

# BearingDistanceToLine
#https://pro.arcgis.com/en/pro-app/latest/tool-reference/data-management/bearing-distance-to-line.htm
arcpy.BearingDistanceToLine_management(input_table, output_line, 'xcoord', 'ycoord', 'dist', 
                                       'KILOMETERS', 'bearing', 'DEGREES','GEODESIC')

#Feature-vertices-to-points.htm
#https://pro.arcgis.com/en/pro-app/2.8/tool-reference/data-management/feature-vertices-to-points.htm
out_point = r"C:\Data\Points\points_data.gdb\bearingPoints"
arcpy.management.FeatureVerticesToPoints(output_line, out_point, "END")       


output  is point feature class also.  

EricEagle
Occasional Contributor III

This is great but I’d like to do it in memory like with GeoPy without incurring overhead of writing feature classes to disk.

0 Kudos
JamesHood
Occasional Contributor

Just for reference I modified my approach to no longer output  intermediates to hard disk. 

import arcpy

# samplept is a point feature class with xcoord, ycoord fields as double, calculated geometry to wgs84,
#also added dist field, and bearing fields both as double, both field calculated with the distance and bearing values.
 
input_table = r"C:\Data\Points\points_data.gdb\samplept"
output_line = r"in_memory\bearingLine"

# BearingDistanceToLine
#https://pro.arcgis.com/en/pro-app/latest/tool-reference/data-management/bearing-distance-to-line.htm
arcpy.BearingDistanceToLine_management(input_table, output_line, 'xcoord', 'ycoord', 'dist', 
                                       'KILOMETERS', 'bearing', 'DEGREES','GEODESIC')

#Feature-vertices-to-points.htm
#https://pro.arcgis.com/en/pro-app/2.8/tool-reference/data-management/feature-vertices-to-points.htm
out_point = r"in_memory\bearingPoints"
arcpy.management.FeatureVerticesToPoints(output_line, out_point, "END")
DanPatterson
MVP Esteemed Contributor
import arcpy
pnt = arcpy.Point(300000, 5000000)
SR = "NAD 1983 CSRS MTM  9"
pnt = arcpy.Point(300000, 5000000)
p_geom = arcpy.PointGeometry(pnt, spatial_reference=SR)

# -- make a new p_geom at 45 degrees from N at 1000m
new_p = p_geom.pointFromAngleAndDistance(angle=45.0, distance=1000, method='GEODESIC')
new_p[0]
 <Point (300707.56992333446, 5000706.502158922, #, #)>

# --- lots more
dir(new_p)
['JSON', 'WKB', 'WKT', '__add__', '__class__', '__cmp__', '__delattr__', 
  '__dict__', '__dir__', '__doc__', '__eq__', '__format__',
...snip
 'angleAndDistanceTo', 'area', 'boundary', 'buffer', 'centroid', 'clip',
 'contains', 'convexHull', 'crosses', 'cut', 'densify', 'difference',
 'disjoint', 'distanceTo', 'equals', 'extent', 'firstPoint', 'generalize',
 'getArea', 'getGeohash', 'getLength', 'getPart', 'hasCurves',
 'hullRectangle', 'intersect', 'isMultipart', 'labelPoint', 'lastPoint',
 'length', 'length3D', 'measureOnLine', 'overlaps', 'partCount',
 'pointCount', 'pointFromAngleAndDistance', 'positionAlongLine', 'projectAs',
 'queryPointAndDistance', 'segmentAlongLine', 'snapToLine',
 'spatialReference', 'symmetricDifference', 'toCoordString', 'touches',
 'trueCentroid', 'type', 'union', 'within']

... sort of retired...
EricEagle
Occasional Contributor III

Nice!  This looks way closer to what I’m looking for.  Will give it a shot and report back, thanks Dan

0 Kudos