From a7fb47d7697e3f7b165d41eea973751b99076bc3 Mon Sep 17 00:00:00 2001 From: Marcus Hughes Date: Sat, 24 Aug 2024 02:00:23 -0600 Subject: [PATCH 1/4] incorporates masking and pv --- thuban/catalog.py | 7 ++++++- thuban/pointing.py | 45 ++++++++++++++++++++++++++++----------------- 2 files changed, 34 insertions(+), 18 deletions(-) diff --git a/thuban/catalog.py b/thuban/catalog.py index e05ec58..c1c8ea8 100644 --- a/thuban/catalog.py +++ b/thuban/catalog.py @@ -1,4 +1,5 @@ import os +from typing import Callable import astropy.units as u import numpy as np @@ -163,7 +164,7 @@ def filter_for_visible_stars(catalog: pd.DataFrame, dimmest_magnitude: float = 6 def find_catalog_in_image( - catalog: pd.DataFrame, wcs: WCS, image_shape: (int, int), + catalog: pd.DataFrame, wcs: WCS, image_shape: (int, int), mask: Callable = None, mode: str = "all" ) -> np.ndarray: """Using the provided WCS converts the RA/DEC catalog into pixel coordinates @@ -195,6 +196,10 @@ def find_catalog_in_image( except NoConvergence as e: xs, ys = e.best_solution[:, 0], e.best_solution[:, 1] bounds_mask = (0 <= xs) * (xs < image_shape[0]) * (0 <= ys) * (ys < image_shape[1]) + + if mask is not None: + bounds_mask *= mask(xs, ys) + reduced_catalog = catalog[bounds_mask].copy() reduced_catalog["x_pix"] = xs[bounds_mask] reduced_catalog['y_pix'] = ys[bounds_mask] diff --git a/thuban/pointing.py b/thuban/pointing.py index 42a27e1..9e15d37 100644 --- a/thuban/pointing.py +++ b/thuban/pointing.py @@ -1,9 +1,11 @@ import warnings +from typing import Callable import numpy as np import pandas as pd import sep_pjw as sep from astropy.wcs import WCS, utils +import astropy.units as u from lmfit import Parameters, minimize from thuban.catalog import (filter_for_visible_stars, find_catalog_in_image, @@ -47,19 +49,21 @@ def calculate_pc_matrix(crota: float, cdelt: (float, float)) -> np.ndarray: """ return np.array( [ - [np.cos(crota), np.sin(crota) * (cdelt[0] / cdelt[1])], - [-np.sin(crota) * (cdelt[1] / cdelt[0]), np.cos(crota)], + [np.cos(crota), -np.sin(crota) * (cdelt[0] / cdelt[1])], + [np.sin(crota) * (cdelt[1] / cdelt[0]), np.cos(crota)], ], ) - def _residual(params: Parameters, observed_coords: np.ndarray, catalog: pd.DataFrame, guess_wcs: WCS, n: int, - image_shape: (int, int)): + image_shape: (int, int), + edge: int = 100, + sigma: float = 3.0, + mask: Callable = None): """Residual used when optimizing the pointing Parameters @@ -86,13 +90,13 @@ def _residual(params: Parameters, refined_wcs.wcs.crpix = guess_wcs.wcs.crpix refined_wcs.wcs.crval = (params["crval1"].value, params["crval2"].value) refined_wcs.wcs.pc = calculate_pc_matrix(params["crota"], refined_wcs.wcs.cdelt) - refined_wcs.cpdis1 = guess_wcs.cpdis1 # TODO what if there is no distortion? + refined_wcs.cpdis1 = guess_wcs.cpdis1 refined_wcs.cpdis2 = guess_wcs.cpdis2 + refined_wcs.wcs.set_pv(guess_wcs.wcs.get_pv()) - reduced_catalog = find_catalog_in_image(catalog, refined_wcs, image_shape=image_shape) + reduced_catalog = find_catalog_in_image(catalog, refined_wcs, image_shape=image_shape, mask=mask) refined_coords = np.stack([reduced_catalog['x_pix'], reduced_catalog['y_pix']], axis=-1) - edge = 300 image_bounds = (image_shape[0] - edge, image_shape[1] - edge) refined_coords = np.array([c for c in refined_coords if (c[0] > edge) and (c[1] > edge) and (c[0] < image_bounds[0]) and (c[1] < image_bounds[1])]) @@ -104,10 +108,9 @@ def _residual(params: Parameters, "Try decreasing `num_stars` or increasing `dimmest_magnitude`.", RepeatedStarWarning) - sigma = 3 out = np.array([np.min(np.linalg.norm(observed_coords - coord, axis=-1)) for coord in refined_coords[:n]]) - mean, stdev = np.mean(out), np.std(out) - out[out > mean + sigma * stdev] = 0.0 + median, stdev = np.median(out), np.std(out) + out[out > median + (sigma * stdev)] = 0.0 # TODO: should this be zero? return out @@ -125,10 +128,16 @@ def refine_pointing_wrapper(image, guess_wcs, file_num, observed_coords=None, ca return new_wcs, observed_coords, solution, trial_num, file_num +def extract_crota_from_wcs(wcs: WCS) -> tuple[float, float]: + """Extract CROTA from a WCS.""" + delta_ratio = wcs.wcs.cdelt[1] / wcs.wcs.cdelt[0] + return np.arctan2(wcs.wcs.pc[1, 0]/delta_ratio, wcs.wcs.pc[0, 0]) * u.rad + + def refine_pointing(image, guess_wcs, observed_coords=None, catalog=None, background_width=16, background_height=16, detection_threshold=5, num_stars=30, max_trials=15, chisqr_threshold=0.1, - dimmest_magnitude=6.0, method='leastsq'): + dimmest_magnitude=6.0, method='leastsq', edge=100, sigma=3.0, mask=None): """ Refine the pointing for an image Parameters @@ -159,14 +168,15 @@ def refine_pointing(image, guess_wcs, observed_coords=None, catalog=None, background = sep.Background(image, bw=background_width, bh=background_height) data_sub = image - background objects = sep.extract(data_sub, detection_threshold, err=background.globalrms) - objects = pd.DataFrame(objects).sort_values('flux')[-3*num_stars:] + objects = pd.DataFrame(objects).sort_values('flux') observed_coords = np.stack([objects["x"], objects["y"]], axis=-1) - + if mask is not None: + observed_coords = observed_coords[mask(objects['x'], objects['y'])] + observed_coords = observed_coords[-3*num_stars:] # set up the optimization params = Parameters() - delta_ratio = guess_wcs.wcs.cdelt[0] / guess_wcs.wcs.cdelt[1] - initial_crota = np.arctan2(guess_wcs.wcs.pc[0, 1]/delta_ratio, guess_wcs.wcs.pc[0, 0]) - params.add("crota", value=initial_crota, min=-np.pi, max=np.pi) + initial_crota = extract_crota_from_wcs(guess_wcs) + params.add("crota", value=initial_crota.to(u.rad).value, min=-np.pi, max=np.pi) params.add("crval1", value=guess_wcs.wcs.crval[0], min=-180, max=180, vary=True) params.add("crval2", value=guess_wcs.wcs.crval[1], min=-90, max=90, vary=True) @@ -176,7 +186,7 @@ def refine_pointing(image, guess_wcs, observed_coords=None, catalog=None, while trial_num < max_trials: try: out = minimize(_residual, params, method=method, - args=(observed_coords, catalog, guess_wcs, num_stars, image.shape)) + args=(observed_coords, catalog, guess_wcs, num_stars, image.shape, edge, sigma, mask)) chisqr = out.chisqr result_minimizations.append(out) except IndexError: @@ -195,6 +205,7 @@ def refine_pointing(image, guess_wcs, observed_coords=None, catalog=None, result_wcs.wcs.pc = calculate_pc_matrix(out.params["crota"], result_wcs.wcs.cdelt) result_wcs.cpdis1 = guess_wcs.cpdis1 # TODO: what if there is no known distortion result_wcs.cpdis2 = guess_wcs.cpdis2 + result_wcs.wcs.set_pv(guess_wcs.wcs.get_pv()) result_wcses.append(result_wcs) trial_num += 1 if chisqr < chisqr_threshold: From 8240ab1c66de1b0439450eb6ad4fbae8f743dbcb Mon Sep 17 00:00:00 2001 From: Marcus Hughes Date: Sat, 24 Aug 2024 12:37:53 -0600 Subject: [PATCH 2/4] Update CITATION.cff --- CITATION.cff | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CITATION.cff b/CITATION.cff index d029ceb..3de1da2 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -5,7 +5,7 @@ authors: given-names: "J. Marcus" orcid: "https://orcid.org/0000-0003-3410-7650" title: "thuban" -version: 0.0.2 +version: 0.0.3 doi: 10.5281/zenodo.13340312 date-released: 2024-08-18 url: "https://github.com/punch-mission/thuban" From 4bc781e0fd36270d28337f57dbfe45e95a86325c Mon Sep 17 00:00:00 2001 From: Marcus Hughes Date: Sat, 24 Aug 2024 12:38:03 -0600 Subject: [PATCH 3/4] Update pyproject.toml --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index b1c0805..b2e7fcf 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,7 +5,7 @@ requires = ["setuptools", [project] name = "thuban" -version = "0.0.2" +version = "0.0.3" dependencies = ["numpy", "astropy", "matplotlib", From 75d0ed92eef8914eb88d464ba9eb6b5fdd614462 Mon Sep 17 00:00:00 2001 From: Marcus Hughes Date: Sat, 24 Aug 2024 12:38:14 -0600 Subject: [PATCH 4/4] Update conf.py --- docs/conf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/conf.py b/docs/conf.py index c564c84..84aa6b9 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -9,7 +9,7 @@ project = "thuban" copyright = "2024, PUNCH Science Operations Center" author = "PUNCH Science Operations Center" -release = "0.0.2" +release = "0.0.3" # -- General configuration --------------------------------------------------- # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration