Geometrical Geodesy: An Example

Needs["GeometricalGeodesy`"]

Needs["Graphics`Graphics3D`"]

This example creates the cover image on the box. It exercises many of the features of the package and demonstrates its ease of use.

The cover image depicts the radius of curvature in the prime vertical (RCPV), emphasizing the fact that it does not intersection the Z-axis at the origin. This is due to the eccentricity of the reference ellipsoid. In order to illustrate the point, it is necessary to create a custom datum and ellipsoid, one that has very exaggerated eccentricity compared to any reference ellipsoid used in geodesy, because the point where the RCPV intersects the Z-axis was too close to the origin to see that it was not coincident.

Creating a new datum and reference ellipsoid is easy. The third statement creates the new ellipsoid with values scaled from those for GRS 80, which are pre-programmed into the package. myEllipsoid's semimajor axis has the same length as that for GRS80 but its semiminor axis is only 3/4's that of GRS 80. Once the datum and ellipsoid are defined, they are available for use.

myEllipsoid = "myEllipsoid";
myDatum = "myDatum";
AddEllipsoid[myEllipsoid, semiMajor[GRS80], 
      SemiminorAxis -> 
        0.75  semiMinor[GRS80]]
AddDatum[myDatum, myEllipsoid]

Now, a function is defined that returns the graphics that draw the RCPV, it's location with respect to the XYZ coordinate system and the geodetic latitude, longitude (BLH) system. The function accepts a latitude and longitude in radians, with the name of a datum, as inputs.
1) sys holds a data structure to transform BLH and XYZ positions to and from a local horizontal (ENU) coordinate system that is defined with its origin at the point specificed by the function's arguments.
2) e is a shorthand for the reference ellipsoid associated with the datum specified by the function's argument, d.
3) p1 is a point on or directly above {ϕ, λ}. This is accomplished by creating an ENU point with coordinates {0,0,h} where h is the ellipsoid height of the point. Mathematica draws points in a Cartesian coordinate system, so p1 is transformed to XYZ using the ⇒ operator. This could also have been accomplished just as easily by creating a blh point at {ϕ,λ,h} and transforming it to XYZ.
4) p2 is the point where the RCPV intersects the Z axis. This is found simply by creating an ENU point with coordinates {0,0,-ν} and transforming it to XYZ. This works because ENU is a locally horizontal coordinate system, which implies that the UP direction is parallel with the RCPV. It is also known that ν is the distance from the reference ellipsoid to the Z-axis.
5) p3 is a point in the equatorial plane in the meridian of the point of interest.
6) x, y, and z indicate the point of interest's location in the XYZ coordinate system.
7) n is a line in the equatorial plane from the center of the reference ellipsoid to p3.
8) p is the parallel holding the point of interest
9) m is a section of the meridian holding the point of interest.

DrawPoint[ϕ_?NumericQ, λ_?NumericQ, 
        d_String] := 
Module[
{sys = 
            MakeENUSystem[ϕ °, λ °, 0, 
              d], 
e = 
            referenceEllipsoid[d], 
p1, p2, p3, rcpv, x, y, 
          z, p, m, n}, 
p1 = 
          MakeENU[0, 0, 0.0  semiMajor[e]] ⇒ {xyz, sys} // 
            ToVector; 
p2 = 
          MakeENU[0, 
                0, -ν[ϕ °, 
                    e]] ⇒ {xyz, sys} // 
            ToVector; 
p3 = 
          MakeBLH[0, λ °, d] ⇒ xyz // 
            ToVector; 
rcpv = {Thickness[0.005], 
            Dashing[{0.02, 0.0125}], RGBColor[1, 0, 0], 
            Line[{p1, p2}]}; 
x = {RGBColor[1, 0, 1], 
            Line[{{p1[[1]], 0, 
                  0}, {p1[[1]], 
                  p1[[2]], 
                  0}}]}; 
y = {RGBColor[1, 0, 1], 
            Line[{{0, p1[[2]], 
                  0}, {p1[[1]], 
                  p1[[2]], 
                  0}}]}; 
z = {RGBColor[1, 0, 1], 
            Line[{{p1[[1]], 
                  p1[[2]], 
                  0}, {p1[[1]], 
                  p1[[2]], 
                  p1[[3]]}}]}; 
n = {RGBColor[1, 0, 1], 
            Line[{{0, 0, 0}, p3}]}; 
p = {RGBColor[1, 0, 
              1], Line[
              Table[ToVector[
                  MakeBLH[ϕ °, L °, 
                      d] ⇒ xyz], {L, 0, 360, 10}]]}; 
        m = {RGBColor[ .9, 0,  .35], 
            Line[Table[
                ToVector[
                  MakeBLH[L °, λ °, 
                      d] ⇒ xyz], {L, 0, 90, 
                  10}]]}; 
{rcpv, x, y, z, p, m, 
          n}
]

Now, draw the image. First, story the length of the semimajor axis in a as a shorthand. Then create a circle for the equator and lines for the XYZ axes. rcpv holds the graphics generated by the function. shell is a set of points defining the backdrop reference ellipsoid. Note the simple use of geographic coordinates transformed to XYZ and then to a Mathematica vector for display. The Show statement displays the results.

a = semiMajor[myEllipsoid];
equator = 
        Line[Table[
            ToVector[
              MakeBLH[0 °, λ °, 
                  myDatum] ⇒ xyz], {λ, 0, 360, 
              10}]];
Off[General::spell1];
xAxis = 
        Line[{{-1 a, 0, 0}, {1.1  a, 0, 
              0}}];
yAxis = 
        Line[{{0, -1 a, 0}, {0, 1.1  a, 
              0}}];
zAxis = 
        Line[{{0, 0, -1 a}, {0, 0, 
              1.1  a}}];
On[General::spell1];
rcpv = DrawPoint[55, 135, myDatum];
shell = 
        ListSurfacePlot3D[
          Table[ToVector[
              MakeBLH[ϕ °, λ °, 
                  myDatum] ⇒ xyz], {λ, 180, 360, 
              20}, {ϕ, -90, 90, 10}], 
          DisplayFunction -> Identity];
Show[
        shell, 
Graphics3D[{
{Thickness[
                0.0025], RGBColor[0, 0, 1], xAxis, yAxis, zAxis, 
              equator}, 
rcpv
}], 
        Boxed -> False, ViewPoint -> {1.5, 1.0,  .6}, 
        DisplayFunction -> $DisplayFunction, 
        PlotLabel -> "GeometricalGeodesy"];