http://docs.astropy.org/en/stable/api/astropy.wcs.WCS.html#astropy.wcs.WCS
    naxis=[0,1] is to take x,y axes in case a color photo was input e.g. to astrometry.net cloud solver

def fits2radec(fitsfn: Path, solve: bool=False, args: str=None) -> xarray.Dataset:

    fitsfn = Path(fitsfn).expanduser()

    if fitsfn.suffix == '.fits':
        WCSfn = fitsfn.with_suffix('.wcs')  # using .wcs will also work but gives a spurious warning
    elif fitsfn.suffix == '.wcs':
        WCSfn = fitsfn
    else:
        raise ValueError(f'please convert {fitsfn} to GRAYSCALE .fits e.g. with ImageJ or ImageMagick')

    if solve:
        doSolve(fitsfn, args)

    with fits.open(fitsfn, mode='readonly') as f:
        yPix, xPix = f[0].shape[-2:]

    x, y = meshgrid(range(xPix), range(yPix))  # pixel indices to find RA/dec of
    xy = column_stack((x.ravel(order='C'), y.ravel(order='C')))
# %% use astropy.wcs to register pixels to RA/DEC
    """
    http://docs.astropy.org/en/stable/api/astropy.wcs.WCS.html#astropy.wcs.WCS
    naxis=[0,1] is to take x,y axes in case a color photo was input e.g. to astrometry.net cloud solver
    """
    try:
        with fits.open(WCSfn, mode='readonly') as f:
            # radec = wcs.WCS(hdul[0].header,naxis=[0,1]).all_pix2world(xy, 0)
            radec = wcs.WCS(f[0].header).all_pix2world(xy, 0)
    except OSError:
        raise OSError(f'It appears the WCS solution is not present, was the FITS image solved?  looking for: {WCSfn}')

    ra = radec[:, 0].reshape((yPix, xPix), order='C')
    dec = radec[:, 1].reshape((yPix, xPix), order='C')
# %% collect output
    radec = xarray.Dataset({'ra': (('y', 'x'), ra),
                            'dec': (('y', 'x'), dec), },
                           {'x': range(xPix), 'y': range(yPix)},
                           attrs={'filename': fitsfn})

    return radec