The IAU has a list of constellations and their bounds. I thought it was fascinating that Alan MacRobert's constellation patterns are partially inspired by H.A. Rey's asterisms. This first code snippet concatenates all of the boundaries listed on that site. There are two ser
files for Serpens Caput and Serpens Cauda. The code after this parses the constellation bound data. It is slightly tricky for bounds that wrap about the 24h-0hr boundary or at the pole. The algorithms break up the the wrapping constellations and linearly interpolate between these discontinuities and add the necessary points. You can plot single constellation bounds or everything, and there is functionality for grouping and coloring the constellations.
%pylab inline
output = '../data/1930_boundaries_iau.dat'
def pretty_hours(h, pos=None):
if h % 1.0 == 0.0:
return '{:.0f}h'.format(h)
else:
return '{:.2g}h'.format(h)
def pretty_degrees(d, pos=None):
return u'{}°'.format(d)
from matplotlib.ticker import FuncFormatter
hours_formatter = FuncFormatter(pretty_hours)
degrees_formatter = FuncFormatter(pretty_degrees)
Populating the interactive namespace from numpy and matplotlib
import urllib2
fn_base = 'http://www.iau.org/static/public/constellations/txt/'
abbrevs = 'And , Ant, Aps, Aqr, Aql, Ara, Ari, Aur, Boo, Cae, Cam, Cnc, CVn, CMa, CMi, Cap, Car, Cas, Cen, Cep, Cet, Cha, Cir, Col, Com, CrA, CrB, Crv, Crt, Cru, Cyg, Del, Dor, Dra, Equ, Eri, For, Gem, Gru, Her, Hor, Hya, Hyi, Ind, Lac, Leo, LMi, Lep, Lib, Lup, Lyn, Lyr, Men, Mic, Mon, Mus, Nor, Oct, Oph, Ori, Pav, Peg, Per, Phe, Pic, Psc, PsA, Pup, Pyx, Ret, Sge, Sgr, Sco, Scl, Sct, Ser1, Ser2, Sex, Tau, Tel, Tri, TrA, Tuc, UMa, UMi, Vel, Vir, Vol, Vul'
abbrevs = abbrevs.split(',')
f = open(output,'wb')
for i in range(len(abbrevs)):
abbrevs[i] = abbrevs[i].strip().lower()
textfilename = fn_base+abbrevs[i]+'.txt'
try:
response = urllib2.urlopen(textfilename)
if(response.getcode()!=200):
print('Didnt get bounds for {} at {}'.format(abbrevs[i],textfilename))
else:
data = response.read()
data = data.strip()
f.write(data)
if(data[-1]!='\n'):
f.write('\n')
except urllib2.URLError as e:
print('URLError ({0}): {1}\n{2}'.format(e.errno,e.strerror,textfilename))
except:
print('I wasnt expecting that error.')
raise
f.close()
found_consts = []
all_ras = []
all_decs = []
with open(output) as f:
for line in f:
l = line.split('|')
if(len(l)>1):
a_const = l[2].strip().lower()
found_consts.append(a_const)
(h,m,s)= l[0].split()
this_ra = float(h)+float(m)/60.+float(s)/3600.
all_ras.append(this_ra)
all_decs.append(float(l[1]))
found_consts_min = sorted(list(set(found_consts)))
print('confirmed {} constellations found'.format(len(found_consts_min)))
confirmed 89 constellations found
all_pts = zip(all_ras,all_decs,found_consts)
all_consts = {}
for pt in all_pts:
try:
all_consts[pt[2]].append((pt[0],pt[1]))
except:
all_consts[pt[2]]=[(pt[0],pt[1])]
The following algorithm is my first pass at dealing with constellations with wrapped coordinates. I could simply hard-code the constellations, but there's no fun in that. I'm on the lookout for the smart and probably already solved-a-hundred-times-over solution that cartographers have been using, but don't yet know the jargon to find it (yet).
My problem with this solution is with the very "sparse" asterisms that for argument's sake could be used. By sparse I mean covering a large swath of space, but having few descriptive nodes for the border where these may not be in the first and last quadrants. This technique won't work for them, because I am using a percent of the left-most and right-most bounds that the constellation coordinates must mutually share.
In the following example we see that Ursa Minor is actually one such edge case and can't be detected this way.
bounds = [0.0,24.0,-90.0,90.0]
def get_wrapping_ra(p=.25):
wrapping_consts = []
plot_center = ((bounds[0]+bounds[1])/2.,((bounds[2]+bounds[3])/2.))
lsb_x = (bounds[0],(bounds[0]+bounds[1])*p) # lower segment bounds ra
usb_x = ((bounds[0]+bounds[1])*(1.-p),bounds[1]) # upper segment bounds ra
for const_name in all_consts.keys():
min_ra = max_ra = plot_center[0]
for pt in all_consts[const_name]:
min_ra = min(min_ra,pt[0])
max_ra = max(max_ra,pt[0])
if(min_ra<lsb_x[1] and max_ra > usb_x[0]):
print('{} is a wrapping constellation'.format(const_name))
wrapping_consts.append(const_name)
return wrapping_consts
print('25 percent: lower and upper quadrants')
wcl = wcl0 = get_wrapping_ra(.25)
print('\n33 percent: thirds')
wcl1 = get_wrapping_ra(.33333)
print('\n40 percent: finally get Ursa Minor,\nbut get some non-wrapping constellations too..')
wcl2 = get_wrapping_ra(.4)
25 percent: lower and upper quadrants scl is a wrapping constellation oct is a wrapping constellation psc is a wrapping constellation cet is a wrapping constellation cep is a wrapping constellation peg is a wrapping constellation tuc is a wrapping constellation cas is a wrapping constellation and is a wrapping constellation phe is a wrapping constellation 33 percent: thirds scl is a wrapping constellation oct is a wrapping constellation psc is a wrapping constellation cet is a wrapping constellation cep is a wrapping constellation peg is a wrapping constellation tuc is a wrapping constellation cas is a wrapping constellation and is a wrapping constellation phe is a wrapping constellation 40 percent: finally get Ursa Minor, but get some non-wrapping constellations too.. scl is a wrapping constellation oct is a wrapping constellation psc is a wrapping constellation cet is a wrapping constellation cep is a wrapping constellation hya is a wrapping constellation uma is a wrapping constellation peg is a wrapping constellation umi is a wrapping constellation tuc is a wrapping constellation dra is a wrapping constellation cas is a wrapping constellation cam is a wrapping constellation and is a wrapping constellation phe is a wrapping constellation
Simply detecting wrapping about the right ascension bounds isn't enough to get every constellation. We got lucky with Octans, but not with Ursa Minor when looking at the upper and lower right ascension quadrant sections. This is actually a good thing, because we can make a fullproof methodology to begin with instead of waiting for such a case later.
Let's also detect constellations that are at the poles. For the moment, I'm going to cave and hardcode Octans and Ursa Minor into this because I want to see the pretty end-result but I think a better solution in the future would be to convert the coordinates to spherical and see if the poles are within the area of these shapes.
def get_pole_consts():
pole_consts = ['oct','umi']
return pole_consts
pcl = get_pole_consts()
# remove pole constellations from wrapping
# constellations list so dont duplicate items in the groups
for c in pcl:
try:
wcl.remove(c)
except:
pass
acme = all_consts.copy() # all constellations minus edge cases
wc = {} # wrapping about ra edge cases
ac = {} # amended (split) wrapping constellations
pc = {} # polar edge cases
for const_name in acme.keys():
if(const_name in pcl):
pc[const_name]=acme[const_name]
del acme[const_name]
if(const_name in wcl):
wc[const_name]=acme[const_name]
del acme[const_name]
def rotate(l,n):
return l[n:] + l[:n]
def linear_interpolate(x0,y0,x1,y1,x=None,y=None):
# if not passed any x or y, assume want bisection pt
if(x==None and y==None):
x=(x0+x1)/2.
y=(y0+y1)/2.
elif(y==None):
y=y0+(y1-y0)*(x-x0)/(x1-x0)
elif(x==None):
x=x0+(x1-x0)*(y-y0)/(y1-y0)
return x,y
# amend pole cases
# this adds points to the constellation bounds that
# are based on the edges of the plotted field
for const_name in pc.keys():
# the pts have to be in order from one discontinuity pt to
# ending at the other continuity pt to use pc[const_name][0]
# and pc[const_name][-1] as end-points. since in the above data
# we see if we print these pts that they are NOT the discontinuity
# pts, we could take the largest jump between right ascensions
# between any two coordinates. Or alternately, can take any coordinate pair
# that the distance is shorter when wrapping than direct. I will do the later.
# this too could have errors if asterisms were given complicated bounds
# at the greatest and least right ascension bounds, but thankfully that
# isnt the case here. We'll have to keep an eye out for these
# fringe cases in future asterisms. The solution I'll use at that point
# is to always arrange the pts with wrapping polar discontinuities at
# beggining and end when I have control over the input.
# The following rotates the constellation pts until first and last indeed
# line up like this.
num_coords = len(pc[const_name])
for i in range(num_coords):
try:
ra_one = pc[const_name][i][0]
ra_two = pc[const_name][i+1][0]
second_index = i+1
except:
ra_one = pc[const_name][i][0]
ra_two = pc[const_name][0][0]
second_index = 0
delta_ra_close = abs(ra_two-ra_one)
delta_ra_wrap = (bounds[1]-ra_one)+(ra_two-bounds[0])
if(delta_ra_wrap<delta_ra_close):
#print('Found the discontinuity in {} ({}<{}).\nPt{}{} and Pt{}{}'.format(const_name.upper(),delta_ra_wrap,delta_ra_close,i,pc[const_name][i],second_index,pc[const_name][second_index]))
nr = second_index
rotated = rotate(pc[const_name],nr)
(estimated_ra2,estimated_dec2) = \
linear_interpolate(rotated[-1][0]-bounds[1],rotated[-1][1],rotated[0][0],rotated[0][1],x=bounds[0])
(estimated_ra1,estimated_dec1) = (bounds[0],estimated_dec2)
est_pt1 = (estimated_ra1,estimated_dec1)
est_pt2 = (estimated_ra2+bounds[1],estimated_dec2)
if(pc[const_name][0][1]>0):
polar_pt1 = (bounds[0],bounds[3])
polar_pt2 = (bounds[1],bounds[3])
else:
polar_pt1 = (bounds[0],bounds[2])
polar_pt2 = (bounds[1],bounds[2])
rotated.insert(0,est_pt1)
rotated.insert(0,polar_pt1)
rotated.append(est_pt2)
rotated.append(polar_pt2)
pc[const_name] = rotated
# amend wrapping cases
# this splits the constellations and linearly interpolates between
# the discontinuity pts to add amending edge pts.
for const_name in wc.keys():
num_coords = len(wc[const_name])
# NOTE: an assumption I am making is that none of the edge constellations are
# wide enough to cross the 12h vertical line except the polar constellations.
plot_center = ((bounds[0]+bounds[1])/2.,((bounds[2]+bounds[3])/2.))
# We will find the discontinuities similarly to the above polar method.
# We'll rotate the constellation pts until we can make a clean cut at the discontinuities
# and then interpolate to the border and add the required extra pts.
# This technique is currently also assuming that the constellations do not cross
# back and forth between the borders. We'll have to account for that
# when it comes up in the future.
nr = []
for i in range(num_coords):
try:
ra_one = wc[const_name][i][0]
ra_two = wc[const_name][i+1][0]
second_index = i+1
except:
ra_one = wc[const_name][i][0]
ra_two = wc[const_name][0][0]
second_index = 0
if(ra_one>ra_two):
delta_ra_close = abs(ra_two-ra_one)
delta_ra_wrap = (bounds[1]-ra_one)+(ra_two-bounds[0])
else:
delta_ra_close = abs(ra_one-ra_two)
delta_ra_wrap = (bounds[1]-ra_two)+(ra_one-bounds[0])
if(delta_ra_wrap<delta_ra_close):
#print('Found a discontinuity in {} ({}<{}).\nPt{}{} and Pt{}{}'.format(const_name.upper(),delta_ra_wrap,delta_ra_close,i,wc[const_name][i],second_index,wc[const_name][second_index]))
nr.append(second_index)
#print('Total of {} discontinuities found in {} at indexes {}'.format(len(nr),const_name.upper(),nr))
if(len(nr)%2!=0):
print('WARNING: there should be an even amount of discontinuities in {}.'.format(const_name.upper()))
rotated = rotate(wc[const_name],nr[0])
# If more than two discontinuities, would recursively break up the constellation
# but confirmed only two discontinuities for all of these IAU bounds.
# One discontinuity is now across the 0th and -1th indexes at beginning and
# end of list. The other is at the discontinuity index, di.
di = nr[1]-nr[0]
l = rotated
if(l[0][0]>plot_center[0]):
left_slice = l[di:]
right_slice = l[:di]
else:
right_slice = l[di:]
left_slice = l[:di]
est_pt1 = linear_interpolate(right_slice[-1][0]-bounds[1],right_slice[-1][1],left_slice[0][0],left_slice[0][1],x=bounds[0])
est_pt2 = linear_interpolate(right_slice[0][0]-bounds[1],right_slice[0][1],left_slice[-1][0],left_slice[-1][1],x=bounds[0])
left_slice.insert(0,est_pt1)
left_slice.append(est_pt2)
est_pt3 = (bounds[1],est_pt2[1])
est_pt4 = (bounds[1],est_pt1[1])
right_slice.insert(0,est_pt3)
right_slice.append(est_pt4)
# amended constellations
ac[const_name+'1'] = right_slice
ac[const_name+'2'] = left_slice
del wc[const_name]
# write over all_consts by merging the separated and amended constellations
all_consts = dict(ac.items()+pc.items()+acme.items())
If we uncomment the discontinuity print line, we see that the wrapping constellation detection is much more certain than using the equivalent for the polar coordinates. Ursa Minor was sortof close to failing our criteria, compared to all wrapping cases being extremely certain. We also see that we must deal with two discontinuities at minimum and whether the pts are going from left-to-right or right-to-left across the discontinuity from higher-to-lower RA.
from matplotlib import pyplot
from shapely.geometry import LineString,Point,LinearRing
def plot_line(ax, ob):
x, y = ob.xy
ax.plot(x, y, color=BLUE, linewidth=3, solid_capstyle='round', zorder=1)
def plot_pt(ax, ob):
x, y = ob.xy
ax.plot(x, y, 'o',color=BASE02, zorder=1)
def set_colors():
# I am using the solarized color scheme out of personal preference.
# If anyone has ideas on a scientific scaling of color, I would love
# to hear.
global BASE03, BASE02, BASE01, BASE00, BASE0, BASE1, BASE2, BASE3, YELLOW, ORANGE, RED, MAGENTA, VIOLET, BLUE, CYAN, GREEN
BASE03 = '#002b36'
BASE02 = '#073642'
BASE01 = '#586e75'
BASE00 = '#657b83'
BASE0 = '#839496'
BASE1 = '#93a1a1'
BASE2 = '#eee8d5'
BASE3 = '#fdf6e3'
YELLOW = '#b58900'
ORANGE = '#cb4b16'
RED = '#dc322f'
MAGENTA = '#d33682'
VIOLET = '#6c71c4'
BLUE = '#268bd2'
CYAN = '#2aa198'
GREEN = '#859900'
def const_color(name=None,group=0,return_labels=False):
# If you have ideas for alternate groupings I would love to hear in the
# issue tracker.
families = {
BASE01:('zodiac',['aqr','cap','sgr','sco','lib','vir','leo','cnc','gem','tau','ari','psc']),
BLUE:('ursa major',['dra','uma','umi','lyn','lmi','cvn','crb','com','aql','boo','cam']),
RED:('perseus',['cep','peg','aur','per','and','tri','cas','cet','lac']),
GREEN:('hercules',['cyg','lyr','her','sge','ser1','cru','ser2','hya','cra','ara','sct','oph','lup','tra','cen','sex','crt','crv','vul']),
ORANGE:('orion',['cma','cmi','mon','lep','ori']),
CYAN:('heavenly waters',['del','equ','vel','car','pup','pyx','eri','col','psa']),
YELLOW:('bayer',['tuc','gru','ind','pav','aps','mus','cha','vol','dor','hyi','phe']),
VIOLET:('la caille',['mic','oct','nor','tel','men','ret','for','scl','cae','hor','ant','pic','cir']),
}
origin = {
BASE02:('ancient',['psa','equ','del','and','aql','aqr','ara','ari',
'cap','cas','cen','cep','cet','cma','cmi','cnc','cra','crb','crt',
'crv','cyg','dra','eri','gem','her','hya','leo','lep','lib','lup',
'lyr','oph','ori','peg','per','psc','sco','ser1','ser2','sge','sgr',
'tau','tri','uma','umi','vir','aur','boo']),
MAGENTA:('1592',['col']),
RED:('1603',['com','vol','tuc','phe','hyi','dor','cha','cru','mus','tra','aps','pav','ind','gru']),
YELLOW:('1613',['cam','mon',]),
GREEN:('1692',['lmi','lyn','sex','vul','cvn','lac','sct']),
BLUE:('1763',['oct','mic','for','scl','cae','hor','pic','ret','men','ant','nor','cir','tel','vel','pyx','pup','car',]),
}
if(group==0):
color_by=families
elif(group==1):
color_by=origin
elif(group==2):
# BY DENSITY
#color_by=color_by_density()
pass
else:
color_by=families
if(return_labels):
for color in color_by.keys():
try:
label_list.append((color,color_by[color][0]))
except:
label_list = [(color,color_by[color][0])]
return label_list
else:
for color in color_by.keys():
for const in color_by[color][1]:
# use startswith because it will be tolerant of splitting constellations if need be
if(name.startswith(const)):
return color
return BASE03
def plot_const(name='all',show_scatter=False,constellations=all_consts,plot_pts=False,color=False,group=0,legend=False,labels=False):
set_colors()
if(show_scatter):
scatter(all_ras,all_decs)
axis(bounds)
else:
fig = pyplot.figure(figsize=(12,9))
ax = fig.add_subplot(111)
for const_name in constellations.keys():
if(const_name==name or name=='all'):
ring = LinearRing(constellations[const_name])
if(color==True):
polygon = Polygon(ring, fc=const_color(const_name,group=group), ec=const_color(const_name,group=group), alpha=0.5, zorder=2)
ax.add_patch(polygon)
else:
plot_line(ax,ring)
if(plot_pts):
plot_pt(ax,ring)
if(name=='all'):
ax.set_title('All IAU Constellation Bounds')
axis(bounds)
plt.xticks([0,3,6,9,12,15,18,21,24])
plt.yticks([-90,-60,-30,0,30,60,90])
else:
ax.set_title('{} Bounds'.format(name.upper()))
if(legend):
color_label_list = const_color(group=group,return_labels=True)
for i in range(len(color_label_list)):
plt.plot(None,None,color=color_label_list[i][0], linewidth=4.0, linestyle="-", label=color_label_list[i][1])
if(len(color_label_list)>0):
plt.legend(loc='upper left')
if(labels):
plt.xlabel('Declination')
plt.ylabel('Right Ascension')
gca().invert_xaxis()
gca().xaxis.set_major_formatter(hours_formatter)
gca().yaxis.set_major_formatter(degrees_formatter)
gca().xaxis.grid(True)
gca().yaxis.grid(True)
You can plot individual constellation bounds by passing the abbreviation of the constellation to plot_const()
.
plot_const('ori')
When plotting everything, the constellations are broken up into three groups: all-constellations-minus-edges, wrapping-constellations, and polar-constellations.
plot_const(constellations=acme,labels=True)
plot_const(constellations=pc,labels=True)
plot_const(constellations=ac,labels=True)
plot_const(color=True,group=0,legend=True)
#savefig('bounds_family.png')
plot_const(color=True,group=1,legend=True)
#savefig('bounds_time.png')
Pierre Barbier brings up some interesting history on these bounds when looking for the original data from 1930. He found that the IAU listed file bound_18.dat is close to the original but converted, truncated, and vertices were added to aid plotting. He has made the original data available on his website compiled directly from Delporte's Délimitation Scientifique des Constellations.
The bound_18.dat
file mentioned by Pierre Barbier refers to the file associated with A Catalogue of Constellation Boundary Data (CDS Archive, PDF) by A.C. Davenhall and S.K. Leggett in 1990. They give the boundaries delineated by Delporte in 1930 which is oriented to the 1875 epoch: bound_18.dat. Then constbnd.dat includes the boundaries for B1875 with adjacent constellations and constbnd.txt gives a thorough explanation of the file. They also give a third file with the boundaries computed for equator and equinox 2000: bound_20.dat. O
= original, I
= interpolated.
In equator and equinox 1875 the lines joining the corners of the constellations were great circles of right ascension and parallels of declination. However, precession to equator and equinox 2000 distorts the boundaries so that they no longer lie along lines of constant right ascension or declination. Thus, prior to precession to equator and equinox 2000, points were interpolated at one degree intervals along the boundaries in order that they should continue to enclose the same area of sky (and thus the same set of stars).