# Understanding geodesic envelope?

I've been looking at some various ways to calculate geodesics. One particularly good site collecting load of helpful information, is this one. There you can also find a link to a map that shows a circle whose geodesic radius is accurately calculated by the WGS84 ellipsoidal approximation. In the resulting map, there are some additional info presented (as red arcs), called `envelope`. I'm struggling to understand what these are and what exactly they represent. The only partially useful resource I can find is from this website, but looking at the faint lines drawn on 2D image of a stone, does not enlighten anyone, especially as you cannot rotate the images, to see where those lines goes, and intersects. OR from another calculation: From the website, the description only say:

The geodesic envelopes as red curves; all the geodesics emanating from lat1, lon1 are tangent to the envelopes (providing they are extended far enough). The number of solutions to the inverse problem changes depending on whether lat2, lon2 lies inside the envelopes. For example, there are four (resp. two) approximately hemispheroidal geodesics if this point lies inside (resp. outside) the inner envelope (only one of which is a shortest path).

Q: How can I better understand what those red arcs represent?

1. Why are they useful on this picture?
2. How are they generated?

I'm more of a visual person, seeking a visual explanation, rather than a purely verbal description.

It looks like you've figured things out OK. In case it helps, I give here the MATLAB/Octave code that I used to generate this Wikipedia figure showing the envelope: The envelope is formed by the light blue lines which are geodesics emanating from a single point on the other side of the ellipsoid (the flattening of the ellipsoid is 1/10). If this were a sphere, all the geodesics would go through the antipodal point and the envelope would collapse to a point. The green curves are geodesic circles and the red curve is the cut locus which is a portion of the circle of latitude through the antipodal point.

This just shows the first envelope obtained by letting the geodesic go slightly more than 1/2 way round the ellipsoid. Bigger envelopes are generated if the geodesics go slight more than 1 time around (envelope is then centered on the originating point), or 1-1/2 times around (centered on the antipodal point), etc.

This code presumes that you have installed MATLAB Central packge 50605. I'm sorry that it's not very pretty; it was just intended to "get the job done".

``````function geod_envelope
thick=1;
thin=0.3;
red=[179,27,27]/255;
white=[1,1,1];
black=[0,0,0];
blue=[0,19,56]/100;
green=[9,45,27]/100;
xsize=4*0.66;ysize=4*0.66;
name='geod-envelope';

figure(1);

