Fast Algorithms to Draw the Curve of Alinements, the Great Ellipse, the Normal Section, and the Loxodrome

Thomas H. Meyer
Dept. Natural Resources & Environment
University of Connecticut

Join the Conversation #WolframTechConf

The Direct Problem of Geodesy

Geodesy is a branch of science concerned with spatial quantities, such as distances and directions, that pertain to the entire Earth.

Geometric geodesy has two main problems: the direct problem and the indirect problem. The direct problem is the subject today

Given

the geodetic coordinates (λ, ϕ) of a starting point A

an azimuth α (horizontal angle from North, positive Eastwards)

a distance

Calculate the implied point B

The direct problem is used for

land surveying over widely separated places

the construction of regional, national, and global coordinate systems, cadastres, and geographic reference frames

navigation

I’m here going in direction α at a speed of s m/s. Where will I be in ten minutes?

drawing curves in a GIS

In a GIS, polygons and polylines are given as a sequence of points. The direct problem is used to draw the implied curve between them.

The Exact Solution from Geodesics

The direct problem takes a distance for one of its inputs.

The distance is to be minimal, so the required curve is a geodesic.

Geodesy employs an oblate spheroid (reference ellipsoid) for an Earth model, and computing the geodesic on such a spheroid is complicated.

The geodesic direct problem involves this (incomplete) elliptic integral of the third kind:

P (n, (Subscript[θ, 0]),m)=(∫_0^(Subscript[θ, 0]))((dθ)/((1-n (sin^2)θ)(Sqrt[1-m (sin^2)θ])))

This integral can be used to provide arc length along the geodesic, so the direct problem involves a root-finding operation to find the correct upper limit of integration.

Sometimes, the exact solution is not required and a more efficient solution is preferable.

Common Geodetic Curves

Loxodrome

A loxodrome is a curve of constant azimuth.

Great Ellipse

A great ellipse is a section curve containing A and B and the center of the spheroid.

Normal Section

A normal section is a section curve containing A, the spheroid normal at A, and B.

Curve of Alinements

A curve of alinements is a curve containing A and B such that the curve’s normal vector is coincident with the spheroid’s normal vector everywhere along the curve.

Curve of Alinements

Connect A and B with a straight line segment.
N.B. The segment penetrates through the spheroid

Imagine a sleeve that is free to slide along the line segment

Imagine the sleeve has a riser attached to it that is perpendicular to the line segment

Imagine the riser is kept perpendicular to the spheroid

The curve of alinements is locus of points where the riser penetrates the spheroid from A to b

Curves of alinements are...

the same from A to B as from B to A.

not planar (not sections)

Curve of Alinements: Algorithm

Given A and B in geocentric Cartesian coordinates (XYZ)

v=B-A

p(t)=A+t v

Let γ(t) denote the conversion from geocentric (XYZ) coordinates to geodetic (λ, ϕ, h) coordinates

Key intuition: the points are inside the spheroid, so they have negative h values. 
Setting h to zero sets them on the spheroid in the direction normal to its surface

The curve of alinements are all the points γ(p(t)) for 0≤t≤1

CoaFromXYZ[Axyz_(*XYZ*),Bxyz_(*XYZ*),RE_,t_]:=
ggToXYZ[
ReplacePart[
ggToλϕh[Axyz+t×(Bxyz-Axyz),RE],
3->0],
RE]

Curve of Alinements: Example

A={165.0°,-40.0°,0};
Axyz=ggToXYZ[A,wgs84];
Aλϕh=ggToλϕh[Axyz,wgs84];

B={0.0°,45.0°,0};
Bxyz=ggToXYZ[B,wgs84];
Bλϕh=ggToλϕh[Bxyz,wgs84];

Show[
bluMrbl,
ParametricPlot3D[CoaFromXYZ[Axyz,Bxyz,wgs84,t],{t,0,1},PlotStyle->Red],
GeodesicPlot[Aλϕh,Bλϕh,wgs84]
]

