Projections for star charting and lunar mapping

The following is a combination of a reply I made (in August 1999) to an inquiry about projections for star charts, plus a reply made (in April 2002) to an inquiry on the Guide users mailing list about projections suitable for use in lunar mapping, with a few links for better navigation.

 Projections best suited for star charts:  
   Generally speaking,  you'll want to stick to azimuthal projections
on star charts.  There are plenty of good reasons to use other
projections on terrestrial charts (the Mercator renders lines
of constant bearing on the earth as straight lines;  various
cylindrical projections give you good accuracy on most of the earth
at the expense of trashing the poles;  and so on).  But few of
those reasons apply when you're dealing with star charts.

   Guide,  and SkyMap Pro,  and a lot of other programs,  use
a simple stereographic projection for most purposes,  because
_for star charts_, it makes a lot of sense.  You can plot a full
hemisphere and get something like the chart in the center pages
of _Sky & Tel_ or _Astronomy_:  the shapes are preserved over
small areas (it's a "conformal" projection) and the scale doesn't
change too badly as you move out from the center of the chart.

   Meteor observers like the gnomonic projection,  because if you
plot meteor paths on such a chart,  they're straight lines and they
converge at the radiant.  (The distortion on a chart covering a big
chunk of the sky is pretty ugly,  though.)

   Neither of the above can show you a full-sky,  both-hemispheres
field of view.  The equidistant azimuthal will do that.  I think
it's uglier than sin myself,  but people asked me to put it into
Guide,  and I swallowed my aesthetic horror and put it in there
(along with a few other projections I personally consider ugly,
like the orthographic... Guide _does_ use this for terrestrial
charts,  though,  and it looks just fine in that use.)

   For star charts,  if you implement the stereographic,  gnomonic,
and possibly equidistant,  you'll probably be pretty happy.

   As far as I know,  Dave Chandler's _Deep Space_ program is the
winner and champion for implementing bazillions of chart projections.
He includes,  as I recall,  a cylindrical projection in galactic
coordinates (maybe Mercator,  I couldn't really tell.)  This is not
such a bad idea.  If you're plotting data that is clustered around
the galactic equator (novae,  etc.), it can make a lot of sense.

   There are two cases in astronomy where conic projections might
be helpful:  (a) if you wanted to make a whole set of charts with
the same two standard parallels, then they would match each other
left to right;  (b) if you wanted to make a Uranometria "clone".

  Math behind some simple azimuthal projections:  Almost
all projections (and all azimuthal ones) assume that you begin with
a central point,  (ra0,  dec0),  which happens to match the center
of your chart.  In an azimuthal projection,  if you measure an
azimuth (or,  to use the celestial rather than the cartographic
term,  a position angle) from this point,  you get a correct
result. Distortions in the rest of the sky mean that this is the
_only_ point for which this holds true.

   In addition to this central point,  you've got a point (ra, dec)
that you want to project,  resulting in a "plane" point (x, y).
Here goes:

delta_ra = ra - ra0;          /* determine difference in RA */
x1 = cos( dec) * sin( delta_ra);
y1 = sin( dec) * cos( dec0) - cos( dec) * cos( delta_ra) * sin( dec0);

Two side points:
   (1) You could set x = x1,  y = y1,  do no more work,  and have an
      orthographic projection.  In fact,  you might want to experiment
      with use of the orthographic projection,  make sure you've gotten
      everything right,  and then come back and add the bits given
      below for the stereographic projection.
   (2) For a given chart,  you'll run through bazillions of (ra,  dec)
      values,  while (ra0,  dec0) stays constant.  So it can pay to
      save math by computing the sines and cosines of ra0 and dec0
      "up front",  before you plot a single star.  (This used to be
      absurdly important in the Bad Old Days before math chips became
      common.  It's not so important now,  but it's still a good idea.)

z1 = sin( dec) * sin( dec0) + cos( dec) * cos( dec0) * cos( delta_ra);

   Physically speaking,  we're treating the celestial sphere as if it
were a ball held at infinity,  and we've put Cartesian axes in it.
We're staring down at (ra0, dec0),  and the Z axis points from the
center of the ball,  up through (ra0,  dec0),  and through us.  So
a point at z1 = -1 would be on the exact far side of the ball,  at
(ra0 - 180 degrees,  -dec0);  a point at z1 = 0 would be on the
"horizon" (edge) of the ball,  as seen by us;  and a point at z1 = 1
would be on the near side of the ball,  as seen by us.

   The result is that (x1,  y1,  z1) is a unit vector.  (ra0,  dec0)
gets mapped to (0, 0, 1).  (ra0 + 180, -dec0) gets mapped to (0, 0, -1).
(ra0 + 90, 0) gets mapped to (1, 0, 0).  And so on.

   So our (x1, y1) corresponds to an orthographic projection;  we
can get any other azimuthal projection with the right scaling factor.
Here's the scaling factor d for a stereographic projection :

if( z1 < -.9)
   d = 20. * sqrt(( 1. - .81) / ( 1.00001 - z1 * z1));
   d = 2. / (z1 + 1.);

   This is going to be hard to explain... mathematically speaking,
all we ought to do is that last line,  d = 2. / (z1 + 1.).  But if
z1 < -.9,  it means that we're really close to the far side of the
earth,  which is a singularity in all azimuthal projections and
is mapped out at infinity.  There's an old science fiction story
titled "Singularities Make Me Nervous" (it's about black holes),
and it's usually a motto I live by.

   To sort of duck around the singularity,  I use the

if( z1 < -.9)
   d = 20. * sqrt(( 1. - .81) / ( 1.00001 - z1 * z1));

   to cause such points to be mapped out very far from the center,
but not at infinity.  That helps evade overflow and divide-by-zero
problems.  Depending on your prejudices,  it's an elegant engineering
solution,  or it's an Ugly Kludge.  But I digress.

   Finally,  we get:

x = x1 * d
y = y1 * d

   in which (x, y) is the projected point using the stereographic
projection,  and (as described above) (x1, y1) is the projected point
using the orthographic projection.

   Congratulations!  You've projected a point!

   If you wish to play around with a still simpler projection...
here's the simplest one of all:  an equidistant cylindrical,  a.k.a.
"plate-carre",  a.k.a. "Cassini":

x = (ra - ra0) / sin( dec0)
y = (dec - dec0)

   It's really not good for star charting.  Its terrestrial counterpart
gets a lot of use,  because it's simple and because no one lives near
the poles (as with the Mercator projection,  the poles are singularities).
But for star charting,  people look at Polaris and Sigma Octans every day.

   Anyway,  all this was supposed to be a lead-in for your question about
determining scale.  For star mapping,  the above x and y values are
dimensionless quantities.  For a terrestrial projection,  they would be
in units of the earth's radius.  You could multiply them by the radius of
the earth (6378140 meters) to get a value in meters.  And if you were
making,  say,  a 1:1000000-scale map,  you would divide your value by a
million and figure out how many meters to use on a plain paper map.

   In our case,  we do something a little different.  Let's take the
case you describe:

   "...Say if I wanted to plot an area from 0h - 4h RA and 90* -
50* Dec on a 640x480 area..."

   In this case,  you are attempting to put an area that is 40 degrees
high into a screen that's 480 pixels high.  Thus,  there are
480*(180/pi)/40,  or 12*180/pi,  pixels per radian,  in the scale you're
aiming for.  So:

screen_height_in_degrees = 40;
screen_height_in_pixels = 480;
screen_width_in_pixels = 640;
scale = screen_height_in_pixels * (180 / pi) / screen_height_in_degrees;
pixel_x = (screen_width_in_pixels / 2) +/- (x * scale);
pixel_y = (screen_height_in_pixels / 2) +/- (y * scale);

   OK,  now I have to explain all that mess... you can probably see
where the scale came from.  The screen width and height divided by two
are put there to ensure that the origin is shifted to center of the
screen. The '+ or -' is there because we've not really defined which way
the axes go.  On most computers,  y pixel values increase as you go from
top to bottom,  totally contrary to what anyone who has taken an
introductory math course would expect.  And RA increases "backwards",
too.  The best guide is:  if your chart shows up flipped left/right,
invert the sign in the 'pixel_x' equation.  If it's flipped top to
bottom,  do the same for the sign in the 'pixel_y' equation.

   You may want to be able to flip both signs at will,  anyway,  to
simulate the view achieved in assorted image-inverting optical systems.

   In the above,  we used the vertical distances to compute scale.
Using the horizontal distances,  we run into two problems.  First off,
due to aspect,  the "distance" from 0h to 4h RA isn't actually 60 degrees;
it's somewhat less than that.  Let's fudge and call it about 30 degrees,
which is about right (we _have_ to fudge to some extent on this,  since
it's not a clearly defined problem.)  We could get a decent result by
multiply the 60-degree "pure" difference in RA by the sine of the
central declination,  70 degrees.  That would give us something a little
less than 30 degrees.

   Anyway,  we run the math through and see our second problem:  we
get some scale,  call it "x_scale",  that doesn't match the scale we
got when we did it vertically,  call it "y_scale".  Run in circles,
scream and shout...

   There are two solutions to this.  Solution one:  take the lower
of the two scale values.  This way,  you'll get the whole area in
there.  You might get more range in RA than you wanted,  or more
range in dec than you wanted.  But you'll get the chunk of sky you
wanted.  Solution two:  set scale = 1,  and project all the various
points you want to be absolutely sure to have in the chart.  Maybe
you want to be sure to have the centers of the sides of the rectangle
in the chart.  Maybe you've got the RA/dec of the corners of a
constellation,  and you want to be sure the whole constellation makes
in on the chart.  Maybe you've got four corners in RA/dec that you
want to get on the chart.

   In the process of projecting everybody,  you find out that one point
is 20% off the screen area.  Everybody else is less than that.  So you
divide the scale by 1.2,  and you're all set.

   Maybe you found that the "outermost" point was only 73% of the way to
the screen edge.  Fine;  divide the scale by .73.

   The solution I actually used in Guide was this,  far simpler,  one.  If
you tell Guide you want a 20-degree field of view,  and the screen is
640x480,  it looks at the shortest (480-pixel) axis.  It sets the scale
to be 480 pixels per 20 degree,  or 24 pixels per degree.  End of story.

   For really wide fields of view,  you _do_ still run into a problem,
even with this simpler scheme.  The scale may be 24 pixels/degree at
the center of the screen,  but it might rise to 48 pixels/degree at the
edge of the field of view.  My advice:  be aware of this,  then ignore it.

Projections for lunar cartography (selenography?):    There are
a couple of projections somewhat more suited for the unique conditions
of lunar cartography.  Those "conditions" are that we usually care most
about one hemisphere,  i.e.,  the one facing us,  with perhaps some
extra charting of the libration zones.

   Antonín Rükl's  Atlas of the Moon tackled the display
of the near side of the moon by showing it with an orthographic
projection. This matches exactly the way the moon would look if
seen at zero libration from infinity,  and is quite easy to work
out;  if you look above,  you'll see that the orthographic is the
simplest of the azimuthal projections.  Just set the "ra0, dec0"
values to correspond to libration terms,  and you can get a version
that is "customized" to the current lunar orientation.

   In this projection,  round craters will become ellipses,  with
the minor axis pointing through the center of the chart.  There won't
be any change in the major axis (scale is preserved in that direction),
but the minor axis should be equal to (major axis * z1).  Near the
center of the near side,  z1 is close to 1 and there is nearly no
distortion.  Near the limbs,  z1 is closer to 0 and the "circular"
craters become extremely flattened.

   Rükl provides eight charts of the libration area.  Four use
a polar stereographic (i.e.,  dec0 = +90 or -90) projection,  and
four use a Mercator projection.  In a perfect world,  he probably
would have wanted to use a transverse Mercator projection,  but
that would have been a bit of a challenge in the pre-computer era.
Still,  if you use the above formulae for a stereographic projection,
with (ra0, dec0) moved to the libration area,  you can get a decent
chart of a chunk of the libration zone.

   The scaling of craters will have to change a bit.  One pleasant
aspect of the stereographic projection is that craters stay round;
they just get scaled up by a factor of 2 / (1 + z1).  Thus,  craters
at the center of the projection don't get scaled at all;  those near
90 degrees from the center get scaled up 2:1.  (If you look at the
"all-sky" charts in the center of each of _Sky & Tel_ or _Astronomy_,
you can see this effect.  Constellations near the horizon are a bit
larger than they ought to be,  compared to those near the zenith.
This may be most apparent if you get two issues three months apart,
so you can find two constellations that go from "one near the zenith,
one rising" to "one setting,  the other near the zenith".)

   There's also a possible compromise projection:  one which looks
superficially much like the orthographic projection,  but in which
the "edge" of the projection is a bit more than 90 degrees from the
center of the near side,  allowing a bit more than half of the moon
to be shown.  At some point,  I may come back and document this beast.
It requires a little more math,  but nothing too unpleasant.