code: def wcs_to_coords(w, shape):
      
      if w.naxis == 2:
          y, x = wcs_to_axes(w, shape)
      elif w.naxis == 3:
          z, y, x = wcs_to_axes(w, shape)
      else:
          raise Exception("Wrong number of WCS axes %i" % w.naxis)
  
      x = 0.5 * (x[1:] + x[:-1])
      y = 0.5 * (y[1:] + y[:-1])
  
      if w.naxis == 2:
          x = np.ravel(np.ones(shape) * x[:, np.newaxis])
          y = np.ravel(np.ones(shape) * y[np.newaxis, :])
          return np.vstack((x, y))
  
      z = 0.5 * (z[1:] + z[:-1])
      x = np.ravel(np.ones(shape) * x[:, np.newaxis, np.newaxis])
      y = np.ravel(np.ones(shape) * y[np.newaxis, :, np.newaxis])
      z = np.ravel(np.ones(shape) * z[np.newaxis, np.newaxis, :])
  
      return np.vstack((x, y, z))
docs: Generate an N x D list of pixel center coordinates where N is
      the number of pixels and D is the dimensionality of the map.