GeometricalGeodesy

A Mathematica Package for High-Accuracy Geodetic Computations

Thomas H. Meyer

Dept. Natural Res. Mgt. & Engr.
University of Connecticut

Preliminaries...

myPath=Switch[$MachineID,
"6139-69034-18747","C:/Documents and Settings/tmeyer/My Documents/presentations/Conferences/Wolfram Tech 07/"
];

If[Head[exDataPoints]===Symbol,exDataPoints=Import[ToFileName[$GeometricalGeodesyDataDirectory,"worldpoints.txt"],"CSV"];thinDataPoints=Select[exDataPoints,RandomInteger[{1,5}]==1&];]

intersectLinePlane[plane:{{{x1_,y1_,z1_},{x2_,y2_,z2_},{x3_,y3_,z3_}}},line:{{x4_,y4_,z4_},{x5_,y5_,z5_}}]:=-((Det[(({{1, 1, 1, 1}, {x1, x2, x3, x4}, {y1, y2, y3, y4}, {z1, z2, z3, z4}}))])/(Det[(({{1, 1, 1, 0}, {x1, x2, x3, x5-x4}, {y1, y2, y3, y5-y4}, {z1, z2, z3, z5-z4}}))]))

tangentPlane:=Graphics3D[{Yellow,Opacity[0.3],Polygon[{{-1,-1,1},{-1,1,1},{1,1,1},{1,-1,1}}]}]

normalPlane:=Graphics3D[{Red,Opacity[0.3],Polygon[{{0,0,1},{1,0,1},{1,0,0},{0,0,0}}]}]

earth:=Graphics3D[{Green,Sphere[{{0,0,0}}]}]

zenithLine:=Graphics3D@Line[{{0,0,0},{0,0,1.5}}]

ssPnt:=Graphics3D[Point[{0,0,1}]]
poi:=Graphics3D[{Red,Point[{Sin[40°],0,Cos[40°]}]}]
projpoi:=Graphics3D[{Blue,Point[{Sin[40°],0,1}]}]
truepoi:=Graphics3D[{Green,Point[t{Sin[40°],0,Cos[40°]}]}]
t=intersectLinePlane[{{{0,0,1},{0,1,1},{1,0,1}}},{{0,0,0},{Sin[40°],0,Cos[40°]}}]
rad:=Graphics3D[{Red,Line[{{0,0,0},t{Sin[40°],0,Cos[40°]}}]}]
prp:=Graphics3D[{Blue,Line[{{Sin[40°],0,Cos[40°]},{Sin[40°],0,1}}]}]

flatEarth={tangentPlane,earth,ssPnt,poi,rad,prp,projpoi,truepoi,normalPlane,zenithLine};

Import["C:/Documents and Settings/tmeyer/My Documents/My Pictures/figures/geodesy_bridge.jpg"]

Import["C:/Documents and Settings/tmeyer/My Documents/My Pictures/figures/fender bender.jpg"]

Geodesy...

is the scientific discipline that deals with the measurement and representation of the earth, its gravitational field, and other geodynamic phenomena, such as crustal motion, oceanic tides, and polar motion. (Wikipedia)

Audience

({{Navigation, Surveying, Mapping and Charting}, {Geographic Information 
Systems, Geophysics, Geodynamics}, {Emergency Management, War Fighting, Remote Sensing}})

Why Study Geodesy?

Why Study Geodesy?

Why Study Geodesy?

Why Study Geodesy?

This is what the runway looks like during a white-out snow storm...

Fields of Geodesy

Physical Geodesy...

is primarily concerned with the earth's gravity field for the purposes of creating heighting systems and vertical datums.

Geometrical Geodesy...

is primarily concerned with the mathematics needed for "horizontal" positioning while accounting for the earth's round shape.

Products of Geometrical Geodesy

Geodetic Datums

A geodetic datum is a system by which coordinates are assigned to places on, in, or above the earth.

Positioning

Geodetic coordinates are spherical (latitude and longitude) and linear (ellipsoid height).

Charting, Mapping, Surveying

Geodetic datums underlie all high-accuracy and globally-defined cartographic products.

Navigation

Global navigation depends on well-defined latitude/longitude graticules.

Geodetic Datums

Horizontal Geodetic Datums

This is the system of latitude and longitude upon which all mapping, charting, and navigation activities depend

Vertical Geodetic Datums

Height

Why do we need geodesy?

