# 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() 

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()