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)); else 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.