Earth as a sphere

This page shows several examples using routines from the JHU/APL/S1R IDL Library (called the S1R library below) that deal with the earth as a sphere.

Several examples are given below


Return to the Spheres main page

Plotting coastlines

A file containing world coastline data is available in the S1R IDL library. This example shows how to access and use this data.

The file is known as a res file (for results file). This format was developed to contain results of data processing in a simple self describing unformatted file. Such files are handled using the res*.pro routines in the S1R IDL library.

The world coastline data file is in XDR format and contains 3 arrays:

  Latitudes in degrees (south < 0)
  Longitudes in degrees (west < 0)
  Pencode
The coordinate arrays contain multiple curves. The pencode is 1 except at the start of each curve, where it is 0. This allows the individual curves to be accessed.

The world coastline data may be read as follows:

resopen,'world_xdr.res',/xdr
resget,'long',lng
resget,'lat',lat
resget,'pen',pen
resclose
The path to the file must be correctly given. It is suggested that you define a system variable, IDL_IDLUSR, that contains the name of the directory containing the S1R IDL library. Then, in UNIX for example, the file may be opened as
resopen,'$IDL_IDLUSR/world_xdr.res',/xdr.

An example plot using the world coastline data:

window, xs=200, ys=200
!p.position = [0,0,1,1]
set_isoxy,-1.1,1.1,-1.1,1.1  
erase
sphinit, rad=1, lat=40, long=-80, pa=30
sphplot, lng, lat, 1, pen
sphinit, rad=1

Other sph* routines may be used along with the coastline data.



Sunclock

The S1R IDL library routine sunclock will show maps of the earth with daylight, night, and twilights indicated. Do sunclock,/help for more details. The built-in IDL map software is used in sunclock. Here is an example plot (reduced in size).

The point on the earth directly beneath the sun may be found as follows:

sun, y, m, d, ut, app_ra=ra, app_dec=dec
st = lmst(jd,ut/24.,0)*24
lat = dec
lng = 15.0*(ra-st)

where y,m,d = year, month, and day numbers,
      ut = universal time in hours,
      ra, dec are the sun's returned RA and Dec,
      jd = Julian Day Number (jd=ymd2jd(y,m,d)),
      st = returned local sidereal time,
      lat, lng = resulting subsolar point.
See the time page for more details.

CRRES example

The 1991 CRRES satellite was a very interesting project. It released chemicals in space that created color patches of light in the sky. A post to the news group sci.space described this project and gave the coordinates of a number of chemical releases.

Some releases were at fairly low altitude and some at high altitude. A figure can be made to show the geometry of both a low and high altitude release. Two potential release times from the net post are:

 DATE    TIME      LAT.     LONG.   ALT.
                            (WEST)  (KM)
JAN 18   0503      7.4      61.0    33565
JAN 25   0843     16.6 S    38.9     5799
If a figure is to show the altitude extremes it must have the earth oriented to put the subsatellite point for both cases on the visible limb of the earth. This is done in IDL as follows.

Enter subsatellite points

lat1 = 7.4      ; High release.
lng1 = -61.0
lat2 = -16.6    ; Low release.
lng2 = -38.9
Use ll2rb to find bearing from point 1 to point 2

The built-in help for ll2rb is (ll2rb,/help):

 From latitude, longitude compute range, bearing.
 ll2rb, lng0, lat0, lng1, lat1, dist, azi
   lng0, lat0 = long, lat of reference point (deg).    in
   lng1, lat1 = long, lat of point of interest (deg).  in
   dist = range to point point of interest (radians).  out
   azi = azimuth to point of interest (degrees).       out
 Notes: A unit sphere is assumed, thus dist is in radians
   so to get actual distance multiply dist by radius.
   Useful constants:
   Radius of Earth (mean) = 6371.23 km = 3958.899 miles.

ll2rb,lng1,lat1,lng2,lat2,dist,azi
dist = dist*6371.23   ; Convert to km (dist not used below).

help, dist, azi
FLOAT           FLOAT     =       3610.04
AZI             FLOAT     =       137.803
Use rb2ll to find earth orientation

If the plot of the earth is oriented to center the point which is 90 degrees away along the surface from both subsatellite points then both points will appear on the earth's limb. There are two points that may be centered to do this. The most northern one will be selected.

The built-in help for rb2ll is (rb2ll,/help):

 From range, bearing compute latitude, longitude .
 rb2ll, lng0, lat0, dist, azi, lng1, lat1
   lng0, lat0 = long, lat of starting point (deg).     in
   dist = range to point of interest in RADIANS.       in
   azi = azimuth to point of interest (degrees).       in
   lng1, lat1 = long, lat of point of interest (deg).  out
 Notes: A unit sphere is assumed, thus dist is in radians.
   Useful constants:
   Radius of Earth (mean) = 6371.23 km = 3958.899 miles.
   Distance to horizon from height H above surface:
     For small H: dist = sqrt(2*H/R) in Radians
     For large H: dist = acos(R/(R+H)) in Radians
   To plot horizon from lat0, lng0, H:
     rb2ll,lng0,lat0,dist,makex(0,360,1),plng,plat
     plots,plng-360.,plat,psym=3

The azimuth of the desired centered point from subsatellite point 1 is: azi-90.
The range of the desired centered point from subsatellite point 1 is: !pi/2 = 90 degrees.
Note that the same central point may be obtained using the second subsatellite point as the reference point.

rb2ll, lng1, lat1, !pi/2., azi-90., lng0, lat0

help, lng0, lat0
LNG0            FLOAT     =       35.6604
LAT0            FLOAT     =       41.7657
Plot the earth

Use the method described above to read in world coastline data. This data consists of 3 arrays:

   lng = longitude
   lat = latitude
   pen = pencode
The following values are also available:
   re = earth radius = 6371.23 km (mean)
   r1 = Satellite max radius = 33565 + re
   r2 = Satellite min radius = 5799 + re
The coordinates are set up as follows:
!p.position=[0,0,1,1]
set_isoxy,-r1*1.1, re*1.5, -re*1.5, re*1.5
The figure is plotted approximately as follows:
sphinit, rad=re, lat=lat0, long=lng0   ; Orient earth.
sphlat, 0, re                          ; Draw equator.
sphplot, lng, lat, re, pen             ; Plot coastlines.
sphinit, rad=re                        ; Redo outline.
sphrad, lng1, lat1, 0, r1, maxrad=re   ; Max altitude.
sphplot, lng1, lat1, r1, psym=2        ; Mark with a symbol.
sphrad, lng2, lat2, 0, r2, maxrad=re   ; Min altitude.
sphplot, lng2, lat2, r2, psym=2        ; Mark with a symbol.
Colors were used to make the figure below.

Another view is available.