a = 1;
r = 10;
ell = [a, flat2ecc(1/r)];
lat1=-30;
lon1=0;
lat0=-lat1;lon0=180;
hold off;
for azi1=[-180:3:180],
[lat,lon] = geodreckon(lat1,lon1,[120:2:240]',azi1,ell,1);
xyz=proj([lat,lon],lat0,lon0,a,r);
xyf=frontx(xyz);
plot(xyf(:,1),xyf(:,2),'-',...
'LineWidth',thin,'Color',0.3*blue+0.7*white);
if azi1 == -180,
hold on;
end
end

[latc,lonc,~, ~, ~, ~, ~, s12c] = ...
geodreckon(lat1,lon1,180,90,ell,1);
loncut=180+[-8:8]'/8*(180-lonc);
latcut=latc+0*loncut;
xyz=proj([latcut,loncut],lat0,lon0,a,r);
xyf=frontx(xyz);
plot(xyf(:,1),xyf(:,2),'-','LineWidth',thick,'Color',red);
s12d = geoddistance(lat1,lon1,-lat1,lon1+180,ell);
ds=(s12d-s12c)/2;
for s=[0:6],
s12=s12c+ds*s;
[lat,lon] = geodreckon(lat1,lon1,s12,[-180:2:180]',ell);
xyz=proj([lat,lon],lat0,lon0,a,r);
xyf=frontx(xyz);
plot(xyf(:,1),xyf(:,2),'-','Color',green);
end

hold off;

axis image;
axis(0.25*[-1,1,-1.109,0.891]);
axis off;
end

function xyz = proj(ll, lat0, lon0, a, f)
%proj convert ll to ortho from pov of lat0 lon0
if f>1
f=1/f;
end
e2=f*(2-f);
e2m=(1-f)^2;
degree=pi/180;
phi=[ll(:,1);lat0]*degree;
lam=[ll(:,2);lon0]*degree;
sphi=sin(phi);
cphi=cos(phi);
n = a ./ sqrt(1 - e2 * sphi.^2);
slam = sin(lam);
clam = cos(lam);
z = e2m * n .* sphi;
x = n .* cphi;
y = x .* slam;
x = x .* clam;
east0 = [-slam(end)          ;  clam(end)          ; 0        ];
north0 =[-clam(end)*sphi(end); -slam(end)*sphi(end); cphi(end)];
up  =   [ clam    .*cphi     ,  slam    .*cphi     , sphi     ];
up0 =   up(end,:)';

pts = [x(1:end-1)-x(end), y(1:end-1)-y(end), z(1:end-1)-z(end)];
x = pts * east0;
y = pts * north0;
z = up(1:end-1,:) * up0;
xyz = [x, y, z];
end

function xy = frontx(xyz),
c=xyz(:,3)>=-0.01;
c=1+0./c;
xy = [xyz(:,1).*c, xyz(:,2).*c];
end
``````
• Many thanks to you! Not only for this answer but also to all the work you have done and presented in the field. For someone who's not affiliated to the GIS field, it really take a lot of effort to wrap your head around all the details. Your code and examples have been really helpful in doing this. Feb 4 '18 at 18:05

Since posting my question, I have made some progress in better understanding the use of the `envelope`. To illustrate my findings, let's first consider what happens to the concept of a geodesic radius in the two cases of (1) a Sphere, and (2) an 1-Ellipsoid. (I call it a "1-ellipsoid" as it is a spheroid only distorted along one axis and not both. For distortions along both axes, I would call it a 2-Ellipsoid.)

### Working Example

As a starting point (A) we chose Wellington (New Zealand) and try to find it's approximate geographic opposite (B) on the planet, we will end up slightly outside Salamanca (Spain).

``````(A) Wellington:           -41.3286847, 174.8068018
(B) Salamanca (Alaejos):   41.3046662, -5.1947402
``````

Here we will use the direct problem to find an approximate location for B. The inverse problem would find the exact position.

Spherical geodesics

But let's take a step back. If we first assume earth is a perfect sphere, then we can select any point (A) on that surface, and draw a geodesic to some radial distance, R, to make a big circle (with radius R) on that sphere. As we increase the radius to more than 1/4 way around the sphere, the circumference of the the circle will become larger than the "circumference of the sphere" which is its diameter. Now, if we start looking at the sphere from the other side, that of the diametrically opposite of point A, which we call point B. We see that radius (approaching us) is still a circle, but now shrinking. (Imagine a rubber band being stretched around a ball to the other side.)

There are two facts to observe in this case:

• Drawing any geodesic of any length >PiR will always pass through the opposite "pole".
• The geodesic is a closed circle.

Ellipsoidal geodesics

However, Earth is actually an ellipsoid, which in its simplest form, is a sphere that has been flattened on its poles. Therefore, a geodesic radius used to define a "circle" will almost never actually look like a circle, unless its center is at one of the poles. That is because the geodesic distance changes depending on which direction (bearing) you draw the radius, except at the polar symmetry axis. The two things to note are:

• Geodesics are not closed paths
• The mathematical circles generated by the geodesic radius are not Cartesian circles

Looking at the following picture series of geodesic radii, we find that the radial "circle" changes shape as it expands around the ellipsoid. Of course mathematically it is always a perfect circle, but by human perception, it is distorted from our flat world experience.

As we progress through the picture series we see that, as the radius (in green) is getting closer to point B, the circle is getting more and more distorted. At the point when it first reaches the inner `envelope` (in red), it develops a cusp. This cusp changes shape until the radius exits the other side of inner envelope. At that point, it once again resembles the shape of a distorted circle or ellipse. Mathematically, I suppose this is equivalent to having inverted the circle.

In the last picture, we have increased the radius to loop around the planet, one more time, to end up near B. However, now see the meaning of the envelope. Because ellipsoidal geodesics are not closed loops, the envelope astroids represent the geodesic difference between each non-closed loop.

I'm sure someone else can write and make a much better explanation, so please feel free to comment and revise.