4

I would like to know how to find a set of coordinates on a small circle some distance from another point. circles

For example, lets say A = (39.73, -104.98) #Denver

B = (39.83, -106.06) #point approximately 50 nautical miles away from A along the great circle track to C

C = (40.75, -111.88) #Salt Lake City, approximately 323 nautical miles away from A

How do I get the set of orange coordinates when I don't know the angle but I do know the cross track distance?

For concreteness, lets say the outer points near B have a cross track distance of +- 30 nautical miles and the inner points have a cross track of +- 15nm.

And the points near C have a cross track of +- 25 nm.

Also I have an initial bearing from A to C of -1.34 #radians

|improve this question
  • you had to mention what software you can use, Qgis, arcgis .. ? – geogeek Dec 20 '13 at 18:10
  • 1
    I am using Qgis but would prefer a python implementation. Or just the math formula. – Neck Beard Dec 20 '13 at 18:14
  • the easiest way is to make a geoprocessing script using python. – geogeek Dec 20 '13 at 18:20
  • but what should be the inputs, they are A , B and C ? – geogeek Dec 20 '13 at 18:21
  • 2
    How much accuracy do you need? Formulas for the sphere are relatively simple and will generally err by a small fraction of 1%. Formulas for a spheroid can be extremely accurate but typically require expansions using a few terms in a power series or trigonometric series. Don't let that discourage you--the code is simple--but the risk of making an error in precomputing and transcribing the coefficients is large and so such solutions should only be undertaken when they can be extensively tested. – whuber♦ Dec 20 '13 at 19:39
1

If the unknown points are on the specified circles, and if you assume a simple circular reference surface (a sphere, not an ellipsoid), then it is easy to get the angle at A from AC to the unknown points:

angle = ctd / scd, where

ctd is your so-called cross-track distance (what I'd call arc distance), and

scd is your small-circle distance away from A

Then you "just" calculate the new coordinates of unknown point, U, say:

coordsU = traverse (coordsA, azimAC + angle, scd), where

coordsA are coordinates at A,

azimAC is azimuth AC, (conveniently, you already have this in radians),

angle and scd are as above,

traverse is the appropriate function for the so-called "direct" problem of spherical trigonometry -- for now, left as a separate exercise for you to research, or already available from some library.

|improve this answer
  • I believe this is offset or perpendicular rather than arc. I am not sure here is the formula. xtd = asin(sin(d1) * sin(theta1 - theta2)). Where d1 is the distance to an orange point ~50 nm and theta1 is orange point bearing theta2 initial bearing. – Neck Beard Dec 20 '13 at 20:00
  • 1
    @NeckBeard -- If you've already got the the orange point bearing, then the xtc distance is irrelevant! Just jump straight to the traverse() problem. Is that your real question? How to "traverse" (that's one of the technical terms) from a known point, given a bearing and distance, to an unknown point, on the sphere? – Martin F Dec 20 '13 at 20:22
  • I don't have the bearing that is the formula I would use to calculate it if I did. – Neck Beard Dec 20 '13 at 22:12
2

