Kapteyn Institute Kapteyn Package

Navigation

All sky plots and graticules

All Sky plots

An all sky plot is a plot where the range in longitude is [0,360> and the range in latitude is [-90,90>. There are many examples in Calabretta’s article Representations of celestial coordinates in FITS We tried to reproduce these figures both to prove that the modules in the Kapteyn Package have the functionality to do it and to facilitate users who want to set up an all-sky plot with module maputils. For this purpose we created for each figure the minimal required FITS headers. Header and other code is listed below the examples. In the HTML documentation, click on the hires link to get a plot which shows more details. With the information in this document, it should be easy to compose a Python program that creates just a single plot which then can be enhanced to fit your needs.

The first plot is a stand alone version. The others are generated with different Python scripts and the service module service.py (see also the source code at the end of this document).

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
from kapteyn import maputils
from numpy import arange
from matplotlib import pyplot as plt

dec0 = 89.9999999999   # Avoid plotting on the wrong side
header = {'NAXIS'  :  2,
          'NAXIS1' :  100, 'NAXIS2': 80,
          'CTYPE1' : 'RA---TAN',
          'CRVAL1' :  0.0, 'CRPIX1' : 50, 'CUNIT1' : 'deg',
          'CDELT1' : -5.0, 'CTYPE2' : 'DEC--TAN',
          'CRVAL2' :  dec0, 'CRPIX2' : 40,
          'CUNIT2' : 'deg', 'CDELT2' : 5.0,
         }
X = arange(0,360.0,15.0)
Y = [20, 30,45, 60, 75]

fig = plt.figure(figsize=(7,6))
frame = fig.add_axes((0.1,0.1,0.8,0.8))
f = maputils.FITSimage(externalheader=header)
annim = f.Annotatedimage(frame)
grat = annim.Graticule(wylim=(20.0,90.0), wxlim=(0,360),
                       startx=X, starty=Y)
grat.setp_gratline(color='0.75')
lon_world = range(0,360,30)
lat_world = [20, 30, 60, 90]
grat.setp_lineswcs1(20, color='g', linestyle='--')

# Plot labels inside the plot
lon_constval = None
lat_constval = 20
lon_kwargs = {'color':'r', 'fontsize':15}
lat_kwargs = {'color':'b', 'fontsize':10}
grat.Insidelabels(wcsaxis=0,
                  world=lon_world, constval=lat_constval,
                  fmt="Dms",
                  **lon_kwargs)
grat.Insidelabels(wcsaxis=1,
                  world=lat_world, constval=lon_constval,
                  fmt="Dms",
                  **lat_kwargs)
annim.plot()
# Set title for Matplotlib
titlepos = 1.02
title = r"Gnomonic projection (TAN) diverges at $\theta=0^\circ$. (Cal. fig.8)"
t = frame.set_title(title, color='g')
t.set_y(titlepos)
plt.show()

(Source code, png, hires.png, pdf)

_images/allsky_single.png

The recipe

The next session shows a gallery of all sky plots, all based on the same recipe.

  • One starts with a self-made header which ensures a complete coverage of the sky by stretching the values of the CDELT‘s.
  • Then an object from class maputils.Annotatedimage.Graticule is created with explicit limits for the world coordinates in both directions.
  • For these plots we don’t have intersections of the graticule with an enclosing rectangle so we cannot plot standard axis labels for the coordinates. Instead we use method wcsgrat.Graticule.Insidelabels() to plot labels inside the plot. In the plots we show different examples how one can manipulate these labeling.

Source code of the service module

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
from kapteyn import maputils, tabarray
import numpy
import sys
from matplotlib import pyplot as plt


__version__ = '1.91'

epsilon = 0.0000000001

def radians(a):
   return a*numpy.pi/180.0

def degrees(a):
   return a*180.0/numpy.pi

def cylrange():
   X = numpy.arange(0,400.0,30.0); 
   # Replace last two (dummy) values by two values around 180 degrees   
   X[-1] = 180.0 - epsilon
   X[-2] = 180.0 + epsilon
   return X

def polrange(): 
   X = numpy.arange(0,380.0,15); 
   # Replace last two (dummy) values by two values around 180 degrees   
   X[-1] = 180.0 - epsilon
   X[-2] = 180.0 + epsilon
   return X


def getperimeter(grat):
   # Calculate perimeter of QUAD projection 
   xlo, y = grat.gmap.topixel((-45.0-epsilon, 0.0))
   xhi, y = grat.gmap.topixel((315+epsilon, 0.0))
   x, ylo = grat.gmap.topixel((180, -45))
   x, yhi = grat.gmap.topixel((180, 45))
   x1, y = grat.gmap.topixel((45-epsilon, 0.0))
   x, ylolo = grat.gmap.topixel((0, -135+epsilon))
   x, yhihi = grat.gmap.topixel((0, 135-epsilon))
   perimeter = [(xlo,ylo), (x1,ylo), (x1,ylolo), (xhi,ylolo), (xhi,yhihi),
                (x1,yhihi), (x1,yhi), (xlo,yhi), (xlo,ylo)]
   return perimeter


def plotcoast(fn, frame, grat, col='k', lim=100, decim=5, plotsym=None, sign=1.0):
   coasts = tabarray.tabarray(fn, comchar='s')  # Read two columns from file
   for segment in coasts.segments:
      coastseg = coasts[segment].T
      xw = sign * coastseg[1]; yw = coastseg[0] # First one appears to be Latitude
      xs = xw; ys = yw                          # Reset lists which store valid pos.
      if 1:
         # Mask arrays if outside plot box
         xp, yp = grat.gmap.topixel((numpy.array(xs),numpy.array(ys)))
         # Be sure you understand 
         # the operator precedence: (a > 2) | (a < 5) is the proper syntax 
         # because a > 2 | a < 5 will result in an error due to the fact 
         # that 2 | a is evaluated first.
         xp = numpy.ma.masked_where(numpy.isnan(xp) |
                              (xp > grat.pxlim[1]) | (xp < grat.pxlim[0]), xp)
         yp = numpy.ma.masked_where(numpy.isnan(yp) |
                              (yp > grat.pylim[1]) | (yp < grat.pylim[0]), yp)
         # Mask array could be of type numpy.bool_ instead of numpy.ndarray
         if numpy.isscalar(xp.mask):
            xp.mask = numpy.array(xp.mask, 'bool')
         if numpy.isscalar(yp.mask):
            yp.mask = numpy.array(yp.mask, 'bool')
         # Count the number of positions in this list that are inside the box
         xdc = []; ydc = []
         for i in range(len(xp)):
            if not xp.mask[i] and not yp.mask[i]:
               if not i%decim:
                  xdc.append(xp.data[i])
                  ydc.append(yp.data[i])
         if len(xdc) >= lim:
            if plotsym == None:
               frame.plot(xdc, ydc, color=col)
            else:
               frame.plot(xdc, ydc, '.', markersize=1, color=col)



# Set defaults which can be overwritten by the allskyfxx.py scripts
title = ''
titlepos = 1.02
dec0 = 89.9999999999
fsize = 10
figsize = (7,6)
drawgrid = False
grat = None
smallversion = False
plotbox = (0.1,0.05,0.8,0.8)
markerpos = "120 deg 60 deg"


def doplot(frame, fignum, annim, grat, title, 
           lon_world=None, lat_world=None, 
           lon_constval=None, lat_constval=None,
           lon_fmt=None, lat_fmt=None,
           markerpos=None, 
           plotdata=False, perimeter=None, drawgrid=None, 
           smallversion=False, addangle0=0.0, addangle1=0.0, 
           framebgcolor=None, deltapx0=0.0, deltapy0=0.0,
           deltapx1=0.0, deltapy1=0.0,
           labkwargs0={'color':'r'}, labkwargs1={'color':'b'}):
# Apply some extra settings
   
   if framebgcolor != None:
      frame.set_axis_bgcolor(framebgcolor)

   # Plot coastlines if required
   if plotdata:
      if fignum == 36:
         plotcoast('WDB/world.txt', frame, grat, col='k', lim=100)
      else:
         plotcoast('WDB/world.txt', frame, grat, col='r', lim=50, 
                    decim=20, plotsym=',', sign=-1)


   if lon_constval == None:
      lon_constval = 0.0    # Reasonable for all sky plots
   if lat_constval == None:
      lat_constval = 0.0    # Reasonable for all sky plots
   if lon_fmt == None:
      lon_fmt = 'Dms'
   if lat_fmt == None:
      lat_fmt = 'Dms'
   # Plot labels inside graticule if required
   #labkwargs0.update({'fontsize':fsize})
   #labkwargs1.update({'fontsize':fsize})
   ilabs1 = grat.Insidelabels(wcsaxis=0, 
                        world=lon_world, constval=lat_constval, 
                        deltapx=deltapx0, deltapy=deltapy0, 
                        addangle=addangle0, fmt=lon_fmt, **labkwargs0)
   ilabs2 = grat.Insidelabels(wcsaxis=1, 
                        world=lat_world, constval=lon_constval, 
                        deltapx=deltapx1, deltapy=deltapy1, 
                        addangle=addangle1, fmt=lat_fmt, **labkwargs1)

   # Plot just 1 pixel c.q. marker
   if markerpos != None:
      annim.Marker(pos=markerpos, marker='o', color='red')

   if drawgrid:
      pixellabels = annim.Pixellabels(plotaxis=(2,3))

   # Plot the title
   title = "Fig. %d: "%fignum + title
   if smallversion:
      t = frame.set_title(title, color='g', fontsize=10)
   else:
      t = frame.set_title(title, color='g', fontsize=13, linespacing=1.5)
   t.set_y(titlepos)
   annim.plot()
   
   # Plot alternative borders. Do this after the graticule is plotted
   # Only then you know the frame of the graticule and plotting in that
   # frame will overwrite graticule lines so that the borders look better
   if perimeter != None:
      p = plt.Polygon(perimeter, facecolor='#d6eaef', lw=2)
      frame.add_patch(p)     # Must be in frame specified by user
      Xp, Yp = zip(*perimeter)
      grat.frame.plot(Xp, Yp, color='r')

   annim.interact_toolbarinfo()

   plt.show()

Table Of Contents

Please cite the Kapteyn Package if you use it in the preparation of a publication. These citations help us to justify the resources spent on this software.

Navigation

© Copyright 2010-2013, Kapteyn Astronomical Institute. Last updated on Jan 04, 2015. Created using Sphinx 1.2.1.