iconEuler Examples

Geometry on the Earth

by R. Grothmann

In this notebook, we want to do some spherical computations. The
functions are contained in the file "spherical.e" in the examples
folder. 

See: Reference

We need to load that file first.
>load "spherical.e";
To enter a geographical position, we use a vector with two coordinates
in radians (north and east, negative values for south and west). The
following are the coordinates for the Campus of the Catholic
University in Eichstätt.
>EI=[rad(48,53.371),rad(11,11.330)]
[ 0.853283035542  0.195281981124 ]
You can print this position with sposprint (spherical position
print).
>sposprint(EI)
N 48°53.371' E 11°11.330'
Let us add two more towns, Ingolstadt and Neuburg.
>IN=[rad(48,45.854),rad(11,25.055)]; ND=[rad(48,44.248),rad(11,10.712)];
>sposprint(IN), sposprint(ND),
N 48°45.854' E 11°25.055'
N 48°44.248' E 11°10.712'
First we compute the vector from one to the other on an ideal
ball. This vector is [heading,distance] in radians. To compute
the distance on the earth, we multiply with the earth radius at
a latitude of 48°.
>br=svector(EI,IN); deg(br[1]), br[2]*rearth(48°)->km
129.671788237
21.766496146
This is a good approximation. The following routines use even better
approximations. On such a short distance the result is almost the
same.
>esdist(EI,IN)->" km"
21.7654976508 km
There is a function for the heading, taking the elliptical shape of
the earth into account. Again, we print in an advanced way.
>sdegprint(esdir(EI,IN))
    129.67°
The angle of a triangle exceeds 180° on the sphere.
>asum=sangle(IN,EI,ND)+sangle(EI,IN,ND)+sangle(IN,ND,EI); deg(asum)
180.000207484
This can be used to compute the area of the triangle. Note: For
small triangles, this is not accurate due to the subtraction error
in asum-pi.
>(asum-pi)*rearth(48°)^2->" km^2"
146.771319034 km^2
There is a function for this, which uses the mean latitude of the
triangle to compute the earth radius, and takes care of rounding
errors for very small triangles.
>esarea(IN,EI,ND)->" km^2"
146.757696278 km^2
We can also add vectors to positions. A vector contains the heading
and the distance, both in radians. To get a vector, we use svector. To
add a vector to a position, we use saddvector.
>v=svector(EI,IN); sposprint(saddvector(EI,v)), sposprint(IN),
N 48°45.854' E 11°25.055'
N 48°45.854' E 11°25.055'
These functions assume an ideal ball. The same on the earth.
>sposprint(esadd(EI,esdir(EI,IN),esdist(EI,IN))), sposprint(IN),
N 48°45.854' E 11°25.055'
N 48°45.854' E 11°25.055'
Let us turn to a larger example, which I found in the Internet,
Brandenburger Tor in Berlin, and the Tejo bridge in Lissabon. For
Lissabon, take a negative value for the west coordinate.
>Berlin=[52.5164°,13.3777°]; Lissabon=[38.692668°,-9.177944°];
>sposprint(Berlin), sposprint(Lissabon)
N 52°30.984' E 13°22.662'
N 38°41.560' W 9°10.677'
According to Google Earth, the distance is 2315.75km. We get a
good approximation.
>esdist(Lissabon,Berlin)->" km"
2313.67485976 km
The heading is the same as the one computed in Google Earth.
>degprint(esdir(Lissabon,Berlin))
41°3'5.08''
However, we do no longer get the exact target position, if we
add the heading and distance to the orginal position. This is
so, since we do not compute the inverse function exactly, but
take an approximation of the earth radius along the path.
>sposprint(esadd(Berlin,esdir(Berlin,Lissabon),esdist(Berlin,Lissabon)))
N 38°41.557' W 9°10.680'
The error is not large, however.
>sposprint(Lissabon),
N 38°41.560' W 9°10.677'
Of course, we cannot sail with the same heading from one destination
to another, if we want to take the shortest path. Imagine, you fly NE
starting at any point on the earth. Then you will spiral to the north
pole. Great circles do not follow a constant heading!

The following computation shows that we are way off the correct
destination, if we use the same heading during our travel.
>dist=esdist(Berlin,Lissabon); hd=esdir(Berlin,Lissabon);
Now we add 10 times one-tenth of the distance, using the heading
to Lissabon, we got in Berlin.
>p=Berlin; loop 1 to 10; p=esadd(p,hd,dist/10); end;
The result is far off.
>sposprint(p), skmprint(esdist(p,Lissabon))
N 41°0.566' W 12°5.348'
   357.810km
As another example, let us take two points on the earth at the same
lattitude.
>P1=[30°,10°]; P2=[30°,50°];
The shortest path from P1 to P2 is not the circle of lattitude 30°,
but a shorter path starting 10° further north at P1.
>sdegprint(esdir(P1,P2))
     79.69°
But, if we follow this compass reading, we will spiral to the north
pole! So we must adjust our heading along the way. For rough purposes,
we adjust it at 1/10 of the total distance.
>p=P1;  dist=esdist(P1,P2); ...
   loop 1 to 10; dir=esdir(p,P2); sdegprint(dir), p=esadd(p,dir,dist/10); end;
     79.69°
     81.67°
     83.71°
     85.78°
     87.89°
     90.00°
     92.12°
     94.22°
     96.29°
     98.33°
The distances are not right, since we will add a bit off error, if we
follow the same heading for too long.
>skmprint(esdist(p,P2))
     0.203km
We get a good approximation, if we adjust out heading after each 1/100
of the total distance from Berlin to Lissabon.
>p=Berlin; dist=esdist(Berlin,Lissabon); ...
   loop 1 to 100; p=esadd(p,esdir(p,Lissabon),dist/100); end;
>skmprint(esdist(p,Lissabon))
     0.057km
For navigational purposes, we can get a sequence of GPS position along
the great circle to Lissabon with the function navigate.
>load spherical; v=navigate(Berlin,Lissabon,10); ...
   loop 1 to rows(v); sposprint(v[#]), end;
N 52°30.984' E 13°22.662'
N 51°21.596' E 10°34.058'
N 51°13.643' E 10°15.858'
N 51°3.637' E 9°53.255'
N 50°50.666' E 9°24.433'
N 50°33.189' E 8°46.420'
N 50°8.382' E 7°53.991'
N 49°30.450' E 6°37.042'
N 48°25.413' E 4°33.183'
N 46°9.155' E 0°40.986'
N 38°41.557' W 9°10.680'
We write a function, which plots the earth, the two positions, and the
positions in between.
>function testplot ...
useglobal;
plotearth;
plotpos(Berlin,"Berlin"); plotpos(Lissabon,"Lissabon");
plotposline(v);
endfunction
Now plot everything.
>plot3d("testplot",>own,>user,zoom=4):

Geometry on the Earth

Or use plot3d to get an anglyph view of it. This looks really great
with red/cyan glasses.
>plot3d("testplot",own=1,anaglyph=1,zoom=4):

Geometry on the Earth

Examples Homepage