diff --git a/scripts/OverlayAltitudes.py b/scripts/OverlayAltitudes.py
index 0ca2820..c017f5e 100755
--- a/scripts/OverlayAltitudes.py
+++ b/scripts/OverlayAltitudes.py
@@ -17,7 +17,7 @@
 from pathlib import Path
 from argparse import ArgumentParser
 
-from astrometry_azel.plot import add_image
+from astrometry_azel.plot import wcs_image
 
 from matplotlib.pyplot import figure, show
 
@@ -37,11 +37,11 @@
 if p.subplots:
     axs: typing.Any = fg.subplots(1, len(flist), sharey=True, sharex=True)
     for fn, ax in zip(flist, axs):
-        add_image(fn, "gray", ax)
+        wcs_image(fn, "gray", ax)
 
 else:
     ax = fg.gca()
     for fn, cm in zip(flist, cmaps):
-        add_image(fn, cm, ax, alpha=0.05)
+        wcs_image(fn, cm, ax, alpha=0.05)
 
 show()
diff --git a/scripts/OverlayStars.py b/scripts/OverlayStars.py
index 36b8a6f..ed0c6f3 100644
--- a/scripts/OverlayStars.py
+++ b/scripts/OverlayStars.py
@@ -10,58 +10,87 @@
 """
 
 import argparse
+from pathlib import Path
 
 import numpy as np
 
 from astrometry_azel.io import get_sources
-from astrometry_azel.plot import add_image
+from astrometry_azel.plot import wcs_image, xy_image
 
 from matplotlib.pyplot import figure, show
 
 
+def label_stars(ax, x, y, mag_norm=None):
+    if mag_norm is not None:
+        ax.scatter(
+            x,
+            y,
+            s=100 * mag_norm,
+            edgecolors="red",
+            marker="o",
+            facecolors="none",
+            label="stars",
+        )
+    else:
+        ax.scatter(
+            x,
+            y,
+            s=25,
+            edgecolors="red",
+            marker="o",
+            facecolors="none",
+            label="stars",
+        )
+
+
 p = argparse.ArgumentParser()
-p.add_argument("rdls_file", help="FITS .rdls filename output from solve-field (Astrometry.net)")
-p.add_argument("new_file", help="FITS .new filename output from solve-field (Astrometry.net)")
+p.add_argument(
+    "stem", help="FITS file stem (without .suffix) output from solve-field (Astrometry.net)"
+)
 p = p.parse_args()
 
-source_ra = get_sources(p.rdls_file)
+stem = Path(p.stem)
+new = stem.with_suffix(".new")
+
+rdls = stem.with_suffix(".rdls")
+source_ra = get_sources(rdls)
+
 
-fg = figure()
-ax = fg.gca()
-add_image(p.new_file, "gray", ax)
+xyls = stem.parent / (stem.name + "-indx.xyls")
+source_xy = get_sources(xyls)
 
 if "MAG" in source_ra.columns.names:
     Ntot = source_ra.shape[0]
     Nkeep = 20  # only plot the brightest stars
-    source_ra.sort(axis=0, order="MAG")
-    source_ra = source_ra[:Nkeep]
+
+    i = source_ra.argsort(axis=0, order="MAG")
+
+    source_ra = np.take_along_axis(source_ra, i, axis=0)[:Nkeep]
+    source_xy = np.take_along_axis(source_xy, i, axis=0)[:Nkeep]
 
     # normalize star magnitude to [0,1]
     mag_norm = (source_ra["MAG"] - source_ra["MAG"].min()) / np.ptp(source_ra["MAG"])
-
-    ax.scatter(
-        source_ra["RA"],
-        source_ra["DEC"],
-        s=100 * mag_norm,
-        edgecolors="red",
-        marker="o",
-        facecolors="none",
-        label="stars",
-    )
-    ttxt = f"{p.rdls_file} {source_ra.shape[0]} / {Ntot} stars on {p.new_file}"
+    ttxt = f"{stem}: {source_ra.shape[0]} / {Ntot} stars"
 else:
-    ax.scatter(
-        source_ra["RA"],
-        source_ra["DEC"],
-        s=25,
-        edgecolors="red",
-        marker="o",
-        facecolors="none",
-        label="stars",
-    )
-    ttxt = f"{p.rdls_file} {source_ra.shape[0]} stars on {p.new_file}"
-
-
-ax.set_title(ttxt, wrap=True)
+    mag_norm = None
+    ttxt = f"{stem} {source_ra.shape[0]} stars"
+
+# %% unwarped image
+
+fgy = figure()
+axy = fgy.gca()
+xy_image(new, "gray", axy)
+
+label_stars(axy, source_xy["X"], source_xy["Y"], mag_norm)
+axy.set_title(ttxt, wrap=True)
+
+# %% WCS warped image
+
+fgr = figure()
+axr = fgr.gca()
+wcs_image(new, "gray", axr)
+
+label_stars(axr, source_ra["RA"], source_ra["DEC"], mag_norm)
+axr.set_title(ttxt, wrap=True)
 
 show()
diff --git a/src/astrometry_azel/plot/__init__.py b/src/astrometry_azel/plot/__init__.py
index 4f06d80..57a2a95 100644
--- a/src/astrometry_azel/plot/__init__.py
+++ b/src/astrometry_azel/plot/__init__.py
@@ -175,11 +175,12 @@ def image_stack(img, fn: Path, clim=None):
     fg.savefig(plotFN)
 
 
-def add_image(fn: Path, cm, ax, alpha=1):
+def wcs_image(fn: Path, cmap, ax, alpha=1):
     """
     Astrometry.net makes file ".new" with the image and the WCS SIP 2-D polynomial fit coefficients in the FITS header
 
-    We use DECL as "x" and RA as "y".
+    Warps the image to sky coordinates using the WCS.
+
     pcolormesh() is used as it handles arbitrary pixel shapes.
     Note that pcolormesh() cannot tolerate NaN in X or Y (NaN in C is OK).
 
@@ -201,6 +202,28 @@ def add_image(fn: Path, cm, ax, alpha=1):
     dec = radec[:, 1].reshape((yPix, xPix), order="C")
 
     ax.set_title(fn.name)
-    ax.pcolormesh(ra, dec, img, alpha=alpha, cmap=cm, norm=LogNorm())
+    ax.pcolormesh(ra, dec, img, alpha=alpha, cmap=cmap, norm=LogNorm())
+    ax.set_ylabel("Right Ascension [deg.]")
+    ax.set_xlabel("Declination [deg.]")
+
+
+def xy_image(fn: Path, cmap, ax):
+    """
+    Astrometry.net makes file ".new" with the image and the WCS SIP 2-D polynomial fit coefficients in the FITS header
+
+    We use DECL as "x" and RA as "y".
+    pcolormesh() is used as it handles arbitrary pixel shapes.
+    Note that pcolormesh() cannot tolerate NaN in X or Y (NaN in C is OK).
+
+    https://github.com/scivision/python-matlab-examples/blob/main/PlotPcolor/pcolormesh_NaN.py
+    """
+
+    fn = Path(fn).expanduser().resolve(True)
+
+    with fits.open(fn, mode="readonly", memmap=False) as f:
+        img = f[0].data
+
+    ax.set_title(fn.name)
+    ax.pcolormesh(img, cmap=cmap, norm=LogNorm())
     ax.set_ylabel("Right Ascension [deg.]")
     ax.set_xlabel("Declination [deg.]")