Normal Section

Connect A and B with a straight line segment.
N.B. The segment penetrates through the spheroid

Construct the plane containing the spheroid’s normal vector at A and the connecting segment from A to B.

The normal section is locus of points where the plane intersects the spheroid from A to B

Normal sections are...

generally not the same from A to B as from B to A.

planar (sections)

very nearly the same length as the geodesic (shorter than the curve of alinements)

Normal Section: Algorithm

Given A and B in geocentric Cartesian coordinates (XYZ)

There is a point V on the Z-axis at (0,0,-(Subscript[N, A])) and VA is normal to the spheroid at A.

N=(a/(Sqrt[1-(ϵ^2)(sin^2)ϕ]))

v=B-A

q(t,u)=A+u(p(t) - V)
Any multiple u(p(t) - V) is in the sectioning plane.

Find (SuperStar[u]) such that γ(q(t, (SuperStar[u]))) has h = 0
This is a root-finding problem
a/|p(t)| is a decent starting guess b/c it must be fairly close to the surface of the ellipsoid

The normal section are all the points q(t, (SuperStar[u])) for 0≤t≤1 such that γ(q(t, (SuperStar[u]))) has h = 0

NormalSection[Axyz:{X_,Y_,Z_},Bxyz_(*XYZ*),RE:{a_,ϵSq_},t_]:=
Module[
{
Aλϕh=ggToλϕh[Axyz,RE], (* need ϕ for N *)
V(* point on the Z-axis where normal section plane touches *),
p (* point on line from A to B *),
ustar (* ns parameter *)
},
(* point on the Z-axis where normal section plane touches *)
V={0,0,ggN[Aλϕh[[2]],RE]Sin[Aλϕh[[2]]]-Z};
(* point on line from A to B @ t *)
p=Axyz+t×(Bxyz-Axyz);
ustar=ustar/.FindRoot[ggToλϕh[V+ustar(p-V),RE,"Bowring1985"][[3]],{ustar,a/Norm[p]}];
V+ustar(p-V)
]

Normal Section: Example

Show[
bluMrbl,
ParametricPlot3D[CoaFromXYZ[Axyz,Bxyz,wgs84,t],{t,0,1},PlotStyle->Red],
ParametricPlot3D[NormalSection[Axyz,Bxyz,wgs84,t],{t,0,1},PlotStyle->Blue],
ParametricPlot3D[NormalSection[Bxyz,Axyz,wgs84,t],{t,0,1},PlotStyle->Blue]
]

Great Ellipse

Connect A and B with a straight line segment.
N.B. The segment penetrates through the spheroid

Form the plane containing A, B, and the center (0,0,0).

The great ellipse is locus of points at the intersection of the plane and the spheroid from A to B

Great Ellipses are...

the same from A to B as from B to A.

planar (sections)

very nearly the same length as the geodesic (shorter than the curve of alinements and normal sections)

Great Ellipse: Algorithm

The algorithm is identical to that for the normal section except use (0,0,0) for V.

GreatEllipse[Axyz_(*XYZ*),Bxyz_(*XYZ*),RE:{a_,ϵSq_},t_]:=
Module[
{
V={0,0,0},
p (* point on line from A to B *),
q(* ns parameter *)
},
(* point on line from A to B @ t *)
p=Axyz+t×(Bxyz-Axyz)//N;
q=q/.FindRoot[ggToλϕh[V+q(p-V),RE,"Bowring1985"][[3]],{q,a/Norm[p]}];
V+q(p-V)
]

Great Ellipse: Example

Show[
bluMrbl,
ParametricPlot3D[NormalSection[Axyz,Bxyz,wgs84,t],{t,0,1},PlotStyle->Blue],
ParametricPlot3D[NormalSection[Bxyz,Axyz,wgs84,t],{t,0,1},PlotStyle->Blue],
ParametricPlot3D[GreatEllipse[Bxyz,Axyz,wgs84,t],{t,0,1},PlotStyle->Orange]
]

