title: Include WCS coordinates in a FITS image
body: <p>I have a 2D numpy array 'ZEA_N_sky' initially created as</p>
  
  <pre><code>n=128
  ZEA_N_sky=np.zeros((n,n))
  </code></pre>
  
  <p>and later some values were assigned to every pixel of this array. I can plot this array as FITS image (that can be opened by ds9) using the following code, however I wish to see the WCS coordinates as well. That is, when I roll my mouse cursor over the image, I should be able to see the gal-long, gal-lat;  Right-Asc, Dec , etc. How can I achieve that ? Do i need to do it manually or there is some header() trick ?</p>
  
  <p>The documentation given using WCS package and use all_pix2world() does not seem to work (rather I dont seem to understand it). If somebody can help me please with the code? Thanks !</p>
  
  <pre><code>from astropy.io import fits
  import numpy as np
  from astropy import wcs
  from astropy.table import Table
  
  out_file_name = 'FITS_image.fits'
  
  hdr = fits.Header()
  hdr['Projection'] = "ZEA"
  hdr['nd_size'] = str(nd_size)
  hdr['SCALE']="nd_size/2"
  fits.writeto(out_file_name, ZEA_N_sky, hdr,clobber=True)
  </code></pre>
  
  <hr>
  
  <p>I have made some progress. This is the bit of code I am using right now using which I do obtain the WCS parameters, but the value of longitude that I am getting on my FITS image is wrong ! I think I am making mistake in the ZEA projection parameters values. </p>
  
  <pre><code>NSGP=1 #Parameter for North side projection, -1 for south side
  
  w = wcs.WCS(naxis=2)
  w.wcs.crpix = [SCALE,SCALE] # SCALE is half the value of my pixel range in the image. n/2=128/2=64
  w.wcs.cdelt = np.array([-NSGP*90.0/float(SCALE) * 0.90032, NSGP*90.0/float(SCALE) * 0.90032]) # increments in degrees per pixel
  w.wcs.crval = [NSGP*90.0, NSGP*90.0] #RA and dec values in hours and degrees    
  w.wcs.ctype = ["GLON-ZEA", "GLAT-ZEA"]
  #w.wcs.set_pv([(180, NSGP, float(SCALE))])
  
  out_file_name = 'N_Mateu_ZEA_stream_mask.fits'
  # Now, write out the WCS object as a FITS header
  header = w.to_header()
  hdu = fits.PrimaryHDU(ZEA_N_sky,header=header)
  hdu.writeto(out_file_name, clobber=True)
  </code></pre>