docs: def surveyPixel(lon, lat, nside_pix, nside_subpix = None):
      """
      Return the set of HEALPix pixels that cover the given coordinates at resolution nside.
      Optionally return the set of subpixels within those pixels at resolution nside_subpix
      """
      pix = np.unique(ang2pix(nside_pix, lon, lat))
      if nside_subpix is None:
          return pix
      else:
          subpix_array = []
          for ii in range(0, len(pix)):
              subpix = subpixel(pix[ii], nside_pix, nside_subpix)
              subpix_array.append(subpix)
          return pix, np.array(subpix_array)