From 632b27382673941b5d83d1c9ae5ea206ded42907 Mon Sep 17 00:00:00 2001 From: scivision Date: Mon, 5 Feb 2024 22:39:57 -0500 Subject: [PATCH] plot unwarps sources also --- scripts/OverlayAltitudes.py | 6 +- scripts/OverlayStars.py | 95 ++++++++++++++++++---------- src/astrometry_azel/plot/__init__.py | 29 ++++++++- 3 files changed, 91 insertions(+), 39 deletions(-) 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.]")