Celestial Cartographic Library

Norman Vine (always a source of fun ideas) and I were speculating a few weeks ago about using MapServer to produce maps of the sky. I finally made the time to demonstrate that it is as straight-forward as we had thought. I'm using PCL with a MapServer-based renderer, but everything should be reproducible with more traditional methods. The star data was obtained from Steve Dick's Corrected Bright Star List, 1997 Astronomical Almanac H2-H31: http://aa.usno.navy.mil/publications/AstroAlmanac/bstar97.html.

To begin, I wrangled the list into an OpenOffice spreadsheet and saved it as a 15-field CSV-formatted file. The following Python script processed the spreedsheet into a PCL feature store. Believe it or not, we can use a lat/long coordinate system for objects in the sky. The right ascension is multiplied by 15 because it is in units of hh:mm:ss, and we take the opposite sign of this value because right ascension increases in the opposite sense of our longitude:

import csv
from cartography import data, spatial

properties = (['FB_DES', 'string'],        # Flamsteed/Bayer Designation
              ['MAGNITUDE', 'float'],      # Visual Magnitude
              ['SPECTRAL_TYPE', 'string']) # such as Sol's G2 V

star_srs = spatial.SpatialReference(['+proj=latlong'])
star_store = data.MemoryFeatureStore(properties)

fd = open('bstar97.csv', 'r')
reader = csv.reader(fd)
for row in reader:
    try:
        # values of our properties
        des = ' '.join([row[1].strip() or row[0].strip(), row[2].strip()])
        mag = float(row[11])
        spe = row[14].strip()

        # the location of the star as a Point
        # right ascension (longitude)
        x = -15.0*(float(row[4]) + float(row[5])/60 + float(row[6])/3600)
        x -= 180.0

        # declination (latitude)
        y = float(row[7]) + float(row[8])/60 + float(row[9])/3600
        p = spatial.Point(x, y)

        # add to the data store
        star_store.add(p, [des, mag, spe])
    except:
        pass

fd.close()

A graticule is as good a visual aid in a star map as it is in a map of the Earth. Set up a graticule (from http://sgillies.net/blog/16) and styling rules for its presentation:

from cartography import mapping, query, rendering

grat_store = data.DiskFeatureStore('../grat/graticule.shp')
grat_style = mapping.Style()
line_symb = mapping.LineSymbolizer({'stroke': '#CCCCCC', 'opacity': 0.5})
grat_style.rules.append(mapping.Rule('catch-all', line_symb))

grat_layer = mapping.Layer(grat_store, 'LINE', star_srs)
graticule = (grat_layer, grat_style)

And now the stars themselves. First thing I do is define a sequence that associates filter expressions with marker sizes. The styling rules for different stellar brightness levels will be generated from this 'mags' sequence:

mags = \
[
['f.MAGNITUDE <= 0',       15],
['0 < f.MAGNITUDE <= 1',   11],
['1 < f.MAGNITUDE <= 2',    7],
['2 < f.MAGNITUDE <= 3',    5],
['3 < f.MAGNITUDE <= 4',    3],
['4 < f.MAGNITUDE <= 5',    1]
]

circle = mapping.Mark('circle')
star_style = mapping.ExtrinsicStyle()

for m in mags:
    symb = mapping.PointSymbolizer(mark=circle, size=m[1])
    rule = mapping.Rule(m[0], symb, query.PythonFilter(m[0]))
    star_style.rules.append(rule)

# Text symbolizer to label brightest stars
font = mapping.Font(family='Vera', size=8.0)
text_symb = mapping.TextSymbolizer(label='FB_DES', font=font,
                                       placement='AUTO')
star_style.rules.append(mapping.Rule('m0t', text_symb,
                  query.PythonFilter('f.MAGNITUDE <= 1')))

# Wrap the star data store up in a Layer and associate with the style
star_layer = mapping.Layer(star_store, 'POINT', star_srs)
stars = (star_layer, star_style)

Nothing fancy, just rendering the brighter stars as bigger filled circles. Now we can render a map of the stars of the Northern Hemisphere using a polar Lambert Azimuthal Equal Area projection:

FONTSET = '/home/sean/projects/ms_45/mapserver/tests/fonts.txt'
mapper = rendering.MapRenderer('MAPSERVER', fontset=FONTSET)

srs = spatial.SpatialReference(['+proj=laea', '+lat_0=90.0'])
ul = spatial.Point(-10000000, 10000000)
lr = spatial.Point(10000000, -10000000)
view = mapping.View(srs, ul, lr)
size = [600, 600]
image = mapper.render([graticule, stars], view, 'image/png', size)
f = open('stars_polar.png', 'wb')
f.write(image.getdata())
f.close()

stars_polar.png:

Polaris is right near the center, Cassiopeia above. Left of center is Vega (alpha Lyra), featured in Carl Sagan's novel Contact. Alpha Centauri appears in the lower left corner. Orion is on the right edge, greatly distorted. LAEA is not an ideal projection for a sky map. I also tried Polar Stereographic and Azimuthal Equidistant projections but ran into some issues with projecting the graticule.

Next we zoom in on Orion, my favorite constellation, using an LCC projection:

srs = spatial.SpatialReference(['+proj=lcc', '+lat_1=10.0', '+lon_0=100.0'])
ul = spatial.Point(-2251981, 2350184)
lr = spatial.Point(1611262, -1237944)
view = mapping.View(srs, ul, lr)
size = [400, 400]
image = mapper.render([graticule, stars], view, 'image/png', size)
f = open('stars_orion.png', 'wb')
f.write(image.getdata())
f.close()

The result, stars_orion.png:

From left to right: Betelgeuse (alpha Ori) the 10th brightess star in the sky, the Rigel (beta Ori) of many Star Trek references, and Aldebaran (alpha Tau).

Update: Frank Warmerdam educated me on PROJ4 and its use by MapServer during rendering. His suggestion was to "inline" the graticule as I had the stars. Simply change the grat_style to an instance of mapping.ExtrinsicStyle. It's costly, as any extrinsic style, but works. Here is the Polar Stereographic map that I wanted:

chart_srs = spatial.SpatialReference(['+proj=stere','+lat_0=90.0'])
ul = spatial.Point(-14000000, 14000000)
lr = spatial.Point( 14000000, -14000000)
view = mapping.View(chart_srs, ul, lr)
size = [600, 600]
image = mapper.render([graticule, stars], view, 'image/png', size)
f = open('stars_stere.png', 'wb')
f.write(image.getdata())
f.close()