From a1c62e50f458b7842f6773b1fb350e3557cb00ad Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 19 May 2026 09:33:27 -0400 Subject: [PATCH 1/3] Ruff options and pre-commit. --- .pre-commit-config.yaml | 43 +++++++++++++++++++++++++++++++++++++++++ pyproject.toml | 25 +++++++++++++----------- 2 files changed, 57 insertions(+), 11 deletions(-) create mode 100644 .pre-commit-config.yaml diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 0000000..bff8a2b --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,43 @@ +repos: + - repo: https://github.com/pre-commit/pre-commit-hooks + rev: v6.0.0 + hooks: + - id: check-added-large-files + args: ["--enforce-all", "--maxkb=300"] + exclude: "^(\ + cextern/expat/lib/xmlparse.c|\ + cextern/wcslib/C/flexed/.*|\ + CHANGES.rst|\ + )$" + # Prevent giant files from being committed. + - id: check-case-conflict + # Check for files with names that would conflict on a case-insensitive + # filesystem like MacOS HFS+ or Windows FAT. + - id: check-json + # Attempts to load all json files to verify syntax. + - id: check-merge-conflict + # Check for files that contain merge conflict strings. + - id: check-symlinks + # Checks for symlinks which do not point to anything. + - id: check-toml + # Attempts to load all TOML files to verify syntax. + - id: check-xml + # Attempts to load all xml files to verify syntax. + - id: check-yaml + # Attempts to load all yaml files to verify syntax. + exclude: ".*(.github.*)$" + - id: detect-private-key + # Checks for the existence of private keys. + - id: end-of-file-fixer + # Makes sure files end in a newline and only a newline. + exclude: ".*(data.*|extern.*|licenses.*|_static.*|_parsetab.py)$" + - id: trailing-whitespace + # Trims trailing whitespace. + exclude_types: [python] # Covered by Ruff W291. + exclude: ".*(data.*|extern.*|licenses.*|_static.*)$" + - repo: https://github.com/astral-sh/ruff-pre-commit + rev: v0.13.3 + hooks: + - id: ruff-check + args: ["--fix", "--show-fixes"] + - id: ruff-format diff --git a/pyproject.toml b/pyproject.toml index f7d4b0a..6509d33 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -26,6 +26,13 @@ dependencies = [ ] dynamic = ["version"] +[project.optional-dependencies] +dev = ["build", "pre-commit", "ruff"] +test = ["pytest", "pytest-remotedata", "pytest-cov", "tox"] + +[project.urls] +homepage = "https://catch.astro.umd.edu/" + [build-system] requires = ["setuptools>=80", "setuptools-scm>=8", "wheel"] build-backend = "setuptools.build_meta" @@ -33,23 +40,19 @@ build-backend = "setuptools.build_meta" [tool.setuptools] packages = ["catch_analysis_tools"] -[project.urls] -homepage = "https://catch.astro.umd.edu/" - [tool.setuptools_scm] write_to = "catch_analysis_tools/version.py" -[project.optional-dependencies] -test = [ - "pytest", - "pytest-remotedata", - "pytest-cov", - "tox" +[tool.ruff] +select = [ + "E", # pycodestyle errors + "W", # pycodestyle warnings + "F", # pyflakes. "E" + "W" + "F" + "C90" (mccabe complexity) is equivalent to flake8 + "I", # isort ] -dev = ["build"] [tool.pytest.ini_options] addopts = "--cov=catch_analysis_tools --cov-report=term-missing" markers = [ - "integration: tests that run solve-field and require astrometry.net index files" + "integration: tests that run solve-field and require astrometry.net index files", ] From b5001460931369d257fa4615147955017e28b9b4 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Mon, 15 Jun 2026 08:51:07 -0400 Subject: [PATCH 2/3] Reformat with ruff, initial linting --- catch_analysis_tools/__init__.py | 4 +- catch_analysis_tools/app/app.py | 59 +--- .../download_index_files.py | 12 +- .../get_remote_index_files.py | 12 +- .../prepare_astrometry_data.py | 40 +-- .../start_astrometry_background_check.py | 36 +-- .../app/astrometry_readiness/state.py | 6 +- .../app/handlers/astrometry.py | 7 +- catch_analysis_tools/app/handlers/health.py | 2 +- catch_analysis_tools/app/handlers/reset.py | 2 +- .../app/services/astrometry.py | 241 +++++++++----- .../app/services/photometry.py | 209 ++++++------ .../services/subtract_median_background.py | 20 +- catch_analysis_tools/astrometry.py | 302 +++++++++++------- catch_analysis_tools/background.py | 96 +++--- .../calibration/astrometry.py | 85 ++--- .../calibration/photometry.py | 7 +- catch_analysis_tools/photometry.py | 200 ++++++------ catch_analysis_tools/tests/astrometry_test.py | 143 ++++++--- catch_analysis_tools/tests/background_test.py | 35 +- catch_analysis_tools/tests/photometry_test.py | 87 +++-- pyproject.toml | 2 +- 22 files changed, 913 insertions(+), 694 deletions(-) diff --git a/catch_analysis_tools/__init__.py b/catch_analysis_tools/__init__.py index f4a5185..7335543 100644 --- a/catch_analysis_tools/__init__.py +++ b/catch_analysis_tools/__init__.py @@ -1,4 +1,6 @@ -from importlib.metadata import version as _version, PackageNotFoundError +from importlib.metadata import PackageNotFoundError +from importlib.metadata import version as _version + try: __version__ = _version(__name__) except PackageNotFoundError: diff --git a/catch_analysis_tools/app/app.py b/catch_analysis_tools/app/app.py index 8233fa7..08bf045 100644 --- a/catch_analysis_tools/app/app.py +++ b/catch_analysis_tools/app/app.py @@ -3,19 +3,21 @@ Entry point to Flask-Connexion API """ -import logging - import connexion from connexion.middleware import MiddlewarePosition from starlette.middleware.cors import CORSMiddleware -#from . import __version__ as version -#from .config.logging import get_logger -#from .config.env import ENV -#from .config.exceptions import SBNSISException -#from .services.database_provider import db_session +from .astrometry_readiness.start_astrometry_background_check import ( + start_astrometry_background_check, +) + +# from . import __version__ as version +# from .config.logging import get_logger +# from .config.env import ENV +# from .config.exceptions import SBNSISException +# from .services.database_provider import db_session -#logger: logging.Logger = get_logger() +# logger: logging.Logger = get_logger() app = connexion.FlaskApp(__name__, specification_dir="api/") app.add_middleware( @@ -28,51 +30,16 @@ ) - app.add_api( "openapi.yaml", arguments={ - #"version": str(version), - #"base_href": ENV.BASE_HREF, + # "version": str(version), + # "base_href": ENV.BASE_HREF, }, ) application = app.app -#@application.teardown_appcontext -#def shutdown_db_session(exception: Exception = None) -> None: -# db_session.remove() - - -#@application.errorhandler(SBNSISException) -#def handle_sbnsis_error(error: Exception): -# """Log errors. -# -# The HTTP status code is based on the exception, or 500 if it is not defined. -# -# """ -# -# get_logger().exception("SBS Survey Image Service error.") -# return str(error), getattr(error, "code", 500) - - -#@application.errorhandler(Exception) -#def handle_other_error(error: Exception): -# """Log errors.""" -# -# get_logger().exception("An error occurred.") -# return ( -# "Unexpected error. Please report if the problem persists.", -# getattr(error, "code", 500), -# ) - - if __name__ == "__main__": - # for development - #logger.info("Running " + ENV.APP_NAME) - #logger.info(application.url_map) - from catch_analysis_tools.app.astrometry_readiness.start_astrometry_background_check import ( - start_astrometry_background_check, - ) start_astrometry_background_check() - app.run(host="0.0.0.0", port=8000)#"catch_analysis_tools.app:app")#, host=ENV.API_HOST, port=ENV.API_PORT) + app.run(host="0.0.0.0", port=8000) diff --git a/catch_analysis_tools/app/astrometry_readiness/download_index_files.py b/catch_analysis_tools/app/astrometry_readiness/download_index_files.py index 49a1b12..97cab0f 100644 --- a/catch_analysis_tools/app/astrometry_readiness/download_index_files.py +++ b/catch_analysis_tools/app/astrometry_readiness/download_index_files.py @@ -1,14 +1,10 @@ import os from urllib.parse import urljoin -from catch_analysis_tools.app.astrometry_readiness.constants import INDEX_URL -from catch_analysis_tools.app.astrometry_readiness.count_complete_index_files import ( - count_complete_index_files, -) -from catch_analysis_tools.app.astrometry_readiness.get_index_dir import get_index_dir -from catch_analysis_tools.app.astrometry_readiness.set_astrometry_readiness_status import ( - set_astrometry_readiness_status, -) +from .constants import INDEX_URL +from .count_complete_index_files import count_complete_index_files +from .get_index_dir import get_index_dir +from .set_astrometry_readiness_status import set_astrometry_readiness_status def download_index_files(session, expected_files): diff --git a/catch_analysis_tools/app/astrometry_readiness/get_remote_index_files.py b/catch_analysis_tools/app/astrometry_readiness/get_remote_index_files.py index cc89ec4..39d656e 100644 --- a/catch_analysis_tools/app/astrometry_readiness/get_remote_index_files.py +++ b/catch_analysis_tools/app/astrometry_readiness/get_remote_index_files.py @@ -8,11 +8,13 @@ def get_remote_index_files(session): response = session.get(INDEX_URL, timeout=60) response.raise_for_status() soup = BeautifulSoup(response.text, "html.parser") - files = sorted({ - link["href"] - for link in soup.find_all("a", href=True) - if link["href"].startswith("index-") and link["href"].endswith(".fits") - }) + files = sorted( + { + link["href"] + for link in soup.find_all("a", href=True) + if link["href"].startswith("index-") and link["href"].endswith(".fits") + } + ) if not files: raise RuntimeError(f"No astrometry index files found at {INDEX_URL}") return files diff --git a/catch_analysis_tools/app/astrometry_readiness/prepare_astrometry_data.py b/catch_analysis_tools/app/astrometry_readiness/prepare_astrometry_data.py index 073c247..35cc0ef 100644 --- a/catch_analysis_tools/app/astrometry_readiness/prepare_astrometry_data.py +++ b/catch_analysis_tools/app/astrometry_readiness/prepare_astrometry_data.py @@ -2,32 +2,16 @@ import requests -from catch_analysis_tools.app.astrometry_readiness.acquire_index_download_lock import ( - acquire_index_download_lock, -) -from catch_analysis_tools.app.astrometry_readiness.constants import INDEX_URL -from catch_analysis_tools.app.astrometry_readiness.count_complete_index_files import ( - count_complete_index_files, -) -from catch_analysis_tools.app.astrometry_readiness.download_index_files import ( - download_index_files, -) -from catch_analysis_tools.app.astrometry_readiness.get_astrometry_readiness_status import ( - get_astrometry_readiness_status, -) -from catch_analysis_tools.app.astrometry_readiness.get_index_dir import get_index_dir -from catch_analysis_tools.app.astrometry_readiness.get_remote_index_files import ( - get_remote_index_files, -) -from catch_analysis_tools.app.astrometry_readiness.index_files_complete import ( - index_files_complete, -) -from catch_analysis_tools.app.astrometry_readiness.set_astrometry_readiness_status import ( - set_astrometry_readiness_status, -) -from catch_analysis_tools.app.astrometry_readiness.write_ready_sentinel import ( - write_ready_sentinel, -) +from .acquire_index_download_lock import acquire_index_download_lock +from .constants import INDEX_URL +from .count_complete_index_files import count_complete_index_files +from .download_index_files import download_index_files +from .get_astrometry_readiness_status import get_astrometry_readiness_status +from .get_index_dir import get_index_dir +from .get_remote_index_files import get_remote_index_files +from .index_files_complete import index_files_complete +from .set_astrometry_readiness_status import set_astrometry_readiness_status +from .write_ready_sentinel import write_ready_sentinel def prepare_astrometry_data(force=False): @@ -78,7 +62,9 @@ def prepare_astrometry_data(force=False): set_astrometry_readiness_status( state="downloading", ready=False, - message="Astrometry index files are incomplete. Downloading missing files.", + message=( + "Astrometry index files are incomplete. Downloading missing files." + ), files_present=files_present, error=None, ) diff --git a/catch_analysis_tools/app/astrometry_readiness/start_astrometry_background_check.py b/catch_analysis_tools/app/astrometry_readiness/start_astrometry_background_check.py index c6d52be..d291cc7 100644 --- a/catch_analysis_tools/app/astrometry_readiness/start_astrometry_background_check.py +++ b/catch_analysis_tools/app/astrometry_readiness/start_astrometry_background_check.py @@ -1,14 +1,10 @@ import threading -from catch_analysis_tools.app.astrometry_readiness import state -from catch_analysis_tools.app.astrometry_readiness.background_prepare_astrometry_data import ( - background_prepare_astrometry_data, -) -from catch_analysis_tools.app.astrometry_readiness.constants import INDEX_URL -from catch_analysis_tools.app.astrometry_readiness.get_current_time import ( - get_current_time, -) -from catch_analysis_tools.app.astrometry_readiness.get_index_dir import get_index_dir +from . import state +from .background_prepare_astrometry_data import background_prepare_astrometry_data +from .constants import INDEX_URL +from .get_current_time import get_current_time +from .get_index_dir import get_index_dir def start_astrometry_background_check(force=False): @@ -17,16 +13,18 @@ def start_astrometry_background_check(force=False): if state.worker is not None and state.worker.is_alive(): return dict(state.status) - state.status.update({ - "state": "checking", - "ready": False, - "message": "Astrometry data readiness check has started.", - "index_dir": str(get_index_dir()), - "index_url": INDEX_URL, - "expected_files": None, - "error": None, - "updated_at": get_current_time(), - }) + state.status.update( + { + "state": "checking", + "ready": False, + "message": "Astrometry data readiness check has started.", + "index_dir": str(get_index_dir()), + "index_url": INDEX_URL, + "expected_files": None, + "error": None, + "updated_at": get_current_time(), + } + ) state.worker = threading.Thread( target=background_prepare_astrometry_data, kwargs={"force": force}, diff --git a/catch_analysis_tools/app/astrometry_readiness/state.py b/catch_analysis_tools/app/astrometry_readiness/state.py index ff5026e..050667f 100644 --- a/catch_analysis_tools/app/astrometry_readiness/state.py +++ b/catch_analysis_tools/app/astrometry_readiness/state.py @@ -1,10 +1,6 @@ import threading -from catch_analysis_tools.app.astrometry_readiness.constants import ( - DEFAULT_INDEX_DIR, - INDEX_URL, -) - +from .constants import DEFAULT_INDEX_DIR, INDEX_URL state_lock = threading.RLock() worker = None diff --git a/catch_analysis_tools/app/handlers/astrometry.py b/catch_analysis_tools/app/handlers/astrometry.py index 41924de..0d03f2c 100644 --- a/catch_analysis_tools/app/handlers/astrometry.py +++ b/catch_analysis_tools/app/handlers/astrometry.py @@ -9,20 +9,19 @@ from flask import Response from werkzeug.exceptions import BadRequest -from catch_analysis_tools.app.astrometry_readiness.get_astrometry_readiness_status import ( +from ..astrometry_readiness.get_astrometry_readiness_status import ( get_astrometry_readiness_status, ) -from catch_analysis_tools.app.astrometry_readiness.is_astrometry_ready import ( +from ..astrometry_readiness.is_astrometry_ready import ( is_astrometry_ready, ) -from catch_analysis_tools.app.services.astrometry import ( +from ..services.astrometry import ( AstrometrySolveError, AstrometryValidationError, run_pipeline, validate_and_normalize, ) - logger = logging.getLogger(__name__) diff --git a/catch_analysis_tools/app/handlers/health.py b/catch_analysis_tools/app/handlers/health.py index 654e599..035f0b3 100644 --- a/catch_analysis_tools/app/handlers/health.py +++ b/catch_analysis_tools/app/handlers/health.py @@ -1,4 +1,4 @@ -from catch_analysis_tools.app.astrometry_readiness.get_astrometry_readiness_status import ( +from ..astrometry_readiness.get_astrometry_readiness_status import ( get_astrometry_readiness_status, ) diff --git a/catch_analysis_tools/app/handlers/reset.py b/catch_analysis_tools/app/handlers/reset.py index 0ee0d75..6e68e28 100644 --- a/catch_analysis_tools/app/handlers/reset.py +++ b/catch_analysis_tools/app/handlers/reset.py @@ -1,4 +1,4 @@ -from catch_analysis_tools.app.astrometry_readiness.start_astrometry_background_check import ( +from ..astrometry_readiness.start_astrometry_background_check import ( start_astrometry_background_check, ) diff --git a/catch_analysis_tools/app/services/astrometry.py b/catch_analysis_tools/app/services/astrometry.py index 5d826d4..1f4bcd3 100644 --- a/catch_analysis_tools/app/services/astrometry.py +++ b/catch_analysis_tools/app/services/astrometry.py @@ -1,21 +1,23 @@ -import numpy as np +import base64 +import io import os import subprocess -import pandas as pd -import sep +from copy import deepcopy +from typing import Any + +import calviacat as cvc import fitsio import matplotlib -matplotlib.use("Agg") -import io -import base64 import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import sep +from astropy.coordinates import SkyCoord from astropy.io import fits -from astropy.wcs import WCS from astropy.table import Table -from astropy.coordinates import SkyCoord -from typing import Dict, Any -from copy import deepcopy -import calviacat as cvc +from astropy.wcs import WCS + +matplotlib.use("Agg") class AstrometrySolveError(RuntimeError): @@ -46,7 +48,7 @@ class AstrometryValidationError(ValueError): "output": { "make_plots": False, "write_fits": False, - } + }, } @@ -66,7 +68,7 @@ def run_solve_field(input_fits, output_wcs, wcs_cfg): """ if os.path.exists(output_wcs): return True - + if wcs_cfg["use_ra_dec"]: if wcs_cfg.get("ra") is None or wcs_cfg.get("dec") is None: raise RuntimeError("RA/Dec missing while use_ra_dec=True") @@ -78,39 +80,47 @@ def run_solve_field(input_fits, output_wcs, wcs_cfg): config_file = os.environ.get("ASTROMETRY_CONFIG") if config_file is None: raise RuntimeError( - "ASTROMETRY_CONFIG is not set. " - "This is required to run solve-field." + "ASTROMETRY_CONFIG is not set. This is required to run solve-field." ) command = [ "solve-field", "--overwrite", - "--config", config_file, + "--config", + config_file, "--fits-image", - "--wcs", output_wcs, - "--no-plots", + "--wcs", + output_wcs, + "--no-plots", ] # --- RA/Dec only if enabled --- if wcs_cfg["use_ra_dec"]: command += [ - "--ra", str(wcs_cfg["ra"]), - "--dec", str(wcs_cfg["dec"]), + "--ra", + str(wcs_cfg["ra"]), + "--dec", + str(wcs_cfg["dec"]), ] command += [ - "--scale-units", "arcsecperpix", - "--scale-low", str(scale_low), - "--scale-high", str(scale_high), + "--scale-units", + "arcsecperpix", + "--scale-low", + str(scale_low), + "--scale-high", + str(scale_high), ] if wcs_cfg["use_ra_dec"]: command += [ - "--radius", str(int(wcs_cfg["search_radius"])), + "--radius", + str(int(wcs_cfg["search_radius"])), ] command += [ - "--downsample", "1", + "--downsample", + "1", input_fits, ] @@ -150,26 +160,23 @@ def find_sources(image, det_cfg): sep.set_sub_object_limit(500) sources = sep.extract( - image_sub, - thresh=det_cfg["snr"], - err=bkg.globalrms, - deblend_nthresh=16 + image_sub, thresh=det_cfg["snr"], err=bkg.globalrms, deblend_nthresh=16 ) source_list = pd.DataFrame(sources) flux, flux_err, _ = sep.sum_circle( image_sub, - source_list['x'], - source_list['y'], + source_list["x"], + source_list["y"], det_cfg["aperture_radius"], - err=bkg.globalrms + err=bkg.globalrms, ) - source_list['aperture_sum'] = flux - source_list['aperture_err'] = flux_err + source_list["aperture_sum"] = flux + source_list["aperture_err"] = flux_err - source_list = source_list[source_list['aperture_sum'] > 0].reset_index(drop=True) + source_list = source_list[source_list["aperture_sum"] > 0].reset_index(drop=True) return source_list, image_sub @@ -213,10 +220,10 @@ def retrieve_sources(source_list, wcs_solution): sky_coords : astropy.coordinates.SkyCoord SkyCoord object with celestial coordinates of sources. """ - world = wcs_solution.pixel_to_world(source_list['x'], source_list['y']) - source_list['RA'] = [c.ra.deg for c in world] - source_list['Dec'] = [c.dec.deg for c in world] - sky_coords = SkyCoord(source_list['RA'], source_list['Dec'], unit='deg') + world = wcs_solution.pixel_to_world(source_list["x"], source_list["y"]) + source_list["RA"] = [c.ra.deg for c in world] + source_list["Dec"] = [c.dec.deg for c in world] + sky_coords = SkyCoord(source_list["RA"], source_list["Dec"], unit="deg") return source_list, sky_coords @@ -268,7 +275,7 @@ def calibrate_photometry(sky_coords, source_list, phot_cfg): objids, distances = ref.xmatch(sky_coords) - m_inst = -2.5 * np.log10(source_list['aperture_sum'].values) + m_inst = -2.5 * np.log10(source_list["aperture_sum"].values) zp, C, unc, m_cal, color_mags, _ = ref.cal_color( objids, @@ -292,14 +299,7 @@ def calibrate_photometry(sky_coords, source_list, phot_cfg): } -def plot_color_correction( - color_mags, - m, - m_inst, - C, - zp, - color_index: str -): +def plot_color_correction(color_mags, m, m_inst, C, zp, color_index: str): """ Plot the relation between instrumental and calibrated magnitudes. @@ -324,11 +324,11 @@ def plot_color_correction( Matplotlib figure and axis objects for the plot. """ fig, ax = plt.subplots() - ax.scatter(color_mags, m - m_inst, marker='.') + ax.scatter(color_mags, m - m_inst, marker=".") x = np.linspace(0, 1.5, 100) - ax.plot(x, C * x + zp, color='red', label=f'$m = C\\times({color_index}) + ZP$') - ax.set_xlabel(f'${color_index}$ (mag)') - ax.set_ylabel(r'$m - m_{\mathrm{inst}}$ (mag)') + ax.plot(x, C * x + zp, color="red", label=f"$m = C\\times({color_index}) + ZP$") + ax.set_xlabel(f"${color_index}$ (mag)") + ax.set_ylabel(r"$m - m_{\mathrm{inst}}$ (mag)") plt.tight_layout() return fig, ax @@ -355,18 +355,55 @@ def plot_image(telescope_image_sub, source_list, matched_idx, colored_idx): """ fig, ax = plt.subplots() m, s = np.mean(telescope_image_sub), np.std(telescope_image_sub) - im = ax.imshow(telescope_image_sub, interpolation='nearest', origin='lower', cmap='gray') - im.set_clim(vmin=m-s, vmax=m+s) + im = ax.imshow( + telescope_image_sub, interpolation="nearest", origin="lower", cmap="gray" + ) + im.set_clim(vmin=m - s, vmax=m + s) fig.colorbar(im, ax=ax) - ax.plot(source_list['x'], source_list['y'], '+', markersize=5, label='Detected', color='red',) - ax.plot(source_list['x'].iloc[matched_idx], source_list['y'].iloc[matched_idx], 'o', markersize=10, color='blue', markerfacecolor='none', label='Matched') - ax.plot(source_list['x'].iloc[colored_idx], source_list['y'].iloc[colored_idx], 'o', markersize=15, color='yellow', markerfacecolor='none', label='Selected for Color Corr') + ax.plot( + source_list["x"], + source_list["y"], + "+", + markersize=5, + label="Detected", + color="red", + ) + ax.plot( + source_list["x"].iloc[matched_idx], + source_list["y"].iloc[matched_idx], + "o", + markersize=10, + color="blue", + markerfacecolor="none", + label="Matched", + ) + ax.plot( + source_list["x"].iloc[colored_idx], + source_list["y"].iloc[colored_idx], + "o", + markersize=15, + color="yellow", + markerfacecolor="none", + label="Selected for Color Corr", + ) ax.legend() return fig, ax -def create_header(image, wcs_solution, zp, unc, source_list, matched_idx, colored_idx, input_fits, cal_band: str, catalog: str, obj_band: str): +def create_header( + image, + wcs_solution, + zp, + unc, + source_list, + matched_idx, + colored_idx, + input_fits, + cal_band: str, + catalog: str, + obj_band: str, +): """ Create and write a FITS file with calibrated header and source tables. @@ -395,20 +432,34 @@ def create_header(image, wcs_solution, zp, unc, source_list, matched_idx, colore """ image_arr = np.asarray(image) primary_hdu = fits.PrimaryHDU(data=image_arr, header=wcs_solution.to_header()) - primary_hdu.header['ZP'] = zp - primary_hdu.header['ZP_STD'] = unc - primary_hdu.header['SUV_FLT'] = cal_band - primary_hdu.header['REF_CATA'] = catalog - primary_hdu.header['REF_FLT'] = obj_band - primary_hdu.header['CAT_COR'] = f"{cal_band}-{obj_band}" - source_list_clean = source_list.map(lambda x: x.filled(np.nan) if hasattr(x, 'filled') else x) - detected_hdu = fits.BinTableHDU(Table.from_pandas(source_list_clean), name='DETECTED_SOURCES') + primary_hdu.header["ZP"] = zp + primary_hdu.header["ZP_STD"] = unc + primary_hdu.header["SUV_FLT"] = cal_band + primary_hdu.header["REF_CATA"] = catalog + primary_hdu.header["REF_FLT"] = obj_band + primary_hdu.header["CAT_COR"] = f"{cal_band}-{obj_band}" + source_list_clean = source_list.map( + lambda x: x.filled(np.nan) if hasattr(x, "filled") else x + ) + detected_hdu = fits.BinTableHDU( + Table.from_pandas(source_list_clean), name="DETECTED_SOURCES" + ) if not source_list_clean.empty: - matched_hdu = fits.BinTableHDU(Table.from_pandas(source_list_clean.iloc[matched_idx].reset_index(drop=True)), name='SELECTED_STARS') - colored_hdu = fits.BinTableHDU(Table.from_pandas(source_list_clean.iloc[colored_idx].reset_index(drop=True)), name='START_COLOR_CORRECTION') + matched_hdu = fits.BinTableHDU( + Table.from_pandas( + source_list_clean.iloc[matched_idx].reset_index(drop=True) + ), + name="SELECTED_STARS", + ) + colored_hdu = fits.BinTableHDU( + Table.from_pandas( + source_list_clean.iloc[colored_idx].reset_index(drop=True) + ), + name="START_COLOR_CORRECTION", + ) else: - matched_hdu = fits.BinTableHDU(name='SELECTED_STARS') - colored_hdu = fits.BinTableHDU(name='START_COLOR_CORRECTION') + matched_hdu = fits.BinTableHDU(name="SELECTED_STARS") + colored_hdu = fits.BinTableHDU(name="START_COLOR_CORRECTION") hdul = fits.HDUList([primary_hdu, detected_hdu, matched_hdu, colored_hdu]) hdul.writeto(input_fits, overwrite=True) @@ -426,8 +477,18 @@ def cleanup_files(file_base): ------- None """ - extensions = ['.axy', '.corr', '.match', '.new', '.rdls', '.solved', - '-ngc.png', '-objs.png', '-indx.png', '-indx.xyls'] + extensions = [ + ".axy", + ".corr", + ".match", + ".new", + ".rdls", + ".solved", + "-ngc.png", + "-objs.png", + "-indx.png", + "-indx.xyls", + ] for ext in extensions: fname = f"{file_base}{ext}" if os.path.exists(fname): @@ -435,7 +496,6 @@ def cleanup_files(file_base): def run_pipeline(input_fits: str, user_config: dict) -> dict: - config = merge_config(user_config, DEFAULT_CONFIG) wcs_cfg = config["wcs"] @@ -475,8 +535,16 @@ def run_pipeline(input_fits: str, user_config: dict) -> dict: objids = calibration["objids"] # --- Matching indices --- - matched_idx = np.where(~objids.mask)[0] if hasattr(objids, "mask") else np.arange(len(source_list)) - colored_idx = np.where(~color_mags.mask)[0] if hasattr(color_mags, "mask") else np.arange(len(source_list)) + matched_idx = ( + np.where(~objids.mask)[0] + if hasattr(objids, "mask") + else np.arange(len(source_list)) + ) + colored_idx = ( + np.where(~color_mags.mask)[0] + if hasattr(color_mags, "mask") + else np.arange(len(source_list)) + ) # --- Plotting --- plots = {} @@ -496,8 +564,13 @@ def run_pipeline(input_fits: str, user_config: dict) -> dict: # --- Optional FITS --- if out_cfg["write_fits"]: create_header( - image, wcs_solution, zp, unc, - source_list, matched_idx, colored_idx, + image, + wcs_solution, + zp, + unc, + source_list, + matched_idx, + colored_idx, input_fits, phot_cfg["obs_band"], phot_cfg["catalog"], @@ -520,7 +593,7 @@ def run_pipeline(input_fits: str, user_config: dict) -> dict: "center_ra_deg": center_ra, "center_dec_deg": center_dec, "pixel_scale": wcs_cfg["pixel_scale"], - } + }, } if out_cfg["make_plots"]: @@ -529,7 +602,7 @@ def run_pipeline(input_fits: str, user_config: dict) -> dict: return results -def validate_and_normalize(body: Dict[str, Any]) -> Dict[str, Any]: +def validate_and_normalize(body: dict[str, Any]) -> dict[str, Any]: """ Minimal validation + normalization layer. Keeps API robust without adding heavy dependencies. @@ -614,29 +687,24 @@ def get_bool(name, default=False): "meta": { "return_plot": return_plot, "plot_type": plot_type, - } + }, } # ## when testing the code locally if __name__ == "__main__": - # --- Simulated API request body --- body = { "image_url": "Comet_65P_Gunn_LONEOS.fits", # not used in local test "ra": 51.0, "dec": 17.0, "use_ra_dec": True, - "pixel_scale": 2.5, - "snr_threshold": 3.0, "aperture_radius": 7.0, - "catalog": "PanSTARRS1", "obs_band": "g", "cal_band": "r", - "return_plot": True, "plot_type": "color_correction", } @@ -660,9 +728,10 @@ def get_bool(name, default=False): # --- Optional: visualize plot locally --- if cfg["output"]["make_plots"]: import base64 - import matplotlib.pyplot as plt import io + import matplotlib.pyplot as plt + plot_type = cfg["meta"]["plot_type"] img_bytes = base64.b64decode(results["plots"][plot_type]) diff --git a/catch_analysis_tools/app/services/photometry.py b/catch_analysis_tools/app/services/photometry.py index 28905dd..1245266 100644 --- a/catch_analysis_tools/app/services/photometry.py +++ b/catch_analysis_tools/app/services/photometry.py @@ -1,73 +1,76 @@ -import os -from ...photometry import get_image,subpixel_centroid,define_aperture,do_aperture_photometry -from astropy.wcs import WCS import matplotlib -matplotlib.use('Agg') -from matplotlib import pyplot as plt -from astropy.visualization import ZScaleInterval -from astropy.visualization.mpl_normalize import ImageNormalize +from astropy.wcs import WCS + +from ...photometry import ( + define_aperture, + do_aperture_photometry, + get_image, + subpixel_centroid, +) + +matplotlib.use("Agg") from astropy import units as u -from astropy.coordinates import SkyCoord, ICRS +from astropy.coordinates import ICRS, SkyCoord from astropy.io import fits +from astropy.visualization import ZScaleInterval +from astropy.visualization.mpl_normalize import ImageNormalize +from matplotlib import pyplot as plt + -def get_world_coordinates(WCS_file,x,y): +def get_world_coordinates(WCS_file, x, y): """ Accepts an astropy-readable WCS object (a .wcs or fits file ) and outputs the world coordinates of a 0-indexed (x,y) pixel point in decimal degrees. """ try: world_coords = WCS(fits.open(WCS_file)[0].header) - except Exception as e: + except Exception: try: world_coords = WCS(WCS_file) - except Exception as e: + except Exception: try: world_coords = WCS_file - except Exception as e: - raise ValueError("Could not read WCS file. Please ensure that the file is either a .wcs file or a .fits file with a valid WCS solution in the header.") - - loc = world_coords.pixel_to_world(x,y) + except Exception: + raise ValueError( + "Could not read WCS file. Please ensure that the file is either a .wcs file or a .fits file with a valid WCS solution in the header." + ) + + loc = world_coords.pixel_to_world(x, y) print(loc) - transform_results = { - "x":x, - "y":y, - "ra":loc.ra.deg, - "dec":loc.dec.deg - } + transform_results = {"x": x, "y": y, "ra": loc.ra.deg, "dec": loc.dec.deg} return transform_results -def get_pixel_coordinates(WCS_file,ra,dec): + +def get_pixel_coordinates(WCS_file, ra, dec): """ Accepts an astropy-readable WCS object (a .wcs or fits file) and outputs the 0-indexed (x,y) pixel coordinates of an (ra,dec) point in decimal degrees. """ try: world_coords = WCS(fits.open(WCS_file)[0].header) - except Exception as e: + except Exception: try: world_coords = WCS(WCS_file) - except Exception as e: + except Exception: try: world_coords = WCS_file - except Exception as e: - raise ValueError("Could not read WCS file. Please ensure that the file is either a .wcs file or a .fits file with a valid WCS solution in the header.") - - sky_loc = SkyCoord(ICRS(ra=ra*u.deg, dec=dec*u.deg)) + except Exception: + raise ValueError( + "Could not read WCS file. Please ensure that the file is either a .wcs file or a .fits file with a valid WCS solution in the header." + ) + + sky_loc = SkyCoord(ICRS(ra=ra * u.deg, dec=dec * u.deg)) loc = world_coords.world_to_pixel(sky_loc) x = loc[0].item() y = loc[1].item() - transform_results = { - "x":x, - "y":y, - "ra":ra, - "dec":dec - } - return transform_results + transform_results = {"x": x, "y": y, "ra": ra, "dec": dec} + return transform_results + -def centroid(file,target_x,target_y,search_radius): +def centroid(file, target_x, target_y, search_radius): """ - Searches for source nearby to expected ephemeris location in image, assumes that astrometry solution has been rerun and that both the cutout .fits file and the redone .wcs file exist. - - + Searches for source nearby to expected ephemeris location in image, assumes that astrometry solution has been rerun and that both the cutout .fits file and the redone .wcs file exist. + + Parameters ---------- file : string @@ -87,102 +90,118 @@ def centroid(file,target_x,target_y,search_radius): ------- search_results : array_like - - - figname: + + + figname: string Output plot showing the default aperture + annulus extraction onto the cutout image. """ img, header = get_image(file) - cent_pix = subpixel_centroid([target_x,target_y],img,search_radius) - + cent_pix = subpixel_centroid([target_x, target_y], img, search_radius) + search_results = { - "init_guess_x":target_x, - "init_guess_y":target_y, - "search_radius":search_radius, - "cent_x":cent_pix[0], - "cent_y":cent_pix[1], + "init_guess_x": target_x, + "init_guess_y": target_y, + "search_radius": search_radius, + "cent_x": cent_pix[0], + "cent_y": cent_pix[1], } - + norm = ImageNormalize(img, interval=ZScaleInterval()) - - plt.figure(figsize=(8,8)) - plt.imshow(img,norm=norm,cmap='gray_r') - plt.scatter(img.shape[0]/2,img.shape[1]/2,s=100,marker='x',label='Image Center') - plt.scatter(target_x,target_y,s=100,marker='+',label='Initial Guess') - plt.scatter(cent_pix[0],cent_pix[1],s=50,marker='.',c='yellow',label='Centroid Location') - plt.xlim(target_x-20,target_x+20) - plt.ylim(target_y-20,target_y+20) + plt.figure(figsize=(8, 8)) + plt.imshow(img, norm=norm, cmap="gray_r") + plt.scatter( + img.shape[0] / 2, img.shape[1] / 2, s=100, marker="x", label="Image Center" + ) + plt.scatter(target_x, target_y, s=100, marker="+", label="Initial Guess") + plt.scatter( + cent_pix[0], + cent_pix[1], + s=50, + marker=".", + c="yellow", + label="Centroid Location", + ) + plt.xlim(target_x - 20, target_x + 20) + plt.ylim(target_y - 20, target_y + 20) plt.legend() - - - - figname = 'centroid.png' + figname = "centroid.png" plt.savefig(figname) plt.close() - centroid_results = { - "search_results": search_results, - "centroid_figure": figname - } + centroid_results = {"search_results": search_results, "centroid_figure": figname} return centroid_results -#todo: create function that takes target location, background location and outputs aperture objects + +# todo: create function that takes target location, background location and outputs aperture objects # this function can be fed the outputs of centroid_location + def target_extraction(body): """ - Searches for source nearby to expected ephemeris location in image, assumes that astrometry solution has been rerun and that both the cutout .fits file and the redone .wcs file exist. - - + Searches for source nearby to expected ephemeris location in image, assumes that astrometry solution has been rerun and that both the cutout .fits file and the redone .wcs file exist. + + Parameters ---------- body : dict Request body containing file, target_aperture_params, and background_aperture_params. """ - file = body['file'] - target_aperture_params = body['target_aperture_params'] - background_aperture_params = body['background_aperture_params'] + file = body["file"] + target_aperture_params = body["target_aperture_params"] + background_aperture_params = body["background_aperture_params"] img, header = get_image(file) - - - - target_aperture = define_aperture(target_aperture_params) background_aperture = define_aperture(background_aperture_params) norm = ImageNormalize(img, interval=ZScaleInterval()) - - target_flux, target_fluxerr = do_aperture_photometry(img,target_aperture,background_aperture) - + + target_flux, target_fluxerr = do_aperture_photometry( + img, target_aperture, background_aperture + ) + # could put code to filter out frames where targ_loc and targ_cent vary by more than a couple pixels here # (would mean star hit) - - plt.figure(figsize=(8,8)) - plt.imshow(img,norm=norm,cmap='gray_r') - target_aperture.plot(color='blue', lw=1.5, alpha=0.5) - background_aperture.plot(color='yellow',lw=1.5) - plt.scatter(target_aperture_params["position"][0],target_aperture_params["position"][1],s=100,marker='+') - plt.scatter(background_aperture_params["position"][0],background_aperture_params["position"][1],s=50,marker='.',c='yellow') - plt.xlim(target_aperture_params["position"][0]-20,target_aperture_params["position"][0]+20) - plt.ylim(target_aperture_params["position"][1]-20,target_aperture_params["position"][1]+20) - - - - figname = 'aperture_extraction.png' + + plt.figure(figsize=(8, 8)) + plt.imshow(img, norm=norm, cmap="gray_r") + target_aperture.plot(color="blue", lw=1.5, alpha=0.5) + background_aperture.plot(color="yellow", lw=1.5) + plt.scatter( + target_aperture_params["position"][0], + target_aperture_params["position"][1], + s=100, + marker="+", + ) + plt.scatter( + background_aperture_params["position"][0], + background_aperture_params["position"][1], + s=50, + marker=".", + c="yellow", + ) + plt.xlim( + target_aperture_params["position"][0] - 20, + target_aperture_params["position"][0] + 20, + ) + plt.ylim( + target_aperture_params["position"][1] - 20, + target_aperture_params["position"][1] + 20, + ) + + figname = "aperture_extraction.png" plt.savefig(figname) plt.close() extraction_results = { "aperture_flux": target_flux, "aperture_fluxerr": target_fluxerr, - "aperture_figure": figname + "aperture_figure": figname, } return extraction_results - diff --git a/catch_analysis_tools/app/services/subtract_median_background.py b/catch_analysis_tools/app/services/subtract_median_background.py index 801420f..811b6e6 100644 --- a/catch_analysis_tools/app/services/subtract_median_background.py +++ b/catch_analysis_tools/app/services/subtract_median_background.py @@ -1,9 +1,12 @@ import os + import requests from astropy.io import fits + from ...background import global_subtraction from ...photometry import get_image + def download_file(url): """ Downloads a file from a given URL and saves it locally, returning the base filename (without extension). @@ -23,13 +26,14 @@ def download_file(url): else: filename = url.split("/")[-1] cleaned_filename = filename.replace(" ", "") - cleaned_filename = cleaned_filename.replace("\"", "") + cleaned_filename = cleaned_filename.replace('"', "") with open(cleaned_filename, mode="wb") as file: file.write(response.content) print(cleaned_filename) print(f"Downloaded file {filename} as {cleaned_filename}") - - return os.path.splitext(cleaned_filename)[0] + + return os.path.splitext(cleaned_filename)[0] + def perform_median_subtraction(url): """Takes a cutout URL from CATCH, opens the file, performs a median background subtraction @@ -49,15 +53,15 @@ def perform_median_subtraction(url): # Download file and get file base name file_base = download_file(url) # Get image data - data, header = get_image(file_base+'.fits') + data, header = get_image(file_base + ".fits") # Perform global background subtraction data_sub, bkg = global_subtraction(data) - header['BKG_MED'] = bkg.background_median - + header["BKG_MED"] = bkg.background_median + # Save background subtracted image - subtract_fname = file_base+'_subtracted.fits' - + subtract_fname = file_base + "_subtracted.fits" + fits.writeto(subtract_fname, data_sub, header, overwrite=True) return subtract_fname diff --git a/catch_analysis_tools/astrometry.py b/catch_analysis_tools/astrometry.py index 706821e..3f1dc00 100644 --- a/catch_analysis_tools/astrometry.py +++ b/catch_analysis_tools/astrometry.py @@ -1,19 +1,22 @@ +import argparse import os import subprocess -import argparse + +import calviacat as cvc +import fitsio +import matplotlib.pyplot as plt import numpy as np import pandas as pd import sep -import fitsio -import matplotlib.pyplot as plt +from astropy.coordinates import SkyCoord from astropy.io import fits -from astropy.wcs import WCS from astropy.table import Table -from astropy.coordinates import SkyCoord -import calviacat as cvc +from astropy.wcs import WCS -def run_solve_field(input_fits, output_wcs, pixel_scale, Ra_deg, Dec_deg, scale_units="arcsecperpix"): +def run_solve_field( + input_fits, output_wcs, pixel_scale, Ra_deg, Dec_deg, scale_units="arcsecperpix" +): """ Execute the `solve-field` command to compute a WCS solution. @@ -35,27 +38,35 @@ def run_solve_field(input_fits, output_wcs, pixel_scale, Ra_deg, Dec_deg, scale_ """ if os.path.exists(output_wcs): print( - f"Output file '{output_wcs}' already exists. Skipping solve-field execution.") + f"Output file '{output_wcs}' already exists. Skipping solve-field execution." + ) return True config_file = os.environ.get("ASTROMETRY_CONFIG") if config_file is None: raise RuntimeError( - "ASTROMETRY_CONFIG is not set. " - "This is required to run solve-field." + "ASTROMETRY_CONFIG is not set. This is required to run solve-field." ) command = [ "solve-field", "--overwrite", - "--config", config_file, - "--ra", str(Ra_deg), - "--dec", str(Dec_deg), - "--scale-units", scale_units, - "--scale-low", str(pixel_scale * 0.5), - "--scale-high", str(pixel_scale * 2.0), - "--radius", "2", - "--downsample", "1", + "--config", + config_file, + "--ra", + str(Ra_deg), + "--dec", + str(Dec_deg), + "--scale-units", + scale_units, + "--scale-low", + str(pixel_scale * 0.5), + "--scale-high", + str(pixel_scale * 2.0), + "--radius", + "2", + "--downsample", + "1", input_fits, ] @@ -89,23 +100,14 @@ def find_sources(image_sub, bkg_err, snr, aperture_radius=7.0): Background-subtracted image array. """ sep.set_sub_object_limit(500) - sources = sep.extract( - image_sub, - thresh=snr, - err=bkg_err, - deblend_nthresh=16 - ) + sources = sep.extract(image_sub, thresh=snr, err=bkg_err, deblend_nthresh=16) source_list = pd.DataFrame(sources) flux, flux_err, _ = sep.sum_circle( - image_sub, - source_list['x'], source_list['y'], - aperture_radius, - err=bkg_err + image_sub, source_list["x"], source_list["y"], aperture_radius, err=bkg_err ) - source_list['aperture_sum'] = flux - source_list['aperture_err'] = flux_err - source_list = source_list[source_list['aperture_sum'] > 0].reset_index( - drop=True) + source_list["aperture_sum"] = flux + source_list["aperture_err"] = flux_err + source_list = source_list[source_list["aperture_sum"] > 0].reset_index(drop=True) return source_list, image_sub @@ -148,19 +150,20 @@ def retrieve_sources(source_list, wcs_solution): sky_coords : astropy.coordinates.SkyCoord SkyCoord object with celestial coordinates of sources. """ - world = wcs_solution.pixel_to_world(source_list['x'], source_list['y']) - source_list['RA'] = [c.ra.deg for c in world] - source_list['Dec'] = [c.dec.deg for c in world] - sky_coords = SkyCoord(source_list['RA'], source_list['Dec'], unit='deg') + world = wcs_solution.pixel_to_world(source_list["x"], source_list["y"]) + source_list["RA"] = [c.ra.deg for c in world] + source_list["Dec"] = [c.dec.deg for c in world] + sky_coords = SkyCoord(source_list["RA"], source_list["Dec"], unit="deg") return source_list, sky_coords def calibrate_photometry( - sky_coords, - source_list, - catalog: str = 'PanSTARRS1', - obs_band: str = 'obs_band', - cal_band: str = 'g'): + sky_coords, + source_list, + catalog: str = "PanSTARRS1", + obs_band: str = "obs_band", + cal_band: str = "g", +): """ Calibrate instrumental magnitudes against a Pan-STARRS1 catalog. @@ -200,12 +203,12 @@ def calibrate_photometry( except AttributeError: raise ValueError(f"Catalog '{catalog}' not found in calviacat") - ref = CatalogClass('cat.db') + ref = CatalogClass("cat.db") results = ref.search(sky_coords) if len(results[0]) < 500: ref.fetch_field(sky_coords) objids, distances = ref.xmatch(sky_coords) - m_inst = -2.5 * np.log10(source_list['aperture_sum'].values) + m_inst = -2.5 * np.log10(source_list["aperture_sum"].values) zp, C, unc, m_cal, color_mags, _ = ref.cal_color( objids, m_inst, @@ -213,28 +216,21 @@ def calibrate_photometry( color_index, ) return { - 'zp': zp, - 'C': C, - 'unc': unc, - 'm': m_cal, - 'm_inst': m_inst, - 'obs_band': obs_band, - 'cal_band': cal_band, - 'color_mags': color_mags, - 'color_index': color_index, - 'objids': objids, - 'distances': distances, + "zp": zp, + "C": C, + "unc": unc, + "m": m_cal, + "m_inst": m_inst, + "obs_band": obs_band, + "cal_band": cal_band, + "color_mags": color_mags, + "color_index": color_index, + "objids": objids, + "distances": distances, } -def plot_color_correction( - color_mags, - m, - m_inst, - C, - zp, - color_index: str -): +def plot_color_correction(color_mags, m, m_inst, C, zp, color_index: str): """ Plot the relation between instrumental and calibrated magnitudes. @@ -259,12 +255,11 @@ def plot_color_correction( Matplotlib figure and axis objects for the plot. """ fig, ax = plt.subplots() - ax.scatter(color_mags, m - m_inst, marker='.') + ax.scatter(color_mags, m - m_inst, marker=".") x = np.linspace(0, 1.5, 100) - ax.plot(x, C * x + zp, color='red', - label=f'$m = C\\times({color_index}) + ZP$') - ax.set_xlabel(f'${color_index}$ (mag)') - ax.set_ylabel(r'$m - m_{\mathrm{inst}}$ (mag)') + ax.plot(x, C * x + zp, color="red", label=f"$m = C\\times({color_index}) + ZP$") + ax.set_xlabel(f"${color_index}$ (mag)") + ax.set_ylabel(r"$m - m_{\mathrm{inst}}$ (mag)") plt.tight_layout() return fig, ax @@ -291,22 +286,55 @@ def plot_image(telescope_image_sub, source_list, matched_idx, colored_idx): """ fig, ax = plt.subplots() m, s = np.mean(telescope_image_sub), np.std(telescope_image_sub) - im = ax.imshow(telescope_image_sub, interpolation='nearest', - origin='lower', cmap='gray') - im.set_clim(vmin=m-s, vmax=m+s) + im = ax.imshow( + telescope_image_sub, interpolation="nearest", origin="lower", cmap="gray" + ) + im.set_clim(vmin=m - s, vmax=m + s) fig.colorbar(im, ax=ax) - ax.plot(source_list['x'], source_list['y'], '+', - markersize=5, label='Detected', color='red',) - ax.plot(source_list['x'].iloc[matched_idx], source_list['y'].iloc[matched_idx], - 'o', markersize=10, color='blue', markerfacecolor='none', label='Matched') - ax.plot(source_list['x'].iloc[colored_idx], source_list['y'].iloc[colored_idx], 'o', - markersize=15, color='yellow', markerfacecolor='none', label='Selected for Color Corr') + ax.plot( + source_list["x"], + source_list["y"], + "+", + markersize=5, + label="Detected", + color="red", + ) + ax.plot( + source_list["x"].iloc[matched_idx], + source_list["y"].iloc[matched_idx], + "o", + markersize=10, + color="blue", + markerfacecolor="none", + label="Matched", + ) + ax.plot( + source_list["x"].iloc[colored_idx], + source_list["y"].iloc[colored_idx], + "o", + markersize=15, + color="yellow", + markerfacecolor="none", + label="Selected for Color Corr", + ) ax.legend() return fig, ax -def create_header(image, wcs_solution, zp, unc, source_list, matched_idx, colored_idx, input_fits, cal_band: str, catalog: str, obj_band: str): +def create_header( + image, + wcs_solution, + zp, + unc, + source_list, + matched_idx, + colored_idx, + input_fits, + cal_band: str, + catalog: str, + obj_band: str, +): """ Create and write a FITS file with calibrated header and source tables. @@ -334,26 +362,35 @@ def create_header(image, wcs_solution, zp, unc, source_list, matched_idx, colore None """ image_arr = np.asarray(image) - primary_hdu = fits.PrimaryHDU( - data=image_arr, header=wcs_solution.to_header()) - primary_hdu.header['ZP'] = zp - primary_hdu.header['ZP_STD'] = unc - primary_hdu.header['SUV_FLT'] = cal_band - primary_hdu.header['REF_CATA'] = catalog - primary_hdu.header['REF_FLT'] = obj_band - primary_hdu.header['CAT_COR'] = f"{cal_band}-{obj_band}" + primary_hdu = fits.PrimaryHDU(data=image_arr, header=wcs_solution.to_header()) + primary_hdu.header["ZP"] = zp + primary_hdu.header["ZP_STD"] = unc + primary_hdu.header["SUV_FLT"] = cal_band + primary_hdu.header["REF_CATA"] = catalog + primary_hdu.header["REF_FLT"] = obj_band + primary_hdu.header["CAT_COR"] = f"{cal_band}-{obj_band}" source_list_clean = source_list.map( - lambda x: x.filled(np.nan) if hasattr(x, 'filled') else x) - detected_hdu = fits.BinTableHDU(Table.from_pandas( - source_list_clean), name='DETECTED_SOURCES') + lambda x: x.filled(np.nan) if hasattr(x, "filled") else x + ) + detected_hdu = fits.BinTableHDU( + Table.from_pandas(source_list_clean), name="DETECTED_SOURCES" + ) if not source_list_clean.empty: - matched_hdu = fits.BinTableHDU(Table.from_pandas( - source_list_clean.iloc[matched_idx].reset_index(drop=True)), name='SELECTED_STARS') - colored_hdu = fits.BinTableHDU(Table.from_pandas( - source_list_clean.iloc[colored_idx].reset_index(drop=True)), name='START_COLOR_CORRECTION') + matched_hdu = fits.BinTableHDU( + Table.from_pandas( + source_list_clean.iloc[matched_idx].reset_index(drop=True) + ), + name="SELECTED_STARS", + ) + colored_hdu = fits.BinTableHDU( + Table.from_pandas( + source_list_clean.iloc[colored_idx].reset_index(drop=True) + ), + name="START_COLOR_CORRECTION", + ) else: - matched_hdu = fits.BinTableHDU(name='SELECTED_STARS') - colored_hdu = fits.BinTableHDU(name='START_COLOR_CORRECTION') + matched_hdu = fits.BinTableHDU(name="SELECTED_STARS") + colored_hdu = fits.BinTableHDU(name="START_COLOR_CORRECTION") hdul = fits.HDUList([primary_hdu, detected_hdu, matched_hdu, colored_hdu]) hdul.writeto(input_fits, overwrite=True) @@ -371,8 +408,18 @@ def cleanup_files(file_base): ------- None """ - extensions = ['.axy', '.corr', '.match', '.new', '.rdls', '.solved', - '-ngc.png', '-objs.png', '-indx.png', '-indx.xyls'] + extensions = [ + ".axy", + ".corr", + ".match", + ".new", + ".rdls", + ".solved", + "-ngc.png", + "-objs.png", + "-indx.png", + "-indx.xyls", + ] for ext in extensions: fname = f"{file_base}{ext}" if os.path.exists(fname): @@ -380,27 +427,40 @@ def cleanup_files(file_base): if __name__ == "__main__": - parser = argparse.ArgumentParser( - description="Photometric calibration on a background-subtracted image.") + description="Photometric calibration on a background-subtracted image." + ) + parser.add_argument("input_fits", help="Path to background-subtracted FITS image") + parser.add_argument("--Ra", type=float, required=True, help="RA from CATCH") + parser.add_argument("--Dec", type=float, required=True, help="DEC from CATCH") parser.add_argument( - "input_fits", help="Path to background-subtracted FITS image") - parser.add_argument("--Ra", type=float, required=True, - help="RA from CATCH") - parser.add_argument("--Dec", type=float, required=True, - help="DEC from CATCH") - parser.add_argument("--bkg_err", type=float, required=True, - help="Global background RMS (float, required)") - parser.add_argument("--pixel_scale", type=float, default=1.86, - help="Pixel scale in arcsec/pixel (default: 1.86)") - parser.add_argument("--snr", type=float, default=7.0, - help="Detection S/N threshold (default: 7.0)") - parser.add_argument("--catalog", default="PanSTARRS1", - help="Photometric reference catalog (default: PanSTARRS1)") - parser.add_argument("--obs_band", default="obs_band", - help="Observed image bandpass (used for labeling only; default: 'obs_band')") - parser.add_argument("--cal_band", default="r", - help="Reference catalog bandpass (default: r)") + "--bkg_err", + type=float, + required=True, + help="Global background RMS (float, required)", + ) + parser.add_argument( + "--pixel_scale", + type=float, + default=1.86, + help="Pixel scale in arcsec/pixel (default: 1.86)", + ) + parser.add_argument( + "--snr", type=float, default=7.0, help="Detection S/N threshold (default: 7.0)" + ) + parser.add_argument( + "--catalog", + default="PanSTARRS1", + help="Photometric reference catalog (default: PanSTARRS1)", + ) + parser.add_argument( + "--obs_band", + default="obs_band", + help="Observed image bandpass (used for labeling only; default: 'obs_band')", + ) + parser.add_argument( + "--cal_band", default="r", help="Reference catalog bandpass (default: r)" + ) args = parser.parse_args() input_fits = args.input_fits @@ -438,9 +498,9 @@ def cleanup_files(file_base): unc = calibration["unc"] m = calibration["m"] m_inst = calibration["m_inst"] - color_mags = calibration['color_mags'] - color_index = calibration['color_index'] - objids = calibration['objids'] + color_mags = calibration["color_mags"] + color_index = calibration["color_index"] + objids = calibration["objids"] if hasattr(objids, "mask"): matched_idx = np.where(~objids.mask)[0] @@ -451,10 +511,8 @@ def cleanup_files(file_base): else: colored_idx = np.arange(len(source_list)) - fig1, ax1 = plot_color_correction( - color_mags, m, m_inst, C, zp, color_index) - fig2, ax2 = plot_image(telescope_image_sub, - source_list, matched_idx, colored_idx) + fig1, ax1 = plot_color_correction(color_mags, m, m_inst, C, zp, color_index) + fig2, ax2 = plot_image(telescope_image_sub, source_list, matched_idx, colored_idx) plt.show() create_header( diff --git a/catch_analysis_tools/background.py b/catch_analysis_tools/background.py index 340022b..20944db 100644 --- a/catch_analysis_tools/background.py +++ b/catch_analysis_tools/background.py @@ -1,20 +1,15 @@ import numpy as np -from matplotlib import pyplot as plt -from astropy.io import fits -from photutils.segmentation import detect_sources, detect_threshold, make_2dgaussian_kernel, SourceFinder, SourceCatalog, make_2dgaussian_kernel -from astropy.stats import sigma_clipped_stats, SigmaClip +from astropy.stats import SigmaClip, sigma_clipped_stats +from photutils.background import Background2D from photutils.background import SExtractorBackground as SourceExtractorBackground -from photutils.background import Background2D, MedianBackground +from photutils.segmentation import ( + detect_sources, + detect_threshold, +) from photutils.utils import circular_footprint -from astropy.convolution import convolve -from astropy.coordinates import SkyCoord, ICRS -from astropy.wcs import WCS -from photutils.aperture import aperture_photometry, CircularAnnulus, CircularAperture -from photutils.centroids import centroid_quadratic, centroid_sources, centroid_com -from astropy import units as u -def get_background(data): +def get_background(data): """computes and returns a global background subtraction, masking sources using image segmentation IDs. Does NOT return background subtracted data @@ -24,7 +19,7 @@ def get_background(data): data : array_like 2D image array to compute background on - + Returns ------- @@ -36,21 +31,26 @@ def get_background(data): # performs a global background subtraction, masking sources using image segmentation IDs sigma_clip = SigmaClip(sigma=3.0, maxiters=10) - + threshold = detect_threshold(data, nsigma=2.0, sigma_clip=sigma_clip) segment_img = detect_sources(data, threshold, npixels=10) footprint = circular_footprint(radius=10) mask = segment_img.make_source_mask(footprint=footprint) - - + bkg_estimator = SourceExtractorBackground(sigma_clip) - bkg = Background2D(data, (50, 50), filter_size=(3, 3), - sigma_clip=sigma_clip, mask=mask, bkg_estimator=bkg_estimator) - + bkg = Background2D( + data, + (50, 50), + filter_size=(3, 3), + sigma_clip=sigma_clip, + mask=mask, + bkg_estimator=bkg_estimator, + ) return bkg + def global_subtraction(data): """performs a global background subtraction, masking sources using image segmentation IDs @@ -60,26 +60,26 @@ def global_subtraction(data): data : array_like 2D image array to be background subtracted - + Returns ------- data_sub : array_like Background subtracted data array - bkg : - background object returned from get_background(data) + bkg : + background object returned from get_background(data) """ bkg = get_background(data) - data_sub = data-bkg.background_median + data_sub = data - bkg.background_median - return data_sub, bkg + return data_sub, bkg -def calc_bkg(data,background_aperture,method,sigma_clip): - - """ Takes in an image array and aperture for background estimate. Determines background + +def calc_bkg(data, background_aperture, method, sigma_clip): + """Takes in an image array and aperture for background estimate. Determines background estimator per specified method (mean or median) and variance of pixels, with optional sigma clipping. @@ -102,43 +102,45 @@ def calc_bkg(data,background_aperture,method,sigma_clip): ------- bkg_estimator : float estimate of true background level in background aperture, per "method" - + bkg_var : float - Variance of pixels within the background aperture. If sigma clipped, square of + Variance of pixels within the background aperture. If sigma clipped, square of clipped stddev. Otherwise, as square of range of values from 50-16th percentiles """ - - background_mask = background_aperture.to_mask(method='center') - + + background_mask = background_aperture.to_mask(method="center") + background_data = background_mask.multiply(data) background_data_1d = background_data[background_mask.data > 0] - - - #Todo: Add flexibility for the stats for computing background. Include rectangular aperture. + + # Todo: Add flexibility for the stats for computing background. Include rectangular aperture. if sigma_clip is not None: - bkg_mean,bkg_median,bkg_stddev = sigma_clipped_stats(background_data_1d, sigma=sigma_clip, maxiters=10) + bkg_mean, bkg_median, bkg_stddev = sigma_clipped_stats( + background_data_1d, sigma=sigma_clip, maxiters=10 + ) bkg_var = bkg_stddev**2 else: bkg_median = np.nanmedian(background_data_1d) bkg_mean = np.nanmean(background_data_1d) # robust approximation via https://wise2.ipac.caltech.edu/staff/fmasci/ApPhotUncert.pdf - bkg_var = np.square(np.percentile(background_data_1d,50)-np.percentile(background_data_1d,16)) + bkg_var = np.square( + np.percentile(background_data_1d, 50) + - np.percentile(background_data_1d, 16) + ) bkg_stddev = np.sqrt(bkg_var) - - if method=='mean': + + if method == "mean": bkg_estimator = bkg_mean - elif method=='median': + elif method == "median": bkg_estimator = bkg_median else: raise ValueError("Method must be 'mean' or 'median'") - - #plt.figure(figsize=(8,8)) + # plt.figure(figsize=(8,8)) # optional code to check the standard deviation estimate of the background - #plt.axvspan(bkg_estimator-bkg_stddev,bkg_estimator+bkg_stddev,alpha=0.5) - #plt.axvline(bkg_estimator,color='red') - #plt.hist(annulus_data_1d) - - + # plt.axvspan(bkg_estimator-bkg_stddev,bkg_estimator+bkg_stddev,alpha=0.5) + # plt.axvline(bkg_estimator,color='red') + # plt.hist(annulus_data_1d) + return bkg_estimator, bkg_var diff --git a/catch_analysis_tools/calibration/astrometry.py b/catch_analysis_tools/calibration/astrometry.py index d8f8dca..1022b55 100644 --- a/catch_analysis_tools/calibration/astrometry.py +++ b/catch_analysis_tools/calibration/astrometry.py @@ -1,17 +1,20 @@ +import argparse import os import subprocess -import argparse + +import fitsio import numpy as np import pandas as pd import sep -import fitsio +from astropy.coordinates import SkyCoord from astropy.io import fits -from astropy.wcs import WCS from astropy.table import Table -from astropy.coordinates import SkyCoord +from astropy.wcs import WCS -def run_solve_field(input_fits, output_wcs, pixel_scale, Ra_deg, Dec_deg, scale_units="arcsecperpix"): +def run_solve_field( + input_fits, output_wcs, pixel_scale, Ra_deg, Dec_deg, scale_units="arcsecperpix" +): """ Execute the `solve-field` command to compute a WCS solution. @@ -33,27 +36,35 @@ def run_solve_field(input_fits, output_wcs, pixel_scale, Ra_deg, Dec_deg, scale_ """ if os.path.exists(output_wcs): print( - f"Output file '{output_wcs}' already exists. Skipping solve-field execution.") + f"Output file '{output_wcs}' already exists. Skipping solve-field execution." + ) return True config_file = os.environ.get("ASTROMETRY_CONFIG") if config_file is None: raise RuntimeError( - "ASTROMETRY_CONFIG is not set. " - "This is required to run solve-field." + "ASTROMETRY_CONFIG is not set. This is required to run solve-field." ) command = [ "solve-field", "--overwrite", - "--config", config_file, - "--ra", str(Ra_deg), - "--dec", str(Dec_deg), - "--scale-units", scale_units, - "--scale-low", str(pixel_scale * 0.5), - "--scale-high", str(pixel_scale * 2.0), - "--radius", "2", - "--downsample", "1", + "--config", + config_file, + "--ra", + str(Ra_deg), + "--dec", + str(Dec_deg), + "--scale-units", + scale_units, + "--scale-low", + str(pixel_scale * 0.5), + "--scale-high", + str(pixel_scale * 2.0), + "--radius", + "2", + "--downsample", + "1", input_fits, ] @@ -87,23 +98,14 @@ def find_sources(image_sub, bkg_err, snr, aperture_radius=7.0): Background-subtracted image array. """ sep.set_sub_object_limit(500) - sources = sep.extract( - image_sub, - thresh=snr, - err=bkg_err, - deblend_nthresh=16 - ) + sources = sep.extract(image_sub, thresh=snr, err=bkg_err, deblend_nthresh=16) source_list = pd.DataFrame(sources) flux, flux_err, _ = sep.sum_circle( - image_sub, - source_list['x'], source_list['y'], - aperture_radius, - err=bkg_err + image_sub, source_list["x"], source_list["y"], aperture_radius, err=bkg_err ) - source_list['aperture_sum'] = flux - source_list['aperture_err'] = flux_err - source_list = source_list[source_list['aperture_sum'] > 0].reset_index( - drop=True) + source_list["aperture_sum"] = flux + source_list["aperture_err"] = flux_err + source_list = source_list[source_list["aperture_sum"] > 0].reset_index(drop=True) return source_list, image_sub @@ -146,10 +148,10 @@ def retrieve_sources(source_list, wcs_solution): sky_coords : astropy.coordinates.SkyCoord SkyCoord object with celestial coordinates of sources. """ - world = wcs_solution.pixel_to_world(source_list['x'], source_list['y']) - source_list['RA'] = [c.ra.deg for c in world] - source_list['Dec'] = [c.dec.deg for c in world] - sky_coords = SkyCoord(source_list['RA'], source_list['Dec'], unit='deg') + world = wcs_solution.pixel_to_world(source_list["x"], source_list["y"]) + source_list["RA"] = [c.ra.deg for c in world] + source_list["Dec"] = [c.dec.deg for c in world] + sky_coords = SkyCoord(source_list["RA"], source_list["Dec"], unit="deg") return source_list, sky_coords @@ -166,13 +168,24 @@ def cleanup_files(file_base): ------- None """ - extensions = ['.axy', '.corr', '.match', '.new', '.rdls', '.solved', - '-ngc.png', '-objs.png', '-indx.png', '-indx.xyls'] + extensions = [ + ".axy", + ".corr", + ".match", + ".new", + ".rdls", + ".solved", + "-ngc.png", + "-objs.png", + "-indx.png", + "-indx.xyls", + ] for ext in extensions: fname = f"{file_base}{ext}" if os.path.exists(fname): os.remove(fname) + def write_astrometry_output( image, wcs_solution, diff --git a/catch_analysis_tools/calibration/photometry.py b/catch_analysis_tools/calibration/photometry.py index 6ccb709..a7d81d9 100644 --- a/catch_analysis_tools/calibration/photometry.py +++ b/catch_analysis_tools/calibration/photometry.py @@ -1,7 +1,7 @@ +import calviacat as cvc +import matplotlib.pyplot as plt import numpy as np import pandas as pd -import matplotlib.pyplot as plt -import calviacat as cvc from astropy.io import fits from astropy.table import Table @@ -212,6 +212,7 @@ def plot_photometric_matches( return fig, ax + def write_photometric_calibration_output( image, wcs_solution, @@ -286,4 +287,4 @@ def write_photometric_calibration_output( ] ) - hdul.writeto(output_fits, overwrite=True) \ No newline at end of file + hdul.writeto(output_fits, overwrite=True) diff --git a/catch_analysis_tools/photometry.py b/catch_analysis_tools/photometry.py index 0ede353..e960759 100644 --- a/catch_analysis_tools/photometry.py +++ b/catch_analysis_tools/photometry.py @@ -1,17 +1,20 @@ import numpy as np -from matplotlib import pyplot as plt -from astropy.io import fits -from photutils.segmentation import detect_sources, detect_threshold, make_2dgaussian_kernel, SourceFinder, SourceCatalog, make_2dgaussian_kernel -from astropy.stats import sigma_clipped_stats, SigmaClip -from photutils.background import SExtractorBackground as SourceExtractorBackground -from photutils.background import Background2D, MedianBackground -from photutils.utils import circular_footprint from astropy.convolution import convolve -from astropy.coordinates import SkyCoord, ICRS +from astropy.io import fits from astropy.wcs import WCS -from photutils.aperture import aperture_photometry, CircularAnnulus, CircularAperture, RectangularAperture +from photutils.aperture import ( + CircularAnnulus, + CircularAperture, + RectangularAperture, +) from photutils.centroids import centroid_quadratic, centroid_sources -from astropy import units as u +from photutils.segmentation import ( + SourceCatalog, + SourceFinder, + make_2dgaussian_kernel, +) +from photutils.utils import circular_footprint + # here are functions for grabbing the data, doing background subtractions and manipulating source extractions def get_image(url): @@ -23,14 +26,14 @@ def get_image(url): url : string CATCH-generated URL from a query. - + Returns ------- data : array_like This is a 2D array containing the image data returned by the CATCH query - header : - FITS header class, from astropy.io.fits.Header + header : + FITS header class, from astropy.io.fits.Header """ @@ -39,9 +42,9 @@ def get_image(url): data = fits_hdu[0].data header = fits_hdu[0].header return data, header - -def id_good_sources(data,bkg): + +def id_good_sources(data, bkg): """Uses a segmentation image to identify reliable sources in image that can be snapped to. Coincidentally, computes baseline photometry that could be used as a quality comparison user results, @@ -60,38 +63,36 @@ def id_good_sources(data,bkg): ------- cat : Astropy Table class, from SourceCatalog output giving source locations, fluxes - - """ - source_threshold = bkg.background_median + 1.5 * bkg.background_rms + """ + source_threshold = bkg.background_median + 1.5 * bkg.background_rms kernel = make_2dgaussian_kernel(3.0, size=5) # FWHM = 3.0 convolved_data = convolve(data, kernel) finder = SourceFinder(npixels=5, progress_bar=False) segment_map = finder(convolved_data, source_threshold) - - - vmax = np.percentile(np.ndarray.flatten(data),99) - vmin = np.percentile(np.ndarray.flatten(data),1) - + + vmax = np.percentile(np.ndarray.flatten(data), 99) + vmin = np.percentile(np.ndarray.flatten(data), 1) + # make a plot to show the background subtracted frame and the resulting segment map - #fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 12.5)) - #ax1.imshow(data, origin='lower', cmap='Greys_r', vmin=vmin,vmax=vmax) - #ax1.set_title('Original Data') + # fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 12.5)) + # ax1.imshow(data, origin='lower', cmap='Greys_r', vmin=vmin,vmax=vmax) + # ax1.set_title('Original Data') - #ax2.imshow(segment_map, origin='lower', cmap=segment_map.cmap, + # ax2.imshow(segment_map, origin='lower', cmap=segment_map.cmap, # interpolation='nearest') - #ax2.set_title('Segmentation Image') + # ax2.set_title('Segmentation Image') - cat = SourceCatalog(data, segment_map, convolved_data=convolved_data) - + return cat - -def create_user_aperture(position,size): + + +def create_user_aperture(position, size): """Simple placeholder function for making user-selected circular apertures @@ -105,34 +106,44 @@ def create_user_aperture(position,size): Returns ------- - aperture : + aperture : Photutils circular aperture object with the specified size, location parameters - - """ - + + """ + aperture = CircularAperture(position, r=size) return aperture + def define_aperture(aperture_params): """ Defines an aperture from user supplied parameters, used both to create target and background apertures """ - if aperture_params["shape"]=="Circular": - aperture = CircularAperture(aperture_params["position"], r=aperture_params["size"]) - elif aperture_params["shape"]=="Circular_Annulus": - aperture = CircularAnnulus(aperture_params["position"], r_in=aperture_params["inner_r"], r_out=aperture_params["outer_r"]) - elif aperture_params["shape"]=="Rectangular": - aperture = RectangularAperture(aperture_params["position"], aperture_params["w"], aperture_params["h"], theta=aperture_params["theta"]) + if aperture_params["shape"] == "Circular": + aperture = CircularAperture( + aperture_params["position"], r=aperture_params["size"] + ) + elif aperture_params["shape"] == "Circular_Annulus": + aperture = CircularAnnulus( + aperture_params["position"], + r_in=aperture_params["inner_r"], + r_out=aperture_params["outer_r"], + ) + elif aperture_params["shape"] == "Rectangular": + aperture = RectangularAperture( + aperture_params["position"], + aperture_params["w"], + aperture_params["h"], + theta=aperture_params["theta"], + ) else: aperture = 0 return aperture - - -def subpixel_centroid(user_point,data,radius): +def subpixel_centroid(user_point, data, radius): """Takes in a user defined point and returns the location of the brightest pixel within radius pixels @@ -151,23 +162,29 @@ def subpixel_centroid(user_point,data,radius): ------- target_position: array_like [x,y] location of the source that is closest to user_point found in data - - """ + + """ footprint = circular_footprint(radius) - x, y = centroid_sources(data, user_point[0], user_point[1], footprint=footprint, - centroid_func=centroid_quadratic) - - target_position = np.array([x[0],y[0]]) + x, y = centroid_sources( + data, + user_point[0], + user_point[1], + footprint=footprint, + centroid_func=centroid_quadratic, + ) + + target_position = np.array([x[0], y[0]]) return target_position - -def do_aperture_photometry(data,source_aperture, bkg_aperture): - """ Takes in an image, a source aperture, and outputs from the calc_annulus_bkg function - + + +def do_aperture_photometry(data, source_aperture, bkg_aperture): + """Takes in an image, a source aperture, and outputs from the calc_annulus_bkg function + Returns the source flux (background subtracted, per-pixel background median) and the uncertainty as defined at the quoted link - + method='center' means pixels are either in or out, no interpolation to a perfect circle (in other words, areas will be in whole pixels) @@ -176,16 +193,16 @@ def do_aperture_photometry(data,source_aperture, bkg_aperture): data : array_like image data to be used for photometry - source_aperture : + source_aperture : Photutils aperture object containing desired source bkg_median : float - median value of pixel background, - + median value of pixel background, + bkg_var : float variance of pixel background, ideally from output of calc_annulus_background - bkg_aperture : + bkg_aperture : Photutils aperture object containing background, ideally as annulus_aperture output from calc_annulus_background @@ -193,34 +210,37 @@ def do_aperture_photometry(data,source_aperture, bkg_aperture): ------- source_sum : float background subtracted flux of the targeted source - + source_err : float error on source flux """ - aperture_mask = source_aperture.to_mask(method='center') + aperture_mask = source_aperture.to_mask(method="center") aperture_data = aperture_mask.multiply(data) aperture_sum = np.nansum(aperture_data) aperture_area = np.nansum(aperture_mask) - bkg_mask = bkg_aperture.to_mask(method='center') + bkg_mask = bkg_aperture.to_mask(method="center") bkg_area = np.nansum(bkg_mask) bkg_data = bkg_mask.multiply(data) bkg_median = np.nanmedian(bkg_data) bkg_var = np.nanvar(bkg_data[bkg_data != 0]) background = bkg_median * aperture_area - source_sum = aperture_sum-background - + source_sum = aperture_sum - background + # Using uncertainty as derived by https://wise2.ipac.caltech.edu/staff/fmasci/ApPhotUncert.pdf # Setting the gain g=1, N_i = 1. Assumes data has already been converted to e- ### TODO: FIND THE GAINS FOR EACH SURVEY term1 = source_sum - term2 = (aperture_area + (np.pi/2) * (np.square(aperture_area)/bkg_area) )*bkg_var + term2 = ( + aperture_area + (np.pi / 2) * (np.square(aperture_area) / bkg_area) + ) * bkg_var source_err = np.sqrt(term1 + term2) - + return source_sum, source_err + def load_thumbnail(url): """Access a cutout via a call to a CATCH URL, returning WCSobject @@ -230,15 +250,15 @@ def load_thumbnail(url): url : string CATCH-generated URL from a query. - + Returns ------- data : array_like This is a 2D array containing the image data returned by the CATCH query - header : - FITS header class, from astropy.io.fits.Header + header : + FITS header class, from astropy.io.fits.Header img_WCS: Astropy WCS object @@ -249,18 +269,18 @@ def load_thumbnail(url): header = fits_hdu[0].header img_WCS = WCS(header) return data, header, img_WCS - -def source_instr_mag(ap_flux,ap_fluxerr,exposure_time): - """ Quick function to return instrumental magnitudes from a source flux + +def source_instr_mag(ap_flux, ap_fluxerr, exposure_time): + """Quick function to return instrumental magnitudes from a source flux Does not force magnitude uncertainties to be symmetric - + To be used by calibrated_mag() function Parameters ---------- - + ap_flux: float Flux (in counts) of source @@ -278,22 +298,23 @@ def source_instr_mag(ap_flux,ap_fluxerr,exposure_time): as [Magnitude, Upper Magnitude Uncertainty, Lower Magnitude Uncertainty] """ - instr_mag = -2.5*np.log10(ap_flux/exposure_time) - instr_mag_hi = -2.5*np.log10((ap_flux-ap_fluxerr)/exposure_time) - instr_mag_lo = -2.5*np.log10((ap_flux+ap_fluxerr)/exposure_time) - + instr_mag = -2.5 * np.log10(ap_flux / exposure_time) + instr_mag_hi = -2.5 * np.log10((ap_flux - ap_fluxerr) / exposure_time) + instr_mag_lo = -2.5 * np.log10((ap_flux + ap_fluxerr) / exposure_time) + instr_mag_hi_uncert = instr_mag_hi - instr_mag instr_mag_lo_uncert = instr_mag_lo - instr_mag - - instr_mag_array = np.array([instr_mag,instr_mag_hi_uncert,instr_mag_lo_uncert]) - + + instr_mag_array = np.array([instr_mag, instr_mag_hi_uncert, instr_mag_lo_uncert]) + return instr_mag_array -def calibrated_mag(instr_mag_array,zero_point,zero_point_uncert): - """ Takes in the array from source_instr_mag, converts to derived magnitude, + +def calibrated_mag(instr_mag_array, zero_point, zero_point_uncert): + """Takes in the array from source_instr_mag, converts to derived magnitude, propagating uncertainties from both - + Parameters ---------- instr_mag_array: array_like @@ -305,7 +326,7 @@ def calibrated_mag(instr_mag_array,zero_point,zero_point_uncert): zero_point_uncert: float zero point uncertainty (mags) of image (ideally taken from metadata in header given by astrometry solution) - + Returns ------- @@ -314,13 +335,12 @@ def calibrated_mag(instr_mag_array,zero_point,zero_point_uncert): as [Magnitude, Upper Magnitude Uncertainty, Lower Magnitude Uncertainty] """ - calib_mag = zero_point+instr_mag_array[0] + calib_mag = zero_point + instr_mag_array[0] calib_mag_hi = np.sqrt(np.square(zero_point_uncert) + np.square(instr_mag_array[1])) calib_mag_lo = np.sqrt(np.square(zero_point_uncert) + np.square(instr_mag_array[2])) calib_mag_array = { "cal_mag": calib_mag, "cal_mag_hi_uncert": calib_mag_hi, - "cal_mag_lo_uncert": calib_mag_lo + "cal_mag_lo_uncert": calib_mag_lo, } - return calib_mag_array - + return calib_mag_array diff --git a/catch_analysis_tools/tests/astrometry_test.py b/catch_analysis_tools/tests/astrometry_test.py index ae2f597..cc2affe 100644 --- a/catch_analysis_tools/tests/astrometry_test.py +++ b/catch_analysis_tools/tests/astrometry_test.py @@ -1,14 +1,16 @@ -import pytest import shutil import tempfile from pathlib import Path -from unittest.mock import patch, MagicMock +from unittest.mock import MagicMock, patch + +import pytest from ..astrometry import * RA_DEG = 263.0 DEC_DEG = 34.5 + @pytest.fixture def synthetic_image(): image = np.random.normal(loc=0, scale=1.0, size=(100, 100)).astype(np.float32) @@ -21,30 +23,28 @@ def synthetic_image(): @pytest.fixture def synthetic_wcs(): header = fits.Header() - header['NAXIS'] = 2 - header['NAXIS1'] = 100 - header['NAXIS2'] = 100 - header['CTYPE1'] = 'RA---TAN' - header['CTYPE2'] = 'DEC--TAN' - header['CRPIX1'] = 50.0 - header['CRPIX2'] = 50.0 - header['CRVAL1'] = 150.0 - header['CRVAL2'] = 2.0 - header['CD1_1'] = -0.000277778 - header['CD1_2'] = 0.0 - header['CD2_1'] = 0.0 - header['CD2_2'] = 0.000277778 + header["NAXIS"] = 2 + header["NAXIS1"] = 100 + header["NAXIS2"] = 100 + header["CTYPE1"] = "RA---TAN" + header["CTYPE2"] = "DEC--TAN" + header["CRPIX1"] = 50.0 + header["CRPIX2"] = 50.0 + header["CRVAL1"] = 150.0 + header["CRVAL2"] = 2.0 + header["CD1_1"] = -0.000277778 + header["CD1_2"] = 0.0 + header["CD2_1"] = 0.0 + header["CD2_2"] = 0.000277778 return WCS(header) @pytest.mark.integration @pytest.mark.skipif( - shutil.which("solve-field") is None or - "ASTROMETRY_CONFIG" not in os.environ, - reason="solve-field or astrometry config not available" + shutil.which("solve-field") is None or "ASTROMETRY_CONFIG" not in os.environ, + reason="solve-field or astrometry config not available", ) def test_run_solve_field_real(tmp_path): - # Use a real sky FITS image committed to the repo input_fits = Path(__file__).parent / "data" / "Comet_65P_Gunn_LONEOS.fits" assert input_fits.exists(), "test image is missing" @@ -60,16 +60,22 @@ def test_run_solve_field_real(tmp_path): Dec_deg=17.0, ) -@pytest.mark.parametrize("file_exists, should_call_run", [ - (True, False), - (False, True), -]) -def test_run_solve_field_conditional_execution(file_exists, should_call_run): - with patch.dict(os.environ, {"ASTROMETRY_CONFIG": "/fake/astrometry/config"}), \ - patch("catch_analysis_tools.astrometry.os.path.exists", - return_value=file_exists) as mock_exists, \ - patch("catch_analysis_tools.astrometry.subprocess.run") as mock_run: +@pytest.mark.parametrize( + "file_exists, should_call_run", + [ + (True, False), + (False, True), + ], +) +def test_run_solve_field_conditional_execution(file_exists, should_call_run): + with ( + patch.dict(os.environ, {"ASTROMETRY_CONFIG": "/fake/astrometry/config"}), + patch( + "catch_analysis_tools.astrometry.os.path.exists", return_value=file_exists + ) as mock_exists, + patch("catch_analysis_tools.astrometry.subprocess.run") as mock_run, + ): result = run_solve_field( "input.fits", "output.wcs", @@ -88,14 +94,14 @@ def test_run_solve_field_conditional_execution(file_exists, should_call_run): def test_run_solve_field_raises_if_subprocess_fails(): - with patch.dict(os.environ, {"ASTROMETRY_CONFIG": "/fake/astrometry/config"}), \ - patch("catch_analysis_tools.astrometry.os.path.exists", - return_value=False), \ - patch( - "catch_analysis_tools.astrometry.subprocess.run", - side_effect=subprocess.CalledProcessError(1, "solve-field") - ): - + with ( + patch.dict(os.environ, {"ASTROMETRY_CONFIG": "/fake/astrometry/config"}), + patch("catch_analysis_tools.astrometry.os.path.exists", return_value=False), + patch( + "catch_analysis_tools.astrometry.subprocess.run", + side_effect=subprocess.CalledProcessError(1, "solve-field"), + ), + ): with pytest.raises(RuntimeError, match="solve-field failed"): run_solve_field( "input.fits", @@ -119,43 +125,59 @@ def test_load_wcs(synthetic_wcs): with tempfile.NamedTemporaryFile(suffix=".fits", delete=False) as tmp: filename = tmp.name try: - fits.PrimaryHDU(header=synthetic_wcs.to_header()).writeto(filename, overwrite=True) + fits.PrimaryHDU(header=synthetic_wcs.to_header()).writeto( + filename, overwrite=True + ) wcs = load_wcs(filename) assert isinstance(wcs, WCS) - assert wcs.wcs.ctype[0] == 'RA---TAN' + assert wcs.wcs.ctype[0] == "RA---TAN" finally: os.remove(filename) def test_retrieve_sources(synthetic_wcs): - df = pd.DataFrame({'x': [50.0], 'y': [50.0]}) + df = pd.DataFrame({"x": [50.0], "y": [50.0]}) df_out, sky_coords = retrieve_sources(df.copy(), synthetic_wcs) assert "RA" in df_out.columns and "Dec" in df_out.columns assert isinstance(sky_coords, SkyCoord) def test_calibrate_photometry(): - df = pd.DataFrame({'aperture_sum': [1000, 2000, 500]}) - coords = SkyCoord([150.0, 150.01, 150.02], [2.0, 2.01, 2.02], unit='deg') + df = pd.DataFrame({"aperture_sum": [1000, 2000, 500]}) + coords = SkyCoord([150.0, 150.01, 150.02], [2.0, 2.01, 2.02], unit="deg") mock_cat = MagicMock() mock_cat.search.return_value = (np.arange(3),) mock_cat.xmatch.return_value = (np.array([1, 2, 3]), np.array([0.1, 0.2, 0.3])) - mock_cat.cal_color.return_value = (25.0, 0.05, 0.01, np.array([20.1, 20.2, 20.3]), np.array([0.3, 0.2, 0.1]), None) + mock_cat.cal_color.return_value = ( + 25.0, + 0.05, + 0.01, + np.array([20.1, 20.2, 20.3]), + np.array([0.3, 0.2, 0.1]), + None, + ) with patch("catch_analysis_tools.astrometry.cvc") as mock_cvc: mock_cvc.PanSTARRS1.return_value = mock_cat result = calibrate_photometry(coords, df) - assert result['zp'] == 25.0 and len(result['m']) == 3 + assert result["zp"] == 25.0 and len(result["m"]) == 3 def test_plot_color_correction(): - fig, ax = plot_color_correction(np.array([0.1, 0.4]), np.array([18.3, 18.8]), np.array([18.0, 18.5]), 0.1, 25.0, "r-g") + fig, ax = plot_color_correction( + np.array([0.1, 0.4]), + np.array([18.3, 18.8]), + np.array([18.0, 18.5]), + 0.1, + 25.0, + "r-g", + ) assert fig and ax plt.close(fig) def test_plot_image(): image = np.random.normal(100, 10, (50, 50)) - df = pd.DataFrame({'x': [10, 20, 30], 'y': [10, 25, 35]}) + df = pd.DataFrame({"x": [10, 20, 30], "y": [10, 25, 35]}) fig, ax = plot_image(image, df, [0, 1], [1, 2]) assert fig and ax plt.close(fig) @@ -163,18 +185,41 @@ def test_plot_image(): def test_create_header_writes_fits(tmp_path, synthetic_wcs): image = np.random.rand(10, 10) - df = pd.DataFrame({'x': [1, 2, 3], 'y': [4, 5, 6], 'aperture_sum': [100, 200, 300]}) + df = pd.DataFrame({"x": [1, 2, 3], "y": [4, 5, 6], "aperture_sum": [100, 200, 300]}) matched_idx, colored_idx = [0, 1], [1, 2] output_path = tmp_path / "output.fits" - create_header(image, synthetic_wcs, 25.0, 0.02, df, matched_idx, colored_idx, str(output_path), "r", "PanSTARRS1", "g") + create_header( + image, + synthetic_wcs, + 25.0, + 0.02, + df, + matched_idx, + colored_idx, + str(output_path), + "r", + "PanSTARRS1", + "g", + ) assert output_path.exists() with fits.open(output_path) as hdul: - assert hdul[0].header['ZP'] == 25.0 and hdul[1].name == 'DETECTED_SOURCES' + assert hdul[0].header["ZP"] == 25.0 and hdul[1].name == "DETECTED_SOURCES" def test_cleanup_files(tmp_path): file_base = tmp_path / "testfile" - extensions = ['.axy', '.corr', '.match', '.new', '.rdls', '.solved', '-ngc.png', '-objs.png', '-indx.png', '-indx.xyls'] + extensions = [ + ".axy", + ".corr", + ".match", + ".new", + ".rdls", + ".solved", + "-ngc.png", + "-objs.png", + "-indx.png", + "-indx.xyls", + ] for ext in extensions: (file_base.with_name(file_base.name + ext)).write_text("tmp") cleanup_files(str(file_base)) diff --git a/catch_analysis_tools/tests/background_test.py b/catch_analysis_tools/tests/background_test.py index 9053985..f04a485 100644 --- a/catch_analysis_tools/tests/background_test.py +++ b/catch_analysis_tools/tests/background_test.py @@ -1,29 +1,42 @@ import numpy as np -import pytest from pytest import approx -from ..photometry import define_aperture + from ..background import * +from ..photometry import define_aperture + def test_global_subtraction(): - data = np.ones((100,100)) - data[40:50,40:50] = 5*data[40:50,40:50] + data = np.ones((100, 100)) + data[40:50, 40:50] = 5 * data[40:50, 40:50] data_sub, bkg = global_subtraction(data) errors = [] # test some basic values checked manually against above image values assert np.mean(data_sub) == 0.04 assert np.mean(bkg.background) == 1.0 - + def test_get_background(): - data = np.ones((100,100)) - data[40:50,40:50] = 5*data[40:50,40:50] + data = np.ones((100, 100)) + data[40:50, 40:50] = 5 * data[40:50, 40:50] bkg = get_background(data) assert np.mean(bkg.background) == 1.0 + def test_calc_bkg(): import photutils.datasets - noise = photutils.datasets.make_noise_image((100,100), distribution='gaussian', mean=5, stddev=1, seed=1) - bkg_aperture = define_aperture({"shape":"Circular_Annulus","position":[50,50],"size":5,"inner_r":1,"outer_r":40}) - bkg_median, bkg_var= calc_bkg(noise,bkg_aperture,"median",sigma_clip=None) - assert approx([bkg_median,bkg_var]) == [4.990637730701752, 0.9217670741324576] + + noise = photutils.datasets.make_noise_image( + (100, 100), distribution="gaussian", mean=5, stddev=1, seed=1 + ) + bkg_aperture = define_aperture( + { + "shape": "Circular_Annulus", + "position": [50, 50], + "size": 5, + "inner_r": 1, + "outer_r": 40, + } + ) + bkg_median, bkg_var = calc_bkg(noise, bkg_aperture, "median", sigma_clip=None) + assert approx([bkg_median, bkg_var]) == [4.990637730701752, 0.9217670741324576] diff --git a/catch_analysis_tools/tests/photometry_test.py b/catch_analysis_tools/tests/photometry_test.py index e0c572d..a9550ff 100644 --- a/catch_analysis_tools/tests/photometry_test.py +++ b/catch_analysis_tools/tests/photometry_test.py @@ -1,74 +1,103 @@ import numpy as np import pytest -from pytest import approx from photutils.datasets import make_wcs +from pytest import approx -# here are test functions for grabbing the data, doing background subtractions and manipulating source extractions +from ..app.services.photometry import * +from ..background import * +# here are test functions for grabbing the data, doing background subtractions and manipulating source extractions from ..photometry import * -from ..background import * -from ..app.services.photometry import * + + @pytest.mark.remote_data def test_image(): - test_url = 'https://sbnsurveys.astro.umd.edu/api/images/urn%3Anasa%3Apds%3Agbo.ast.neat.survey%3Adata_tricam%3Ap20020121_obsdata_20020121132624c?format=fits&size=10.00arcmin&ra=177.51011&dec=15.25013' + test_url = "https://sbnsurveys.astro.umd.edu/api/images/urn%3Anasa%3Apds%3Agbo.ast.neat.survey%3Adata_tricam%3Ap20020121_obsdata_20020121132624c?format=fits&size=10.00arcmin&ra=177.51011&dec=15.25013" data, header = get_image(test_url) # test some basic values checked manually against above image values - assert approx(data[10,10]) == 1558.4028 - assert approx(header['SECZ']) == 1.055449 + assert approx(data[10, 10]) == 1558.4028 + assert approx(header["SECZ"]) == 1.055449 + @pytest.mark.remote_data def test_id_good_sources(): - test_url = 'https://sbnsurveys.astro.umd.edu/api/images/urn%3Anasa%3Apds%3Agbo.ast.neat.survey%3Adata_tricam%3Ap20020121_obsdata_20020121132624c?format=fits&size=10.00arcmin&ra=177.51011&dec=15.25013' + test_url = "https://sbnsurveys.astro.umd.edu/api/images/urn%3Anasa%3Apds%3Agbo.ast.neat.survey%3Adata_tricam%3Ap20020121_obsdata_20020121132624c?format=fits&size=10.00arcmin&ra=177.51011&dec=15.25013" data, header = get_image(test_url) bkg = get_background(data) - cat = id_good_sources(data,bkg) + cat = id_good_sources(data, bkg) assert len(cat.to_table()) == 125 + def test_create_aperture(): - aperture = create_user_aperture((11,54),12) + aperture = create_user_aperture((11, 54), 12) assert approx(aperture.area) == 452.3893421169302 + def test_define_aperture(): - aperture = define_aperture(({"shape":"Circular","position":[11,54],"size":12})) + aperture = define_aperture( + ({"shape": "Circular", "position": [11, 54], "size": 12}) + ) assert approx(aperture.area) == 452.3893421169302 + def test_subpixel_centroid(): - data = np.zeros((100,100)) - data[40:50,40:50] = 5+data[40:50,40:50] - data[45,45] = 100 - source = subpixel_centroid([41,48],data,20) - assert approx(source) == np.array([45., 45.]) + data = np.zeros((100, 100)) + data[40:50, 40:50] = 5 + data[40:50, 40:50] + data[45, 45] = 100 + source = subpixel_centroid([41, 48], data, 20) + assert approx(source) == np.array([45.0, 45.0]) + def test_do_aperture_photometry(): import photutils.datasets - noise = photutils.datasets.make_noise_image((100,100), distribution='gaussian', mean=0, stddev=1, seed=1) - noise[50,50] = noise[50,50]+ 100 # give us a "source" with flux of 100 to recover - source_aperture = define_aperture({"shape":"Circular","position":[50,50],"size":5}) - bkg_aperture = define_aperture({"shape":"Circular_Annulus","position":[50,50],"size":5,"inner_r":5,"outer_r":10}) - source_sum, source_err = do_aperture_photometry(noise,source_aperture,bkg_aperture) + + noise = photutils.datasets.make_noise_image( + (100, 100), distribution="gaussian", mean=0, stddev=1, seed=1 + ) + noise[50, 50] = ( + noise[50, 50] + 100 + ) # give us a "source" with flux of 100 to recover + source_aperture = define_aperture( + {"shape": "Circular", "position": [50, 50], "size": 5} + ) + bkg_aperture = define_aperture( + { + "shape": "Circular_Annulus", + "position": [50, 50], + "size": 5, + "inner_r": 5, + "outer_r": 10, + } + ) + source_sum, source_err = do_aperture_photometry( + noise, source_aperture, bkg_aperture + ) assert approx([source_sum, source_err]) == [89.2212216869083, 14.261818958632185] + def test_get_world_coordinates(): shape = (100, 100) wcs = make_wcs(shape) - skycoord = get_world_coordinates(wcs,42, 57) - assert approx([skycoord['ra'],skycoord['dec']]) == [197.89278975, -1.36561284] + skycoord = get_world_coordinates(wcs, 42, 57) + assert approx([skycoord["ra"], skycoord["dec"]]) == [197.89278975, -1.36561284] + def test_get_pixel_coordinates(): shape = (100, 100) wcs = make_wcs(shape) - pixcoord = get_pixel_coordinates(wcs,197.89278975, -1.36561284) + pixcoord = get_pixel_coordinates(wcs, 197.89278975, -1.36561284) + + assert approx([np.round(pixcoord["x"]), np.round(pixcoord["y"])]) == [42.0, 57.0] - assert approx([np.round(pixcoord['x']),np.round(pixcoord['y'])]) == [42.0,57.0] def test_source_instr_mag(): - mag = source_instr_mag(10,1,1) + mag = source_instr_mag(10, 1, 1) assert approx(mag[0] + mag[1]) == -2.38560627 + def test_calibrated_mag(): - calib_mag = calibrated_mag(source_instr_mag(10,1,1),22,0.5) - print (calib_mag) + calib_mag = calibrated_mag(source_instr_mag(10, 1, 1), 22, 0.5) + print(calib_mag) assert approx(calib_mag["cal_mag"] + calib_mag["cal_mag_hi_uncert"]) == 20.01291902 - diff --git a/pyproject.toml b/pyproject.toml index 6509d33..f2a2182 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -44,7 +44,7 @@ packages = ["catch_analysis_tools"] write_to = "catch_analysis_tools/version.py" [tool.ruff] -select = [ +lint.select = [ "E", # pycodestyle errors "W", # pycodestyle warnings "F", # pyflakes. "E" + "W" + "F" + "C90" (mccabe complexity) is equivalent to flake8 From 200f706196939c89743d072577401440bfa8c9d9 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 16 Jun 2026 11:00:19 -0400 Subject: [PATCH 3/3] Linting and docstring reformatting. --- .../app/services/astrometry.py | 21 ++- .../app/services/photometry.py | 48 ++++--- .../services/subtract_median_background.py | 15 ++- catch_analysis_tools/astrometry.py | 72 +++++++++- catch_analysis_tools/background.py | 43 +++--- .../calibration/astrometry.py | 41 +++++- catch_analysis_tools/photometry.py | 124 ++++++++++-------- catch_analysis_tools/tests/astrometry_test.py | 22 +++- catch_analysis_tools/tests/background_test.py | 7 +- catch_analysis_tools/tests/photometry_test.py | 23 +++- pyproject.toml | 1 + scripts/Checking_indexfiles.py | 2 +- 12 files changed, 305 insertions(+), 114 deletions(-) diff --git a/catch_analysis_tools/app/services/astrometry.py b/catch_analysis_tools/app/services/astrometry.py index 1f4bcd3..2a856ff 100644 --- a/catch_analysis_tools/app/services/astrometry.py +++ b/catch_analysis_tools/app/services/astrometry.py @@ -136,23 +136,30 @@ def find_sources(image, det_cfg): """ Detect sources in an image using SEP background subtraction and extraction. + Parameters ---------- image_sub : array_like 2D numpy array after background subtraction (cleaned image). + bkg_err : float or array_like Background noise estimate (global RMS or per‐pixel error map). + snr : float Minimum signal-to-noise ratio threshold for source extraction. + aperture_radius : float, optional Radius of the circular aperture in pixels for flux summation (default is 7.0). + Returns ------- source_list : pd.DataFrame Table of detected sources with aperture photometry columns. + image_sub : np.ndarray Background-subtracted image array. + """ bkg = sep.Background(image) image_sub = image - bkg.back() @@ -185,15 +192,18 @@ def load_wcs(output_wcs): """ Load a WCS solution from a FITS file header. + Parameters ---------- output_wcs : str Path to the FITS file containing the WCS header from astrometry.net(). + Returns ------- wcs_solution : astropy.wcs.WCS World coordinate system solution object. + """ if not os.path.exists(output_wcs): raise FileNotFoundError(f"WCS file not found: {output_wcs}") @@ -231,19 +241,26 @@ def calibrate_photometry(sky_coords, source_list, phot_cfg): """ Calibrate instrumental magnitudes against a Pan-STARRS1 catalog. + Parameters ---------- sky_coords : astropy.coordinates.SkyCoord Celestial coordinates of detected sources. + source_list : pd.DataFrame Table of detected sources containing 'aperture_sum'. + catalog : str, optional Name of the photometric catalog class in calviacat (default 'PanSTARRS1'). + obs_band : str, optional - Filter of the observed image (used for labeling and color index only; default: 'obs_band'). + Filter of the observed image (used for labeling and color index only; + default: 'obs_band'). + cal_band : str, optional Reference catalog filter for color term (e.g. 'g', 'r', 'i'; default 'g'). + Returns ------- calibration : dict @@ -259,6 +276,7 @@ def calibrate_photometry(sky_coords, source_list, phot_cfg): - color_index : str, the color string used (e.g. 'r-g') - objids : array_like, matched catalog object IDs - distances : array_like, matching distances + """ catalog = phot_cfg["catalog"] obs_band = phot_cfg["obs_band"] @@ -511,7 +529,6 @@ def run_pipeline(input_fits: str, user_config: dict) -> dict: run_solve_field(input_fits, output_wcs, wcs_cfg) wcs_solution = load_wcs(output_wcs) - header = wcs_solution.to_header() ny, nx = image.shape center_world = wcs_solution.pixel_to_world(nx / 2.0, ny / 2.0) diff --git a/catch_analysis_tools/app/services/photometry.py b/catch_analysis_tools/app/services/photometry.py index 1245266..f48bdd6 100644 --- a/catch_analysis_tools/app/services/photometry.py +++ b/catch_analysis_tools/app/services/photometry.py @@ -19,7 +19,8 @@ def get_world_coordinates(WCS_file, x, y): """ - Accepts an astropy-readable WCS object (a .wcs or fits file ) and outputs the world coordinates of a 0-indexed (x,y) pixel point in decimal degrees. + Accepts an astropy-readable WCS object (a .wcs or fits file ) and outputs + the world coordinates of a 0-indexed (x,y) pixel point in decimal degrees. """ try: world_coords = WCS(fits.open(WCS_file)[0].header) @@ -31,7 +32,8 @@ def get_world_coordinates(WCS_file, x, y): world_coords = WCS_file except Exception: raise ValueError( - "Could not read WCS file. Please ensure that the file is either a .wcs file or a .fits file with a valid WCS solution in the header." + "Could not read WCS file. Please ensure that the file is either a " + ".wcs file or a .fits file with a valid WCS solution in the header." ) loc = world_coords.pixel_to_world(x, y) @@ -42,7 +44,8 @@ def get_world_coordinates(WCS_file, x, y): def get_pixel_coordinates(WCS_file, ra, dec): """ - Accepts an astropy-readable WCS object (a .wcs or fits file) and outputs the 0-indexed (x,y) pixel coordinates of an (ra,dec) point in decimal degrees. + Accepts an astropy-readable WCS object (a .wcs or fits file) and outputs the + 0-indexed (x,y) pixel coordinates of an (ra,dec) point in decimal degrees. """ try: world_coords = WCS(fits.open(WCS_file)[0].header) @@ -54,7 +57,8 @@ def get_pixel_coordinates(WCS_file, ra, dec): world_coords = WCS_file except Exception: raise ValueError( - "Could not read WCS file. Please ensure that the file is either a .wcs file or a .fits file with a valid WCS solution in the header." + "Could not read WCS file. Please ensure that the file is either a " + ".wcs file or a .fits file with a valid WCS solution in the header." ) sky_loc = SkyCoord(ICRS(ra=ra * u.deg, dec=dec * u.deg)) @@ -68,19 +72,23 @@ def get_pixel_coordinates(WCS_file, ra, dec): def centroid(file, target_x, target_y, search_radius): """ - Searches for source nearby to expected ephemeris location in image, assumes that astrometry solution has been rerun and that both the cutout .fits file and the redone .wcs file exist. + Searches for source nearby to expected ephemeris location in image, assumes + that astrometry solution has been rerun and that both the cutout .fits file + and the redone .wcs file exist. Parameters ---------- file : string - Base filename (without .fits or .wcs extension) taken from CATCH-generated query URL. + Base filename (without .fits or .wcs extension) taken from + CATCH-generated query URL. target_x : float Target x pixel location, to be used as initial guess for the object. target_y : float - Target y pixel ephemeris location, to be used as initial guess for the object. + Target y pixel ephemeris location, to be used as initial guess for the + object. search_radius : float Radius in pixels of search area (centered on (x,y) ) for centroiding @@ -88,13 +96,12 @@ def centroid(file, target_x, target_y, search_radius): Returns ------- - search_results : - array_like + search_results : array_like + figname : string + Output plot showing the default aperture + annulus extraction onto the + cutout image. - figname: - string - Output plot showing the default aperture + annulus extraction onto the cutout image. """ img, header = get_image(file) @@ -137,19 +144,24 @@ def centroid(file, target_x, target_y, search_radius): return centroid_results -# todo: create function that takes target location, background location and outputs aperture objects -# this function can be fed the outputs of centroid_location +# todo: create function that takes target location, background location and +# outputs aperture objects this function can be fed the outputs of +# centroid_location def target_extraction(body): """ - Searches for source nearby to expected ephemeris location in image, assumes that astrometry solution has been rerun and that both the cutout .fits file and the redone .wcs file exist. + Searches for source nearby to expected ephemeris location in image, assumes + that astrometry solution has been rerun and that both the cutout .fits file + and the redone .wcs file exist. Parameters ---------- body : dict - Request body containing file, target_aperture_params, and background_aperture_params. + Request body containing file, target_aperture_params, and + background_aperture_params. + """ file = body["file"] target_aperture_params = body["target_aperture_params"] @@ -166,8 +178,8 @@ def target_extraction(body): img, target_aperture, background_aperture ) - # could put code to filter out frames where targ_loc and targ_cent vary by more than a couple pixels here - # (would mean star hit) + # could put code to filter out frames where targ_loc and targ_cent vary by + # more than a couple pixels here (would mean star hit) plt.figure(figsize=(8, 8)) plt.imshow(img, norm=norm, cmap="gray_r") diff --git a/catch_analysis_tools/app/services/subtract_median_background.py b/catch_analysis_tools/app/services/subtract_median_background.py index 811b6e6..758efb5 100644 --- a/catch_analysis_tools/app/services/subtract_median_background.py +++ b/catch_analysis_tools/app/services/subtract_median_background.py @@ -9,15 +9,21 @@ def download_file(url): """ - Downloads a file from a given URL and saves it locally, returning the base filename (without extension). + Downloads a file from a given URL and saves it locally, returning the base + filename (without extension). + + Parameters: ---------- url: str URL to the file to be downloaded + + Returns: ------- file_base: str Base filename (without extension) of the downloaded file + """ response = requests.get(url) if "content-disposition" in response.headers: @@ -36,18 +42,21 @@ def download_file(url): def perform_median_subtraction(url): - """Takes a cutout URL from CATCH, opens the file, performs a median background subtraction - after masking sources. Returns background subtracted image. + """ + Takes a cutout URL from CATCH, opens the file, performs a median background + subtraction after masking sources. Returns background subtracted image. Parameters: ---------- url: str URL to CATCH cutout image + Returns: ------- subt_fname: str Filename of background subtracted image saved locally + """ # Download file and get file base name diff --git a/catch_analysis_tools/astrometry.py b/catch_analysis_tools/astrometry.py index 3f1dc00..cd105af 100644 --- a/catch_analysis_tools/astrometry.py +++ b/catch_analysis_tools/astrometry.py @@ -20,25 +20,32 @@ def run_solve_field( """ Execute the `solve-field` command to compute a WCS solution. + Parameters ---------- input_fits : str Path to the input FITS image. + output_wcs : str Path for the output WCS solution file. + pixel_scale : float Approximate pixel scale (e.g., arcsec/pixel). + scale_units : str, optional Units for pixel scale (default is "arcsecperpix"). + Returns ------- success : bool True if the solve-field command succeeded or file already exists. + """ if os.path.exists(output_wcs): print( - f"Output file '{output_wcs}' already exists. Skipping solve-field execution." + f"Output file '{output_wcs}' already exists. " + "Skipping solve-field execution." ) return True @@ -81,24 +88,32 @@ def find_sources(image_sub, bkg_err, snr, aperture_radius=7.0): """ Detect sources in an image using SEP background subtraction and extraction. + Parameters ---------- image_sub : array_like 2D numpy array after background subtraction (cleaned image). + bkg_err : float or array_like Background noise estimate (global RMS or per‐pixel error map). + snr : float Minimum signal-to-noise ratio threshold for source extraction. + aperture_radius : float, optional Radius of the circular aperture in pixels for flux summation (default is 7.0). + Returns ------- source_list : pd.DataFrame Table of detected sources with aperture photometry columns. + image_sub : np.ndarray Background-subtracted image array. + """ + sep.set_sub_object_limit(500) sources = sep.extract(image_sub, thresh=snr, err=bkg_err, deblend_nthresh=16) source_list = pd.DataFrame(sources) @@ -115,16 +130,20 @@ def load_wcs(output_wcs): """ Load a WCS solution from a FITS file header. + Parameters ---------- output_wcs : str Path to the FITS file containing the WCS header from astrometry.net(). + Returns ------- wcs_solution : astropy.wcs.WCS World coordinate system solution object. + """ + if not os.path.exists(output_wcs): raise FileNotFoundError(f"WCS file not found: {output_wcs}") with fits.open(output_wcs) as hdul: @@ -136,20 +155,26 @@ def retrieve_sources(source_list, wcs_solution): """ Convert pixel coordinates to sky coordinates using a WCS. + Parameters ---------- source_list : pd.DataFrame Table with 'x' and 'y' pixel positions of detected sources. + wcs_solution : astropy.wcs.WCS World coordinate system solution object. + Returns ------- source_list : pd.DataFrame Updated table including 'RA' and 'Dec' columns in degrees. + sky_coords : astropy.coordinates.SkyCoord SkyCoord object with celestial coordinates of sources. + """ + world = wcs_solution.pixel_to_world(source_list["x"], source_list["y"]) source_list["RA"] = [c.ra.deg for c in world] source_list["Dec"] = [c.dec.deg for c in world] @@ -167,18 +192,27 @@ def calibrate_photometry( """ Calibrate instrumental magnitudes against a Pan-STARRS1 catalog. + Parameters ---------- sky_coords : astropy.coordinates.SkyCoord Celestial coordinates of detected sources. + source_list : pd.DataFrame Table of detected sources containing 'aperture_sum'. + catalog : str, optional - Name of the photometric catalog class in calviacat (default 'PanSTARRS1'). + Name of the photometric catalog class in calviacat (default + 'PanSTARRS1'). + obs_band : str, optional - Filter of the observed image (used for labeling and color index only; default: 'obs_band'). + Filter of the observed image (used for labeling and color index only; + default: 'obs_band'). + cal_band : str, optional - Reference catalog filter for color term (e.g. 'g', 'r', 'i'; default 'g'). + Reference catalog filter for color term (e.g. 'g', 'r', 'i'; default + 'g'). + Returns ------- @@ -195,7 +229,9 @@ def calibrate_photometry( - color_index : str, the color string used (e.g. 'r-g') - objids : array_like, matched catalog object IDs - distances : array_like, matching distances + """ + color_index = f"{obs_band}-{cal_band}" try: @@ -234,25 +270,33 @@ def plot_color_correction(color_mags, m, m_inst, C, zp, color_index: str): """ Plot the relation between instrumental and calibrated magnitudes. + Parameters ---------- color_mags : array_like Color indices (obs_band - cal_band) of matched stars. + m : array_like Calibrated magnitudes from reference catalog. + m_inst : array_like Instrumental magnitudes measured. + C : float Color term coefficient. + zp : float Zero-point magnitude. + color_index : str Label for the color axis (e.g. 'r-g'). + Returns ------- fig, ax : tuple Matplotlib figure and axis objects for the plot. + """ fig, ax = plt.subplots() ax.scatter(color_mags, m - m_inst, marker=".") @@ -268,21 +312,28 @@ def plot_image(telescope_image_sub, source_list, matched_idx, colored_idx): """ Overlay detected and matched sources on the background-subtracted image. + Parameters ---------- + telescope_image_sub : np.ndarray Background-subtracted image array. + source_list : pd.DataFrame Table of detected sources with 'x' and 'y' pixel positions. + matched_idx : array_like Indices of matched catalog sources in source_list. + colored_idx : array_like Indices of sources selected for color correction. + Returns ------- fig, ax : tuple Matplotlib figure and axis objects for the plot. + """ fig, ax = plt.subplots() m, s = np.mean(telescope_image_sub), np.std(telescope_image_sub) @@ -338,28 +389,38 @@ def create_header( """ Create and write a FITS file with calibrated header and source tables. + Parameters ---------- image : array_like Original 2D image data array. + wcs_solution : astropy.wcs.WCS WCS solution for the image. + zp : float Zero-point magnitude. + unc : float Uncertainty of the zero-point. + source_list : pd.DataFrame Table of detected sources. + matched_idx : array_like Indices for matched catalog sources. + colored_idx : array_like Indices for color-correction sources. + input_fits : str, optional Filename for the input FITS file. + Returns ------- None + """ image_arr = np.asarray(image) primary_hdu = fits.PrimaryHDU(data=image_arr, header=wcs_solution.to_header()) @@ -399,14 +460,17 @@ def cleanup_files(file_base): """ Remove temporary files generated during the processing pipeline. + Parameters ---------- file_base : str Base filename (without extension) for the files to remove. + Returns ------- None + """ extensions = [ ".axy", diff --git a/catch_analysis_tools/background.py b/catch_analysis_tools/background.py index 20944db..076a550 100644 --- a/catch_analysis_tools/background.py +++ b/catch_analysis_tools/background.py @@ -10,8 +10,9 @@ def get_background(data): - """computes and returns a global background subtraction, masking sources using image segmentation IDs. - Does NOT return background subtracted data + """ + Computes and returns a global background subtraction, masking sources using + image segmentation IDs. Does NOT return background subtracted data Parameters @@ -20,16 +21,15 @@ def get_background(data): 2D image array to compute background on - Returns ------- bkg : background object returned from photutils Background2D - """ - # performs a global background subtraction, masking sources using image segmentation IDs + # performs a global background subtraction, masking sources using image + # segmentation IDs sigma_clip = SigmaClip(sigma=3.0, maxiters=10) threshold = detect_threshold(data, nsigma=2.0, sigma_clip=sigma_clip) @@ -52,7 +52,9 @@ def get_background(data): def global_subtraction(data): - """performs a global background subtraction, masking sources using image segmentation IDs + """ + Performs a global background subtraction, masking sources using image + segmentation IDs. Parameters @@ -61,15 +63,14 @@ def global_subtraction(data): 2D image array to be background subtracted - Returns ------- data_sub : array_like Background subtracted data array + bkg : background object returned from get_background(data) - """ bkg = get_background(data) @@ -79,33 +80,36 @@ def global_subtraction(data): def calc_bkg(data, background_aperture, method, sigma_clip): - """Takes in an image array and aperture for background estimate. Determines background - estimator per specified method (mean or median) and variance of pixels, with optional - sigma clipping. + """ + Takes in an image array and aperture for background estimate. Determines + background estimator per specified method (mean or median) and variance of + pixels, with optional sigma clipping. Parameters ---------- data : array_like - image data to be used for photometry + image data to be used for photometry position : array_like - [x,y] pixel location to define center CircularAnnulus object + [x,y] pixel location to define center CircularAnnulus object inner_r : float - Distance (in pixels) from annulus center to the inner edge + Distance (in pixels) from annulus center to the inner edge outer_r : float - Distance (in pixels) from annulus center to the outer edge + Distance (in pixels) from annulus center to the outer edge + Returns ------- bkg_estimator : float - estimate of true background level in background aperture, per "method" + Estimate of true background level in background aperture, per "method" bkg_var : float - Variance of pixels within the background aperture. If sigma clipped, square of - clipped stddev. Otherwise, as square of range of values from 50-16th percentiles + Variance of pixels within the background aperture. If sigma clipped, + square of clipped stddev. Otherwise, as square of range of values from + 50-16th percentiles. """ @@ -114,7 +118,8 @@ def calc_bkg(data, background_aperture, method, sigma_clip): background_data = background_mask.multiply(data) background_data_1d = background_data[background_mask.data > 0] - # Todo: Add flexibility for the stats for computing background. Include rectangular aperture. + # Todo: Add flexibility for the stats for computing background. Include + # rectangular aperture. if sigma_clip is not None: bkg_mean, bkg_median, bkg_stddev = sigma_clipped_stats( background_data_1d, sigma=sigma_clip, maxiters=10 diff --git a/catch_analysis_tools/calibration/astrometry.py b/catch_analysis_tools/calibration/astrometry.py index 1022b55..8e92985 100644 --- a/catch_analysis_tools/calibration/astrometry.py +++ b/catch_analysis_tools/calibration/astrometry.py @@ -18,25 +18,33 @@ def run_solve_field( """ Execute the `solve-field` command to compute a WCS solution. + Parameters ---------- input_fits : str Path to the input FITS image. + output_wcs : str Path for the output WCS solution file. + pixel_scale : float Approximate pixel scale (e.g., arcsec/pixel). + scale_units : str, optional Units for pixel scale (default is "arcsecperpix"). + Returns ------- success : bool True if the solve-field command succeeded or file already exists. + """ + if os.path.exists(output_wcs): print( - f"Output file '{output_wcs}' already exists. Skipping solve-field execution." + f"Output file '{output_wcs}' already exists. " + "Skipping solve-field execution." ) return True @@ -79,24 +87,32 @@ def find_sources(image_sub, bkg_err, snr, aperture_radius=7.0): """ Detect sources in an image using SEP background subtraction and extraction. + Parameters ---------- image_sub : array_like 2D numpy array after background subtraction (cleaned image). + bkg_err : float or array_like Background noise estimate (global RMS or per‐pixel error map). + snr : float Minimum signal-to-noise ratio threshold for source extraction. + aperture_radius : float, optional Radius of the circular aperture in pixels for flux summation (default is 7.0). + Returns ------- source_list : pd.DataFrame Table of detected sources with aperture photometry columns. + image_sub : np.ndarray Background-subtracted image array. + """ + sep.set_sub_object_limit(500) sources = sep.extract(image_sub, thresh=snr, err=bkg_err, deblend_nthresh=16) source_list = pd.DataFrame(sources) @@ -113,16 +129,20 @@ def load_wcs(output_wcs): """ Load a WCS solution from a FITS file header. + Parameters ---------- output_wcs : str Path to the FITS file containing the WCS header from astrometry.net(). + Returns ------- wcs_solution : astropy.wcs.WCS World coordinate system solution object. + """ + if not os.path.exists(output_wcs): raise FileNotFoundError(f"WCS file not found: {output_wcs}") with fits.open(output_wcs) as hdul: @@ -134,20 +154,26 @@ def retrieve_sources(source_list, wcs_solution): """ Convert pixel coordinates to sky coordinates using a WCS. + Parameters ---------- source_list : pd.DataFrame Table with 'x' and 'y' pixel positions of detected sources. + wcs_solution : astropy.wcs.WCS World coordinate system solution object. + Returns ------- source_list : pd.DataFrame Updated table including 'RA' and 'Dec' columns in degrees. + sky_coords : astropy.coordinates.SkyCoord SkyCoord object with celestial coordinates of sources. + """ + world = wcs_solution.pixel_to_world(source_list["x"], source_list["y"]) source_list["RA"] = [c.ra.deg for c in world] source_list["Dec"] = [c.dec.deg for c in world] @@ -159,15 +185,19 @@ def cleanup_files(file_base): """ Remove temporary files generated during the processing pipeline. + Parameters ---------- file_base : str Base filename (without extension) for the files to remove. + Returns ------- None + """ + extensions = [ ".axy", ".corr", @@ -195,10 +225,12 @@ def write_astrometry_output( """ Write an astrometrically calibrated FITS file. - This writes WCS information and detected source coordinates only. - It intentionally does not write photometric calibration metadata such as - zero point, color term, reference catalog, or reference filter. + This writes WCS information and detected source coordinates only. It + intentionally does not write photometric calibration metadata such as zero + point, color term, reference catalog, or reference filter. + """ + image_arr = np.asarray(image) primary_hdu = fits.PrimaryHDU( @@ -231,6 +263,7 @@ def run_astrometry_calibration( """ Run astrometric calibration only. """ + file_base = os.path.splitext(input_fits)[0] if output_fits is None: diff --git a/catch_analysis_tools/photometry.py b/catch_analysis_tools/photometry.py index e960759..4ab4e0a 100644 --- a/catch_analysis_tools/photometry.py +++ b/catch_analysis_tools/photometry.py @@ -16,7 +16,8 @@ from photutils.utils import circular_footprint -# here are functions for grabbing the data, doing background subtractions and manipulating source extractions +# here are functions for grabbing the data, doing background subtractions and +# manipulating source extractions def get_image(url): """Access a cutout via a call to a CATCH URL @@ -27,15 +28,14 @@ def get_image(url): CATCH-generated URL from a query. - Returns ------- data : array_like This is a 2D array containing the image data returned by the CATCH query + header : FITS header class, from astropy.io.fits.Header - """ fits_hdu = fits.open(url) @@ -45,10 +45,13 @@ def get_image(url): def id_good_sources(data, bkg): - """Uses a segmentation image to identify reliable sources in image that can be snapped to. + """ + Uses a segmentation image to identify reliable sources in image that can be + snapped to. - Coincidentally, computes baseline photometry that could be used as a quality comparison user results, - though this flux isn't always a good comparison as it often underestimates the source size + Coincidentally, computes baseline photometry that could be used as a quality + comparison user results, though this flux isn't always a good comparison as + it often underestimates the source size Parameters @@ -59,12 +62,12 @@ def id_good_sources(data, bkg): bkg : background object returned from get_background() or global_subtraction() + Returns ------- cat : - Astropy Table class, from SourceCatalog output giving source locations, fluxes - - + Astropy Table class, from SourceCatalog output giving source locations, + fluxes """ @@ -75,8 +78,8 @@ def id_good_sources(data, bkg): finder = SourceFinder(npixels=5, progress_bar=False) segment_map = finder(convolved_data, source_threshold) - vmax = np.percentile(np.ndarray.flatten(data), 99) - vmin = np.percentile(np.ndarray.flatten(data), 1) + # vmax = np.percentile(np.ndarray.flatten(data), 99) + # vmin = np.percentile(np.ndarray.flatten(data), 1) # make a plot to show the background subtracted frame and the resulting segment map # fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 12.5)) @@ -98,18 +101,18 @@ def create_user_aperture(position, size): Parameters ---------- - Position : array_like + position : array_like [x,y] location of aperture center - Size : Int + size : Int Radius of aperture, given in pixels + Returns ------- aperture : - Photutils circular aperture object with the specified size, location parameters - - + Photutils circular aperture object with the specified size, location + parameters """ @@ -119,8 +122,10 @@ def create_user_aperture(position, size): def define_aperture(aperture_params): """ - Defines an aperture from user supplied parameters, used both to create target and background apertures + Defines an aperture from user supplied parameters, used both to create + target and background apertures """ + if aperture_params["shape"] == "Circular": aperture = CircularAperture( aperture_params["position"], r=aperture_params["size"] @@ -144,7 +149,9 @@ def define_aperture(aperture_params): def subpixel_centroid(user_point, data, radius): - """Takes in a user defined point and returns the location of the brightest pixel within radius pixels + """ + Takes in a user defined point and returns the location of the brightest + pixel within radius pixels. Parameters @@ -153,16 +160,16 @@ def subpixel_centroid(user_point, data, radius): [x,y] location of desired point, to be snapped from data : array_like - image to be searched for sources to be snapped to + image to be searched for sources to be snapped to radius : Int - number of pixels within which objects can be snapped to from user_point + number of pixels within which objects can be snapped to from user_point + Returns ------- target_position: array_like - [x,y] location of the source that is closest to user_point found in data - + [x,y] location of the source that is closest to user_point found in data """ @@ -180,39 +187,44 @@ def subpixel_centroid(user_point, data, radius): def do_aperture_photometry(data, source_aperture, bkg_aperture): - """Takes in an image, a source aperture, and outputs from the calc_annulus_bkg function + """ + Takes in an image, a source aperture, and outputs from the calc_annulus_bkg + function. - Returns the source flux (background subtracted, per-pixel background median) and the - uncertainty as defined at the quoted link + Returns the source flux (background subtracted, per-pixel background median) + and the uncertainty as defined at the quoted link + + method='center' means pixels are either in or out, no interpolation to a + perfect circle (in other words, areas will be in whole pixels) - method='center' means pixels are either in or out, no interpolation to a perfect circle - (in other words, areas will be in whole pixels) Parameters ---------- data : array_like - image data to be used for photometry + image data to be used for photometry source_aperture : - Photutils aperture object containing desired source + Photutils aperture object containing desired source bkg_median : float - median value of pixel background, + median value of pixel background, bkg_var : float - variance of pixel background, ideally from output of calc_annulus_background + variance of pixel background, ideally from output of + calc_annulus_background bkg_aperture : - Photutils aperture object containing background, ideally as annulus_aperture output from calc_annulus_background + Photutils aperture object containing background, ideally as + annulus_aperture output from calc_annulus_background Returns ------- source_sum : float - background subtracted flux of the targeted source + background subtracted flux of the targeted source source_err : float - error on source flux + error on source flux """ @@ -251,17 +263,16 @@ def load_thumbnail(url): CATCH-generated URL from a query. - Returns ------- data : array_like This is a 2D array containing the image data returned by the CATCH query header : - FITS header class, from astropy.io.fits.Header + FITS header class, from astropy.io.fits.Header img_WCS: - Astropy WCS object + Astropy WCS object """ fits_hdu = fits.open(url) @@ -273,29 +284,28 @@ def load_thumbnail(url): def source_instr_mag(ap_flux, ap_fluxerr, exposure_time): """Quick function to return instrumental magnitudes from a source flux - Does not force magnitude uncertainties to be symmetric + Does not force magnitude uncertainties to be symmetric - To be used by calibrated_mag() function + To be used by calibrated_mag() function Parameters ---------- - ap_flux: float - Flux (in counts) of source + Flux (in counts) of source ap_fluxerr: float - Flux error (in counts) + Flux error (in counts) exposure_time: float - integration time of the frame (s) + integration time of the frame (s) Returns ------- - instr_mag_array: array_like - array containing source instrumental magnitude and uncertainties, - as [Magnitude, Upper Magnitude Uncertainty, Lower Magnitude Uncertainty] + array containing source instrumental magnitude and uncertainties, as + [Magnitude, Upper Magnitude Uncertainty, Lower Magnitude Uncertainty] + """ instr_mag = -2.5 * np.log10(ap_flux / exposure_time) @@ -311,28 +321,34 @@ def source_instr_mag(ap_flux, ap_fluxerr, exposure_time): def calibrated_mag(instr_mag_array, zero_point, zero_point_uncert): - """Takes in the array from source_instr_mag, converts to derived magnitude, - propagating uncertainties from both + """ + Takes in the array from source_instr_mag, converts to derived magnitude, + propagating uncertainties from both Parameters ---------- instr_mag_array: array_like - From source_instr_mag() output: array containing source instrumental magnitude and uncertainties, - as [Magnitude, Upper Magnitude Uncertainty, Lower Magnitude Uncertainty] + From source_instr_mag() output: array containing source instrumental + magnitude and uncertainties, as [Magnitude, Upper Magnitude Uncertainty, + Lower Magnitude Uncertainty] zero_point: float - zero point magnitude of image (ideally taken from metadata in header given by astrometry solution) + zero point magnitude of image (ideally taken from metadata in header + given by astrometry solution) zero_point_uncert: float - zero point uncertainty (mags) of image (ideally taken from metadata in header given by astrometry solution) + zero point uncertainty (mags) of image (ideally taken from metadata in + header given by astrometry solution) + Returns ------- calib_mag_array: array_like - array containing source CALIBRATED magnitude and uncertainties, - as [Magnitude, Upper Magnitude Uncertainty, Lower Magnitude Uncertainty] + array containing source CALIBRATED magnitude and uncertainties, as + [Magnitude, Upper Magnitude Uncertainty, Lower Magnitude Uncertainty] + """ calib_mag = zero_point + instr_mag_array[0] diff --git a/catch_analysis_tools/tests/astrometry_test.py b/catch_analysis_tools/tests/astrometry_test.py index cc2affe..6e807c9 100644 --- a/catch_analysis_tools/tests/astrometry_test.py +++ b/catch_analysis_tools/tests/astrometry_test.py @@ -1,11 +1,29 @@ +import os import shutil +import subprocess import tempfile from pathlib import Path from unittest.mock import MagicMock, patch +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd import pytest - -from ..astrometry import * +from astropy.coordinates import SkyCoord +from astropy.io import fits +from astropy.wcs import WCS + +from ..astrometry import ( + calibrate_photometry, + cleanup_files, + create_header, + find_sources, + load_wcs, + plot_color_correction, + plot_image, + retrieve_sources, + run_solve_field, +) RA_DEG = 263.0 DEC_DEG = 34.5 diff --git a/catch_analysis_tools/tests/background_test.py b/catch_analysis_tools/tests/background_test.py index f04a485..41a05ad 100644 --- a/catch_analysis_tools/tests/background_test.py +++ b/catch_analysis_tools/tests/background_test.py @@ -1,7 +1,11 @@ import numpy as np from pytest import approx -from ..background import * +from ..background import ( + calc_bkg, + get_background, + global_subtraction, +) from ..photometry import define_aperture @@ -10,7 +14,6 @@ def test_global_subtraction(): data[40:50, 40:50] = 5 * data[40:50, 40:50] data_sub, bkg = global_subtraction(data) - errors = [] # test some basic values checked manually against above image values assert np.mean(data_sub) == 0.04 assert np.mean(bkg.background) == 1.0 diff --git a/catch_analysis_tools/tests/photometry_test.py b/catch_analysis_tools/tests/photometry_test.py index a9550ff..4872ae2 100644 --- a/catch_analysis_tools/tests/photometry_test.py +++ b/catch_analysis_tools/tests/photometry_test.py @@ -3,11 +3,24 @@ from photutils.datasets import make_wcs from pytest import approx -from ..app.services.photometry import * -from ..background import * - -# here are test functions for grabbing the data, doing background subtractions and manipulating source extractions -from ..photometry import * +from ..app.services.photometry import ( + get_image, + get_pixel_coordinates, + get_world_coordinates, +) +from ..background import get_background + +# here are test functions for grabbing the data, doing background subtractions +# and manipulating source extractions +from ..photometry import ( + calibrated_mag, + create_user_aperture, + define_aperture, + do_aperture_photometry, + id_good_sources, + source_instr_mag, + subpixel_centroid, +) @pytest.mark.remote_data diff --git a/pyproject.toml b/pyproject.toml index f2a2182..c540131 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -44,6 +44,7 @@ packages = ["catch_analysis_tools"] write_to = "catch_analysis_tools/version.py" [tool.ruff] +line-length = 88 lint.select = [ "E", # pycodestyle errors "W", # pycodestyle warnings diff --git a/scripts/Checking_indexfiles.py b/scripts/Checking_indexfiles.py index 76ed334..5afb715 100644 --- a/scripts/Checking_indexfiles.py +++ b/scripts/Checking_indexfiles.py @@ -1,4 +1,4 @@ -from catch_analysis_tools.app.astrometry_readiness.get_astrometry_readiness_status import ( +from catch_analysis_tools.app.astrometry_readiness.get_astrometry_readiness_status import ( # noqa: E501 get_astrometry_readiness_status, ) from catch_analysis_tools.app.astrometry_readiness.prepare_astrometry_data import (