diff --git a/scripts/OverlayStars.py b/scripts/OverlayStars.py
index a7646c4..36b8a6f 100644
--- a/scripts/OverlayStars.py
+++ b/scripts/OverlayStars.py
@@ -3,10 +3,16 @@
 """
 Overlay stars from a catalog on an image.
 This is like a closed loop check that solve-field output (here, the .rdls file) makes sense.
+
+Use --tag-all option of solve-field to output star magnitude in the .rdls file:
+
+    solve-field src/astrometry_azel/tests/apod4.fits --tag-all
 """
 
 import argparse
 
+import numpy as np
+
 from astrometry_azel.io import get_sources
 from astrometry_azel.plot import add_image
 
@@ -23,8 +29,39 @@
 fg = figure()
 ax = fg.gca()
 add_image(p.new_file, "gray", ax)
-ax.set_title(f"{p.rdls_file} stars on {p.new_file}", wrap=True)
 
-ax.scatter(source_ra["ra"], source_ra["dec"], s=10, c="r", marker=".", label="stars")
+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]
+
+    # 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}"
+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)
 
 show()
diff --git a/src/astrometry_azel/io.py b/src/astrometry_azel/io.py
index 226f480..e37de3a 100755
--- a/src/astrometry_azel/io.py
+++ b/src/astrometry_azel/io.py
@@ -29,7 +29,7 @@
     loadmat = None
 
 
-def get_sources(fn: Path) -> xarray.Dataset:
+def get_sources(fn: Path):
     """
     read source (star) coordinates from .rdls file
 
@@ -42,7 +42,7 @@ def get_sources(fn: Path) -> xarray.Dataset:
     Returns
     -------
 
-    radec: xarray.Dataset
+    radec: astropy.fits.fitsrec.FITS_rec
         RA, Dec of sources
     """