Geodetic coordinate systems and computations account for the Earth's curvature

What direction is up?

Zenith Angle

The zenith angle ζ is the spatial angle from the local vertical down

o={0,0,0}; (* the origin of the ECEF coordinate system *)
Champaign=MakeBLH[42.12°,-88.26°];
%//PPDMS

Bowring's coordinate system transformation

Champaign⇒xyz

full algebra of the GeometricalGeodesy objects

2(Champaign⇒xyz)

convert the object to a Mathematica vector for display

zenithPoint=2(Champaign⇒xyz)//ToVector

Graphics3D[{Sphere[o,semiMajor[]],Thick,Line[{o,zenithPoint}]}]

Construct a local Cartesian coordinate system on the surface of the reference ellipsoid

lcs=MakeENUSystem[Champaign];
%//P

I can easily express the corners of a plane tangent to the Earth here just by giving the coordinates in the local coordinate system in the obvious way then coverting to XYZ

MakeENU[1000km,1000km,0]⇒{xyz,lcs}

These are the four corners in the global Cartesian system

d=3000;
tangentPlane=
Polygon[
{MakeENU[d km,d km,0]⇒{xyz,lcs}//ToVector,
MakeENU[-d km,d km,0]⇒{xyz,lcs}//ToVector,
MakeENU[-d km,-d km,0]⇒{xyz,lcs}//ToVector,
MakeENU[d km,-d km,0]⇒{xyz,lcs}//ToVector}]

Graphics3D[{
Sphere[o,semiMajor[]],
Thick,Line[{o,zenithPoint}],
Blue,
tangentPlane
}]

Create a line of observation from the zenith down to the horizon. The zenith angle is the spatial angle from the black line to the red

r=5000km;
Manipulate[
Graphics3D[{
Sphere[o,semiMajor[]],
Thick,Line[{o,zenithPoint}],
Blue,
tangentPlane,
Red,
Thick,
Line[{MakeENU[0,0,0]⇒{xyz,lcs}//ToVector,MakeENU[0,r Sin[ζ °],r Cos[ζ °]]⇒{xyz,lcs}//ToVector}]
}],
{ζ,0,90}
]

Plane Surveying

The Forward Problem

Given an origin, an orientation, and a distance, provide the coordinates of the implied point

The Flat Earth Assumption

Instruments, such as total stations, are oriented to the local horizontal in which they observe  zenith angles ζ, horizontal angles β, and straightline distances (Subscript[d, slope]).

Compute horizontal distance d from observed zenith angle ζ and straightline distance (Subscript[d, slope]). Orient the survey with horizontal angles referred to NORTH.

This is nothing more than the conversion from polar to Cartesian coordinates

Import[myPath<>"Surveying.gif"]

The Flat Earth Assumption

Assuming a polar/Cartesian coordinate system implies that the surface over which the computations are performed is planar. In fact, the earth is round and this assumption introduces unacceptable positioning errors over surprisingly short distances.

Show[flatEarth,Boxed->False]

Assume the earth is spherical. The forward point is in green. It is well above the earth but, by construction, lays in the tangent plane thus having a projected height of zero. The height error (red line) grows as the square of the station separation. For the earth, it would be roughly one cm at 350 m.

The East-North-Up Coordinate System

In geodesy, tangent plane projections are known as "East-North-Up" coordinate systems. GeometricalGeodesy has almost all coordinate systems used in geodesy and built-in conversions between them, often implementing several methods for the same conversion. Let's compute this "trading height for distance" with some real numbers.

Define a stand station (local origin of the ENU system) in terms of latitude, longitude and ellipsoid height

ss=MakeBLH[DMS[41,33,29.042723,N],DMS[72,18,55.923732,W],0.]

"pretty print" functions are defined to make looking at the coordinates easier

ss//PPDMS

Define a forward point exactly one mile away at an azimuth of 45°. Use the GG Direct[] function

?GeometricalGeodesy`Computations`*

fs=Direct[ss,45.°,1USMile][[1]]
fs//PPDMS

Note the units. GG has an extensive collection of linear and angular units used in geodesy, mapping, and surveying. All computations are performed in metric linear units and radians, however.

Now, convert the forward station to ENU in the local system. First create the local system. There is a pretty print function for it.

enusys=MakeENUSystem[ss];
enusys//P

The conversion is done with Esc␣= > Esc (note the space after the first escape). Note the negative ellipsoid height

fsENU=fs⇒{enu,enusys}

Use the norm operator to compute the planimetric distance

\[LeftBracketingBar]fsENU,NormOrder->2\[RightBracketingBar]/USFt//P

At 50 miles, the horizontal displacement is 0.0044'. Use P[expr,ndigits] to display w/o scientific notation

P[(50-\[LeftBracketingBar](Direct[ss,45.°,50USMile][[1]])⇒{enu,enusys},NormOrder->2\[RightBracketingBar]/USMile)/USFt,2]

d=1000;
uE[(Direct[ss,0 °,d USMile][[1]])⇒{enu,enusys}]/USFt

At 1000 miles we've lost 35 feet horizontally and the quadratic growth of the error is evident

GraphicsGrid[{{Plot[(d-\[LeftBracketingBar](Direct[ss,0 °,d USMile][[1]])⇒{enu,enusys},NormOrder->2\[RightBracketingBar]/USMile)/USFt,{d,0,1000},Frame->True,FrameLabel->{"Separation on the Ellipsoid (Miles)","Horizontal Position Error (Feet)"}],
Plot[uE[(Direct[ss,0 °,d USMile][[1]])⇒{enu,enusys}]/USMile,{d,0,1000},Frame->True,FrameLabel->{"Separation (Miles)","Vertical Position Error (Miles)"}]
}}
]

It is a function of azimuth, which is due to the eccentricity of the reference ellipsoid

With[{d=1000},
Plot[(d-\[LeftBracketingBar](Direct[ss,α °,d USMile][[1]])⇒{enu,enusys},NormOrder->2\[RightBracketingBar]/USMile)/USFt,{α,0,180},AxesLabel->{"Azimuth (degrees)","Horizontal Position Error (Feet)"}]]

The Inverse Problem

In geodesy, the inverse problem is, given two geodetic points, what is the distance and angle between them? If the earth were spherical, we'd use spherical trigonometry. However, the eccentricity of the earth is better modeled with a spheroid, which complicates the mathematics significantly.

We'll use the same stand and forward stations

ss//PPDMS
fs//PPDMS

Notice that the forward and reverse angles do not sum to 180°. The sum exceeds 180° by 32"5 seconds of arc. This is due to the eccentricity of the reference ellipsoid

{dist,forwardAngle,reverseAngle,convergence}=InverseGeodetic[fs,ss];
Print["distance: ",dist/USMile//P," USMile"];
Print["forwardAngle: ",forwardAngle//PPDMS];
Print["reverseAngle: ",reverseAngle//PPDMS];
Print["angular difference: ",forwardAngle-reverseAngle//PPDMS];
Print["convergence: ",convergence];

The Geometrical Geodesy Package

Units

Reference Ellipsoids

Coordinate Systems

Computations

Cartographic Projections and Standardized Coordinate Systems

Units

All Units

?GeometricalGeodesy`Units`*

Conversions are ratios from the given unit to meters

100USFt

Divide by the conversion to go from meters to that unit

900./GuntersChain

Angular Units

Degrees-Minutes-Decimal Seconds, decimal degrees, mil, circle, milli-arcsecond, gon

DMS[55,1,9.2,S]

DMS[155,1,9.2,West]

6400milNATO

1milNATO//PPDMS

Reference Ellipsoids

$ellipsoidNames

{#,semiMajor[#],semiMinor[#],flattening[#],eSq[#],ϵSq[#]}&/@$ellipsoidNames//TableForm

Radii of Curvature

in the meridian, prime vertical and normal section

etc

meridional arcs, mean radius, authalic sphere radius

Coordinate Systems

Geodetic (BLH)

datum specific latitude, longitude, ellipsoid height

Earth-Centered, Earth-Fixed Cartesian (XYZ)

datum specific "geocentric" x, y, z. The x and y axes define the equatorial plane. The z-axis forms a right-handed system and represents an average location of the rotational axis of the earth during some epoch. These are not truly geocentric. Instead, the orgin is at the location that the datum places the origin of its reference ellipsoid.

Local Horizontal (ENU)

datum specific, tangent plane, Cartesian such that the easting axis matches the direction of the projection of the parallel of the local origin at the origin; the northing axis matches the direction of the projection of the meridian of the local origin into the tangent plane, and the up axis forms a right-handed system, positive directed away from the reference ellipsoid.

Grid

Coordinates resulting from a cartographic projection.

Geodetic Surveying with a Total Station

Given the geodetic coordinates of the stand station and suppose a target is observed at a spatial distance of 130.51 m, a zenith angle of 89°42'12'' and an azimuth of 23°02'35'', what are the geodetic coordinates of the forward station?

ss=MakeBLH[DMS[45,12,15.,S],DMS[98,33,10.,E],1300.]
enusys=MakeENUSystem[ss];

The observations form a vector in the ENU coordinate system

d=130.51;
α=DMS[23,2,35.];
ζ=DMS[89,42,12.];
h=d Sin[ζ]

Create the implied vector in the local horizontal system. It's length is as given, as expected

v=MakeENU[h Sin[α],h Cos[α],d Cos[ζ]]
\[LeftBracketingBar]v\[RightBracketingBar]

ENU coordinates of forward point: origin (0,0,0) + v = v; nothing to do.
Next, transform ENU to XYZ to BLH

v⇒{xyz,enusys}

v⇒{xyz,enusys}⇒blh

fs=v⇒{xyz,enusys}⇒blh;
fs//PPDMS

{dist,forwardAngle,reverseAngle,convergence}=InverseGeodetic[ss,fs];
Print["distance: ",dist//P," meter"];
Print["forwardAngle: ",forwardAngle//PPDMS];
Print["reverseAngle: ",reverseAngle//PPDMS];
Print["angular difference: ",forwardAngle-reverseAngle//PPDMS];
Print["convergence: ",convergence];

This can be done very easily in one step. This simple function takes total station observations and produces the geodetic latitude, longitude and ellipsoid height of the observed station. The forward station can serve as a local origin, thus implementing the traverse.

GeodeticTraverse[d_,ζ_,α_,enusys_]:=
With[{h=d Sin[ζ]},
MakeENU[h Sin[α],h Cos[α],d Cos[ζ]]⇒{xyz,enusys}⇒blh
]
GeodeticTraverse[d,ζ,α,enusys]//PPDMS

Computations

?GeometricalGeodesy`Computations`*

liftoffDL[]//PPDMS

?GeometricalGeodesy`Area`*

?GeometricalGeodesy`Clairaut`*

Reductions

?GeometricalGeodesy`Reductions`*

Cartographic Projections

General

Albers Equal Area Conic, Gnomic, Lamberts Conformal Conic (Snyder & NGS), Mercator, Oblique Mercator (Snyder & NGS), Transverse Mercator (Snyder, NGS, & Rollins)

Standard Coordinate Systems

State Plane Coordinate System (NAD 27 & NAD 83), Univeral Transverse Mercator

ps=MakeProjectionSystem[SPCS,"0600"];
ps//PP

SetProjectionSystem[ps];

ss=MakeBLH[42°,-72°];
ProjectBLH[ss]//P

ProjectBLH[ss,IncludeScaleFactor->True]//P

ProjectBLH[ss,IncludeConvergence->True]//P

{{e,n},k,γ}=ProjectBLH[ss,IncludeScaleFactor->True,IncludeConvergence->True];
γ//PPDMS

Here's the LCC graticule. Build two lists of graphics. The first holds meridians (red); the second holds parallels (blue)

SetProjectionSystem[
MakeProjectionSystem[
LambertConformalConic,
Standards->{33°,45°},
CentralParallel->(23°),
CentralMeridian->(-96°)]];
maxLat=80;maxLon=-6;
minLat=0;minLon=-186;
meridians=(ParametricPlot[ProjectBLH[ϕ °,#1 °],{ϕ,minLat,maxLat},PlotStyle->{RGBColor[1,0,0]}]&)/@Range[minLon,maxLon,10.];
parallels=(ParametricPlot[ProjectBLH[#1 °,λ °],{λ,minLon,maxLon},PlotStyle->{RGBColor[0,0,1]}]&)/@Range[minLat,maxLat,10.];
Show[meridians,parallels,PlotRange->All]

Show[meridians,parallels,ListPlot[(ProjectBLH[#1[[2]] °,#1[[3]] °]&)/@Select[thinDataPoints,#1[[2]]<maxLat&&#1[[2]]>minLat&&#1[[3]]<maxLon&&#1[[3]]>minLon&],PlotStyle->PointSize[0.001]],PlotRange->All]

Future work

datum transformations

support for physical geodesy

Thank you

Thomas H. Meyer
thomas.meyer@uconn.edu