def read_partial_map(filenames, column, fullsky=True, **kwargs):
    """
    Read a partial HEALPix file(s) and return pixels and values/map. Can
    handle 3D healpix maps (pix, value, zdim). Returned array has
    shape (dimz,npix).

    Parameters:
    -----------
    filenames     : list of input filenames
    column        : column of interest
    fullsky       : partial or fullsky map
    kwargs        : passed to fitsio.read

    Returns:
    --------
    (nside,pix,map) : pixel array and healpix map (partial or fullsky)
    """
    # Make sure that PIXEL is in columns
    #kwargs['columns'] = ['PIXEL',column]
    kwargs['columns'] = ['PIXEL'] + np.atleast_1d(column).tolist()

    filenames = np.atleast_1d(filenames)
    header = fitsio.read_header(filenames[0],ext=kwargs.get('ext',1))
    data = ugali.utils.fileio.load_files(filenames,**kwargs)

    pix = data['PIXEL']
    value = data[column]
    nside = header['NSIDE']
    npix = hp.nside2npix(nside)

    ndupes = len(pix) - len(np.unique(pix))
    if ndupes > 0:
        msg = '%i duplicate pixels during load.'%(ndupes)
        raise Exception(msg)

    if fullsky and not np.isscalar(column):
        raise Exception("Cannot make fullsky map from list of columns.")
    
    if fullsky:
        shape = list(value.shape)
        shape[0] = npix
        hpxmap = hp.UNSEEN * np.ones(shape,dtype=value.dtype)
        hpxmap[pix] = value
        return (nside,pix,hpxmap.T)
    else:
        return (nside,pix,value.T)