There is no one formula to solve this problem on a spheroid (since it's a partial differential equation which can only be solved by iterative means), but the solution is a trivial implementation of the "direct" problem (given a point, a bearing, and a distance, locate a point). The US National Geodetic Survey has a web site with FORTRAN source (I converted an earlier version of this code to 'C' and Java without difficulty).

You'd need to determine if the cross-track distance should be a chord or an arc when calculating the change in bearing (inverse sine of distance/radius vs. a fraction of the 2*pi*radius perimeter).

EXAMPLE:

Assuming arc distance, the bearing changes 17.18873 degrees for 15nm, 34.37746 degrees for 30nm, and 4.43466 degrees for 25nm, so the problem becomes:

name    lat     lon     bearing distance
B+30    39.73   -104.98 -42.40  50
B+15    39.73   -104.98 -59.59  50
B       39.73   -104.98 -76.78  50
B-15    39.73   -104.98 -93.97  50
B-30    39.73   -104.98 -111.15 50
C+25    39.73   -104.98 -72.34  323
C       39.73   -104.98 -76.78  323
C-25    39.73   -104.98 -81.21  323

which solved to:

name    lat     lon 
B+30    40.344  -105.715
B+15    40.148  -105.917
B       39.916  -106.034
B-15    39.667  -106.057
B-30    39.425  -105.983
C+25    41.167  -111.778
C       40.758  -111.883
C-25    40.345  -111.945

which looks right in this plot:

enter image description here

|improve this answer
  • I suggest "Direct", instead of "Inverse", as the problem. The "Inverse" problem in cogo is to determine bearing and distance between two known points. – Martin F Dec 20 '13 at 19:35
  • Doh. Yes, Direct. Changes the URL as well, though – Vince Dec 20 '13 at 19:36
  • @Vince So basically just a while loop to increment/decrement the bearing and check the cross track. – Neck Beard Dec 20 '13 at 19:57
  • @NeckBeard -- That's one way of doing it i suppose, but not the one i'd want to announce to the world. – Martin F Dec 20 '13 at 20:07
  • The answer i gave is essentially an elaboration of 2nd half of Vince's answer. – Martin F Dec 20 '13 at 20:11
1

Geodesic calculations like these can be done in python using http://pypi.python.org/pypi/geographiclib. Assuming you have this installed, the following python computes the answers you want

import sys
sys.path.append("/usr/local/lib/python/site-packages")
from geographiclib.geodesic import Geodesic,Math
NM=1852
MASK=Geodesic.DISTANCE|Geodesic.AZIMUTH|Geodesic.REDUCEDLENGTH
lata=39.73; lona=-104.98
latb=39.83; lonb=-106.06
latc=40.75; lonc=-111.88

gab=Geodesic.WGS84.Inverse(lata,lona,latb,lonb,MASK)
print 'Points around B'
for dist in [-30,-15,0,15,30]:
    g=Geodesic.WGS84.Direct(lata,lona,
                            gab['azi1']+dist*NM/gab['m12']/Math.degree,
                            gab['s12'])
    print g['lat2'],g['lon2']

gac=Geodesic.WGS84.Inverse(lata,lona,latc,lonc,MASK)
print 'Points around C'
for dist in [-25,0,25]:
    g=Geodesic.WGS84.Direct(lata,lona,
                            gac['azi1']+dist*NM/gac['m12']/Math.degree,
                            gac['s12'])
    print g['lat2'],g['lon2']

Running this I get

Points around B
39.3452857009 -105.943379739
39.5808782474 -106.048518139
39.83 -106.06
40.0707028272 -105.975822882
40.2815616574 -105.802535085
Points around C
40.3357882212 -111.941310951
40.75 -111.88
41.1592654552 -111.775837558

Note that the cross track distance should be converted to a change in azimuth using the reduced length, m12, instead of the distance, s12. It doesn't make much difference in this case; but if A and C were separated by a greater distance it would begin to matter.

|improve this answer
  • What is the definition of reduced length? I tried looking for it link but the documentation isn't very good. Just trying to understand the difference between s12 and m12. – Neck Beard Dec 26 '13 at 22:18
  • Consider a geodesic from A to B of length s12; the azimuth at A is alpha1. Now consider a second geodesic of length s12 starting at A with azimuth alpha1 + dalpha1. The second end of this geodesic is separated from B by a distance m12 * dalpha1 where m12 is the reduced length. On a plane, m12 = s12. For a sphere of radius R, m12 = R * sin(s12/R). For an ellipsoid, Helmert gives a formula. – cffk Jan 1 '14 at 12:21
  • See also the Wikipedia article: en.wikipedia.org/wiki/… – cffk Jan 1 '14 at 13:18
  • Finally(!), I've modified my online geodesic calculator to return the reduced length to allow you to see how s12 and m12 differ. See geographiclib.sourceforge.net/cgi-bin/GeodSolve – cffk Jan 1 '14 at 14:38

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.