Loxodromes

Loxodromes are...

spirals

extremely useful in navigation

drawn as straight-lines on Mercator projection maps

the same from A to B as from B to A.

often very much longer than any of the other curves

Loxodromes: Algorithm

The algorithm is to simply draw a straight-line between A and B in the Mercator mapping plane and convert to geodetic. This inverse mapping requires iteration

Loxodrome[Aλϕh:{Aλ_,Aϕ_,Ah_},Bλϕh:{Bλ_,Bϕ_,Bh_},RE:{a_,ϵSq_},t_]:=
Module[
{
ϵ=(Sqrt[ϵSq]),
xa,ya,xb,yb,x,y
},
xa=a Aλ;
xb=a Bλ;
ya=(a/2)Log[((1+Sin[Aϕ])/(1-Sin[Aϕ]))(((((1-ϵ Sin[Aϕ])/(1+ϵ Sin[Aϕ]))))^ϵ)];
yb=(a/2)Log[((1+Sin[Bϕ])/(1-Sin[Bϕ]))(((((1-ϵ Sin[Bϕ])/(1+ϵ Sin[Bϕ]))))^ϵ)];
{x,y}={xa,ya}+t({xb,yb}-{xa,ya});
ggToXYZ[{x/a,FixedPoint[π/2-2ArcTan[(E^(-y/a))(((((1-ϵ Sin[#])/(1+ϵ Sin[#]))))^(ϵ/2))]&,π/2-2ArcTan[(E^(-y/a))]],0},RE]
]

Loxodromes: Example

Show[
bluMrbl,
ParametricPlot3D[CoaFromXYZ[Axyz,Bxyz,wgs84,t],{t,0,1},PlotStyle->Red],
ParametricPlot3D[Loxodrome[Aλϕh,Bλϕh,wgs84,t],{t,0,1},PlotStyle->Green]
]

Questions?

thomas.meyer@uconn.edu

Initialization code

SetDirectory[NotebookDirectory[]]
Needs["ggConversions`"]
Needs["ggCurvatures`"]
Needs["ellipsoids`"]
Needs["ggGeodesics`"]
a=ellipsoidPart["a",WGS84];
ϵSq=ellipsoidPart["eccSq",WGS84];
wgs84={a,ϵSq}

exDataPoints=Import["worldpoints.CSV"];
thinDataPoints=Select[exDataPoints,RandomInteger[{1,5}]==1&];
wp=ggToXYZ[{#[[1]]°,#[[2]]°,0},wgs84]&/@thinDataPoints[[All,{3,2}]];
pts=ListPointPlot3D[wp,PlotStyle->{Tiny,Gray}];

bluMrbl=Show[
ParametricPlot3D[ggToXYZ[{λ,45,0}°,wgs84],{λ,0,360},PlotStyle->Darker[Green]],
ParametricPlot3D[ggToXYZ[{λ,0,0}°,wgs84],{λ,0,360},PlotStyle->Darker[Green]],
ParametricPlot3D[ggToXYZ[{0,ϕ,0}°,wgs84],{ϕ,-90,90},PlotStyle->Orange],
Graphics3D[{LightBlue,Sphere[{0,0,0},6378000]}],
pts,
PlotRange->All,
Boxed->False,Axes->False
];

GeodesicPlot[Aλϕh:{Aλ_,Aϕ_,Ah_},Bλϕh:{Bλ_,Bϕ_,Bh_},RE:{a_,ϵSq_}]:=Module[
{
sc,c,Δs,α1,α2,
plt
},
Off[FindRoot::precw];
{sc,c,Δs,α1,α2}=ggIndirectGeodesic[{Aλϕh,Bλϕh},RE];
ParametricPlot3D[ggToXYZ[ggDirectGeodesic[Aλϕh,{α1,d},RE][[1;;2]],RE],{d,0,Δs},PlotStyle->Green]
]