From 5d43d7e7856f697ff850393d380a4671cd9a8cff Mon Sep 17 00:00:00 2001 From: yaeji43 Date: Wed, 20 May 2026 16:55:20 -0400 Subject: [PATCH 1/2] Separate astrometry and photometry functionality into modular files --- .../calibration/astrometry.py | 323 ++++++++++++++++++ .../calibration/photometry.py | 289 ++++++++++++++++ 2 files changed, 612 insertions(+) create mode 100644 catch_analysis_tools/calibration/astrometry.py create mode 100644 catch_analysis_tools/calibration/photometry.py diff --git a/catch_analysis_tools/calibration/astrometry.py b/catch_analysis_tools/calibration/astrometry.py new file mode 100644 index 0000000..d8f8dca --- /dev/null +++ b/catch_analysis_tools/calibration/astrometry.py @@ -0,0 +1,323 @@ +import os +import subprocess +import argparse +import numpy as np +import pandas as pd +import sep +import fitsio +from astropy.io import fits +from astropy.wcs import WCS +from astropy.table import Table +from astropy.coordinates import SkyCoord + + +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. + + 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.") + 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." + ) + + 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", + input_fits, + ] + + try: + subprocess.run(command, check=True) + return True + except subprocess.CalledProcessError as e: + raise RuntimeError(f"solve-field failed: {e}") + + +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) + flux, flux_err, _ = sep.sum_circle( + 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) + return source_list, image_sub + + +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: + wcs_solution = WCS(hdul[0].header) + return wcs_solution + + +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] + sky_coords = SkyCoord(source_list['RA'], source_list['Dec'], unit='deg') + return source_list, sky_coords + + +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', '.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, + source_list, + output_fits, +): + """ + 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. + """ + image_arr = np.asarray(image) + + primary_hdu = fits.PrimaryHDU( + data=image_arr, + header=wcs_solution.to_header(), + ) + + 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", + ) + + hdul = fits.HDUList([primary_hdu, detected_hdu]) + hdul.writeto(output_fits, overwrite=True) + + +def run_astrometry_calibration( + input_fits, + ra_deg, + dec_deg, + bkg_err, + pixel_scale=1.86, + snr=7.0, + output_fits=None, +): + """ + Run astrometric calibration only. + """ + file_base = os.path.splitext(input_fits)[0] + + if output_fits is None: + output_fits = f"{file_base}_astrometry.fits" + + output_wcs = f"{file_base}.wcs" + + image = fitsio.read(input_fits).astype(np.float32) + + if run_solve_field(input_fits, output_wcs, pixel_scale, ra_deg, dec_deg): + wcs_solution = load_wcs(output_wcs) + else: + raise RuntimeError("solve-field did not produce a WCS solution.") + + source_list, image_sub = find_sources(image, bkg_err, snr) + source_list, sky_coords = retrieve_sources(source_list, wcs_solution) + + write_astrometry_output( + image=image, + wcs_solution=wcs_solution, + source_list=source_list, + output_fits=output_fits, + ) + + cleanup_files(file_base) + + return { + "output_fits": output_fits, + "output_wcs": output_wcs, + "source_list": source_list, + "sky_coords": sky_coords, + "wcs_solution": wcs_solution, + } + + +if __name__ == "__main__": + parser = argparse.ArgumentParser( + description="Astrometric calibration for a FITS image." + ) + + parser.add_argument( + "input_fits", + help="Path to input FITS image.", + ) + + parser.add_argument( + "--Ra", + type=float, + required=True, + help="Initial RA estimate in degrees.", + ) + + parser.add_argument( + "--Dec", + type=float, + required=True, + help="Initial Dec estimate in degrees.", + ) + + parser.add_argument( + "--bkg_err", + type=float, + required=True, + help="Global background RMS.", + ) + + 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( + "--output_fits", + default=None, + help="Output astrometrically calibrated FITS file.", + ) + + args = parser.parse_args() + + try: + result = run_astrometry_calibration( + input_fits=args.input_fits, + ra_deg=args.Ra, + dec_deg=args.Dec, + bkg_err=args.bkg_err, + pixel_scale=args.pixel_scale, + snr=args.snr, + output_fits=args.output_fits, + ) + + print(f"Astrometric calibration complete: {result['output_fits']}") + + except Exception as exc: + raise SystemExit(f"Astrometric calibration failed: {exc}") diff --git a/catch_analysis_tools/calibration/photometry.py b/catch_analysis_tools/calibration/photometry.py new file mode 100644 index 0000000..6ccb709 --- /dev/null +++ b/catch_analysis_tools/calibration/photometry.py @@ -0,0 +1,289 @@ +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 + + +def calibrate_photometric_zero_point( + sky_coords, + source_list: pd.DataFrame, + catalog: str = "PanSTARRS1", + obs_band: str = "obs_band", + cal_band: str = "g", + catalog_db: str = "cat.db", +): + """ + Calibrate instrumental magnitudes using a reference catalog. + + This performs photometric calibration only: + catalog matching, instrumental magnitude calculation, + zero-point estimation, and color-term correction. + + Parameters + ---------- + sky_coords : astropy.coordinates.SkyCoord + Sky coordinates of detected sources. + + source_list : pandas.DataFrame + Source table containing aperture fluxes in the + ``aperture_sum`` column. + + catalog : str, optional + Name of the reference catalog class in calviacat. + Default is ``"PanSTARRS1"``. + + obs_band : str, optional + Observed image band label. Used to construct the color index. + + cal_band : str, optional + Reference catalog band used for calibration. + + catalog_db : str, optional + Local catalog database path. + + Returns + ------- + dict + Photometric calibration results containing zero point, + color term, calibrated magnitudes, matched object IDs, + and match distances. + """ + if "aperture_sum" not in source_list.columns: + raise ValueError("source_list must contain an 'aperture_sum' column.") + + aperture_sum = source_list["aperture_sum"].to_numpy() + + if np.any(aperture_sum <= 0): + raise ValueError( + "All aperture_sum values must be positive before magnitude calibration." + ) + + color_index = f"{obs_band}-{cal_band}" + + try: + CatalogClass = getattr(cvc, catalog) + except AttributeError as exc: + raise ValueError(f"Catalog '{catalog}' not found in calviacat.") from exc + + ref = CatalogClass(catalog_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(aperture_sum) + + zp, color_term, zp_unc, m_cal, color_mags, _ = ref.cal_color( + objids, + m_inst, + cal_band, + color_index, + ) + + return { + "zp": zp, + "color_term": color_term, + "unc": zp_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 get_matched_indices(objids, source_count: int): + """ + Return indices of sources matched to catalog objects. + """ + if hasattr(objids, "mask"): + return np.where(~objids.mask)[0] + + return np.arange(source_count) + + +def get_color_corrected_indices(color_mags, source_count: int): + """ + Return indices of sources used for color correction. + """ + if hasattr(color_mags, "mask"): + return np.where(~color_mags.mask)[0] + + return np.arange(source_count) + + +def plot_color_correction( + color_mags, + calibrated_magnitude, + instrumental_magnitude, + color_term, + zero_point, + color_index: str, +): + """ + Plot color correction relation for photometric calibration. + """ + fig, ax = plt.subplots() + + ax.scatter( + color_mags, + calibrated_magnitude - instrumental_magnitude, + marker=".", + ) + + x = np.linspace(0, 1.5, 100) + ax.plot( + x, + color_term * x + zero_point, + 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.legend() + + plt.tight_layout() + + return fig, ax + + +def plot_photometric_matches( + image_sub, + source_list: pd.DataFrame, + matched_idx, + color_corrected_idx, +): + """ + Overlay detected, catalog-matched, and color-corrected sources. + """ + fig, ax = plt.subplots() + + mean_value = np.mean(image_sub) + std_value = np.std(image_sub) + + im = ax.imshow( + image_sub, + interpolation="nearest", + origin="lower", + cmap="gray", + ) + im.set_clim(vmin=mean_value - std_value, vmax=mean_value + std_value) + + 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[color_corrected_idx], + source_list["y"].iloc[color_corrected_idx], + "o", + markersize=15, + color="yellow", + markerfacecolor="none", + label="Selected for Color Correction", + ) + + ax.legend() + + return fig, ax + +def write_photometric_calibration_output( + image, + wcs_solution, + source_list: pd.DataFrame, + matched_idx, + color_corrected_idx, + output_fits, + zero_point, + zero_point_uncertainty, + catalog: str, + obs_band: str, + cal_band: str, + color_index: str, + color_term=None, +): + """ + Write a FITS file with photometric calibration metadata and source tables. + + This should be called after astrometry has already produced a WCS solution + and photometric calibration has estimated the zero point and color term. + """ + image_arr = np.asarray(image) + + primary_hdu = fits.PrimaryHDU( + data=image_arr, + header=wcs_solution.to_header(), + ) + + primary_hdu.header["ZP"] = zero_point + primary_hdu.header["ZP_STD"] = zero_point_uncertainty + primary_hdu.header["REF_CATA"] = catalog + primary_hdu.header["OBS_FLT"] = obs_band + primary_hdu.header["REF_FLT"] = cal_band + primary_hdu.header["CAT_COR"] = color_index + + if color_term is not None: + primary_hdu.header["COLORTRM"] = color_term + + 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", + ) + + color_hdu = fits.BinTableHDU( + Table.from_pandas( + source_list_clean.iloc[color_corrected_idx].reset_index(drop=True) + ), + name="COLOR_CORRECTION_STARS", + ) + else: + matched_hdu = fits.BinTableHDU(name="SELECTED_STARS") + color_hdu = fits.BinTableHDU(name="COLOR_CORRECTION_STARS") + + hdul = fits.HDUList( + [ + primary_hdu, + detected_hdu, + matched_hdu, + color_hdu, + ] + ) + + hdul.writeto(output_fits, overwrite=True) \ No newline at end of file From 2f0d3b2327379be2d0b03e60886193c893f5b7c0 Mon Sep 17 00:00:00 2001 From: yaeji43 Date: Fri, 26 Jun 2026 18:24:13 -0400 Subject: [PATCH 2/2] Split astometry calibration API into astrometry and photometry endpoints under service/calibration forlder --- .gitignore | 1 + catch_analysis_tools/app/api/openapi.yaml | 167 ++---- .../app/services/astrometry.py | 526 +++--------------- .../app/services/calibration/__init__.py | 0 .../app/services/calibration/astrometry.py | 147 +++++ .../app/services/calibration/photometry.py | 426 ++++++++++++++ catch_analysis_tools/astrometry.py | 474 ---------------- .../calibration/astrometry.py | 100 +++- .../calibration/photometry.py | 27 +- 9 files changed, 821 insertions(+), 1047 deletions(-) create mode 100644 catch_analysis_tools/app/services/calibration/__init__.py create mode 100644 catch_analysis_tools/app/services/calibration/astrometry.py create mode 100644 catch_analysis_tools/app/services/calibration/photometry.py delete mode 100644 catch_analysis_tools/astrometry.py diff --git a/.gitignore b/.gitignore index 52c8dbf..5e1f421 100644 --- a/.gitignore +++ b/.gitignore @@ -185,3 +185,4 @@ terraform.tfstate.* Thumbs.db catch_analysis_tools/version.py +.venv/ diff --git a/catch_analysis_tools/app/api/openapi.yaml b/catch_analysis_tools/app/api/openapi.yaml index 268899b..5ac7ca9 100644 --- a/catch_analysis_tools/app/api/openapi.yaml +++ b/catch_analysis_tools/app/api/openapi.yaml @@ -45,18 +45,15 @@ paths: application/json: schema: type: object + /astrometry: post: tags: - Astrometry - summary: Perform astrometry and optionally return a plot image. - operationId: catch_analysis_tools.app.handlers.astrometry.do_astrometry - + summary: Perform geometric plate solving and update FITS headers with WCS data. + operationId: catch_analysis_tools.app.services.calibration.astrometry.run_astrometry_pipeline requestBody: required: true - description: > - Provide parameters for astrometry solving, source detection, - and photometric calibration. content: application/json: schema: @@ -66,110 +63,94 @@ paths: properties: image_url: type: string - description: > - URL of FITS image OR local file path (if supported). example: "https://uxzqjwo0ye.execute-api.us-west-1.amazonaws.com/api/images/urn:nasa:pds:gbo.ast.loneos.survey:data_augmented:060104_1a_099_fits?ra=51.11064&dec=17.38378&size=12.60arcmin&format=fits" - ra: type: number - nullable: true - default: 0.0 - minimum: 0 - maximum: 360 - description: Right Ascension (deg). Used only if use_ra_dec=true. - + default: 51.11064 dec: type: number - nullable: true - default: 0.0 - minimum: -90 - maximum: 90 - description: Declination (deg). Used only if use_ra_dec=true. - + default: 17.38378 use_ra_dec: type: boolean default: false - description: > - If true, constrain astrometry solve around (ra, dec). - If false, perform blind solve. - pixel_scale: type: number - nullable: true default: 2.5 - description: Pixel scale (arcsec/pixel). Used if scale_low/high not provided. - scale_low: type: number nullable: true default: null - description: Lower bound of pixel scale (arcsec/pixel). - scale_high: type: number nullable: true default: null - description: Upper bound of pixel scale (arcsec/pixel). - search_radius: type: number default: 2.0 - minimum: 0 - description: Search radius (degrees) when using RA/Dec constraint. - snr_threshold: type: number default: 3.0 - minimum: 0 - description: Signal-to-noise threshold for source detection. + responses: + "200": + description: WCS geometric calculation successful. Returns updated FITS location. + content: + application/json: + schema: + type: object + properties: + wcs_image_url: + type: string + description: Filepath or URL to the updated FITS file containing WCS headers. + sources_detected: + type: integer + center_ra_deg: + type: number + center_dec_deg: + type: number + pixel_scale: + type: number + "422": + description: Astrometry engine could not solve the submitted field image. + /photometry: + post: + tags: + - Photometry + summary: Perform zero-point flux calibration from an existing WCS-solved FITS file. + operationId: catch_analysis_tools.app.services.calibration.photometry.run_photometry_pipeline + requestBody: + required: true + content: + application/json: + schema: + type: object + required: + - wcs_image_url + properties: + wcs_image_url: + type: string + example: "/var/folders/p7/l1v8bpgn10g_vgyg0kzfv6g00000gn/T/catch_astrometry_m9q8renb/input_astrometry.fits" + description: The WCS-solved FITS path returned from the /astrometry endpoint. aperture_radius: type: number default: 7.0 - minimum: 0 - description: Aperture radius (pixels) for photometry. - catalog: type: string enum: [PanSTARRS1, Gaia, SDSS] default: PanSTARRS1 - description: Reference catalog for photometric calibration. - - obs_band: - type: string - default: g - description: Observed band (label only). - cal_band: type: string default: r - description: Calibration band from catalog. - + color_term: + type: string + default: "g-r" return_plot: type: boolean default: false - description: If true, return PNG instead of JSON. - plot_type: type: string enum: [color_correction, image_overlay] default: color_correction - description: Plot type to return when return_plot=true. - - example: - image_url: "https://uxzqjwo0ye.execute-api.us-west-1.amazonaws.com/api/images/urn:nasa:pds:gbo.ast.loneos.survey:data_augmented:060104_1a_099_fits?ra=51.11064&dec=17.38378&size=12.60arcmin&format=fits" - ra: 51.11064 - dec: 17.38378 - 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 - responses: "200": description: JSON results or PNG plot image. @@ -187,28 +168,14 @@ paths: type: number uncertainty: type: number - sources: type: object properties: - detected: - type: integer matched: type: integer - - astrometry: - type: object - properties: - center_ra_deg: - type: number - center_dec_deg: - type: number - pixel_scale: - type: number - plots: type: object - description: Base64 encoded PNG plots (optional) + description: Base64 encoded PNG plots, optional for JSON responses. properties: color_correction: type: string @@ -220,29 +187,10 @@ paths: type: string format: binary - "400": - description: Bad request (invalid input or FITS error) - - "422": - description: Astrometry could not solve the submitted image. - content: - application/json: - schema: - type: object - - "503": - description: Astrometry index files are not ready yet. - content: - application/json: - schema: - type: object - - "500": - description: Internal server error /subtract_median_background: post: tags: - - Background Estimation + - Background Estimation operationId: catch_analysis_tools.app.services.subtract_median_background.perform_median_subtraction summary: "Compute the median background with source masking, saving the background-subtracted image as a new FITS file." parameters: @@ -265,10 +213,11 @@ paths: background_subtracted_image_path: description: Path to the background-subtracted FITS file. type: string + /get_world_coordinates: post: tags: - - Target Search + - Target Search operationId: catch_analysis_tools.app.services.photometry.get_world_coordinates summary: "Retrieve the RA and Dec in decimal degrees from an (x,y) pixel location, using a supplied WCS file." parameters: @@ -338,10 +287,11 @@ paths: dec: description: Dec value of point in decimal degrees type: number + /get_pixel_coordinates: post: tags: - - Target Search + - Target Search operationId: catch_analysis_tools.app.services.photometry.get_pixel_coordinates summary: "Retrieve the pixel (x,y) coordinates from RA and Dec world coordinates, using a supplied WCS file." parameters: @@ -387,10 +337,11 @@ paths: dec: description: Dec value of point in decimal degrees type: number + /target_photometry: post: tags: - - Photometry + - Photometry operationId: catch_analysis_tools.app.services.photometry.target_extraction summary: "Perform aperture photometry on a target using specified target and background aperture configurations." requestBody: @@ -491,10 +442,11 @@ paths: aperture_figure: description: File path to the output plot showing apertures on the cutout image. type: string + /centroid: post: tags: - - Target Search + - Target Search operationId: catch_analysis_tools.app.services.photometry.centroid summary: "Find nearest target centroid to a user-specified point, with user-specified search radius." parameters: @@ -558,5 +510,4 @@ paths: properties: centroid_figure: description: Location of centroid plot. - type: string - + type: string \ No newline at end of file diff --git a/catch_analysis_tools/app/services/astrometry.py b/catch_analysis_tools/app/services/astrometry.py index 5d826d4..bc679db 100644 --- a/catch_analysis_tools/app/services/astrometry.py +++ b/catch_analysis_tools/app/services/astrometry.py @@ -1,22 +1,24 @@ import numpy as np import os -import subprocess -import pandas as pd -import sep +import io +import base64 import fitsio +import sep import matplotlib matplotlib.use("Agg") -import io -import base64 import matplotlib.pyplot as plt -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 catch_analysis_tools.calibration.astrometry import run_astrometry_calibration +from catch_analysis_tools.calibration.photometry import ( + calibrate_photometric_zero_point, + get_matched_indices, + get_color_corrected_indices, + plot_color_correction, + plot_photometric_matches, + write_photometric_calibration_output +) class AstrometrySolveError(RuntimeError): pass @@ -60,382 +62,7 @@ def merge_config(user_config, default_config): return result -def run_solve_field(input_fits, output_wcs, wcs_cfg): - """ - Execute the `solve-field` command to compute a WCS solution. - """ - 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") - - pixel_scale = wcs_cfg["pixel_scale"] - scale_low = wcs_cfg["scale_low"] or pixel_scale * 0.5 - scale_high = wcs_cfg["scale_high"] or pixel_scale * 2.0 - - 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." - ) - - command = [ - "solve-field", - "--overwrite", - "--config", config_file, - "--fits-image", - "--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"]), - ] - - command += [ - "--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"])), - ] - - command += [ - "--downsample", "1", - input_fits, - ] - - subprocess.run(command, check=True) - if not os.path.exists(output_wcs): - raise AstrometrySolveError( - f"solve-field did not produce a WCS solution: {output_wcs}" - ) - return True - - -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() - - sep.set_sub_object_limit(500) - - sources = sep.extract( - 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'], - det_cfg["aperture_radius"], - err=bkg.globalrms - ) - - 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 - - -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: - wcs_solution = WCS(hdul[0].header, relax=True) - return wcs_solution - - -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] - sky_coords = SkyCoord(source_list['RA'], source_list['Dec'], unit='deg') - return source_list, sky_coords - - -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'). - cal_band : str, optional - Reference catalog filter for color term (e.g. 'g', 'r', 'i'; default 'g'). - - Returns - ------- - calibration : dict - Dictionary with keys: - - zp : float, zero-point magnitude - - C : float, color-term coefficient - - unc : float, uncertainty of zero-point - - m : array_like, calibrated magnitudes in the observed band - - m_inst : array_like, instrumental magnitudes - - obs_band : str, label of the observed band - - cal_band : str, same as input cal_band - - color_mags : array_like, color indices (obs_band - cal_band) - - 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"] - cal_band = phot_cfg["cal_band"] - - color_index = f"{obs_band}-{cal_band}" - - CatalogClass = getattr(cvc, catalog) - 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) - - zp, C, unc, m_cal, color_mags, _ = ref.cal_color( - objids, - m_inst, - cal_band, - 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, - } - - -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='.') - 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)') - plt.tight_layout() - return fig, ax - - -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) - 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.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): - """ - 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()) - 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') - else: - 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) - - -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', '.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 run_pipeline(input_fits: str, user_config: dict) -> dict: - config = merge_config(user_config, DEFAULT_CONFIG) wcs_cfg = config["wcs"] @@ -445,28 +72,43 @@ def run_pipeline(input_fits: str, user_config: dict) -> dict: file_base = os.path.splitext(input_fits)[0] image = fitsio.read(input_fits).astype(np.float32) - output_wcs = f"{file_base}.wcs" - # --- WCS solution --- - run_solve_field(input_fits, output_wcs, wcs_cfg) + # Compute global background RMS needed by your colleagues' updated astrometry code + bkg = sep.Background(image) + bkg_err = bkg.globalrms + + # --- 1. Run Modular Astrometry Calibration --- + astrom_res = run_astrometry_calibration( + input_fits=input_fits, + ra_deg=wcs_cfg.get("ra"), + dec_deg=wcs_cfg.get("dec"), + bkg_err=bkg_err, + pixel_scale=wcs_cfg.get("pixel_scale", 1.86), + snr=det_cfg.get("snr", 7.0), + output_fits=None + ) - wcs_solution = load_wcs(output_wcs) - header = wcs_solution.to_header() + source_list = astrom_res["source_list"] + sky_coords = astrom_res["sky_coords"] + wcs_solution = astrom_res["wcs_solution"] + # Calculate center coordinates using the returned WCS ny, nx = image.shape center_world = wcs_solution.pixel_to_world(nx / 2.0, ny / 2.0) center_ra = float(center_world.ra.deg) center_dec = float(center_world.dec.deg) - # --- Source detection --- - source_list, telescope_image_sub = find_sources(image, det_cfg) - source_list, sky_coords = retrieve_sources(source_list, wcs_solution) - - # --- Photometry --- - calibration = calibrate_photometry(sky_coords, source_list, phot_cfg) + # --- 2. Run Modular Photometry Calibration --- + calibration = calibrate_photometric_zero_point( + sky_coords=sky_coords, + source_list=source_list, + catalog=phot_cfg["catalog"], + obs_band=phot_cfg["obs_band"], + cal_band=phot_cfg["cal_band"] + ) zp = calibration["zp"] - C = calibration["C"] + color_term = calibration["color_term"] unc = calibration["unc"] m = calibration["m"] m_inst = calibration["m_inst"] @@ -474,42 +116,49 @@ def run_pipeline(input_fits: str, user_config: dict) -> dict: color_index = calibration["color_index"] 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)) + # Filter source indices using modular photometry filters + matched_idx = get_matched_indices(objids, len(source_list)) + colored_idx = get_color_corrected_indices(color_mags, len(source_list)) - # --- Plotting --- + # --- 3. Dynamic Plotting Layer --- plots = {} if out_cfg["make_plots"]: - fig1, _ = plot_color_correction(color_mags, m, m_inst, C, zp, color_index) + fig1, _ = plot_color_correction(color_mags, m, m_inst, color_term, zp, color_index) buf1 = io.BytesIO() fig1.savefig(buf1, format="png", dpi=150, bbox_inches="tight") plots["color_correction"] = base64.b64encode(buf1.getvalue()).decode() plt.close(fig1) - fig2, _ = plot_image(telescope_image_sub, source_list, matched_idx, colored_idx) + # Get the background subtracted image array for the stars overlay + image_sub = image - bkg.back() + fig2, _ = plot_photometric_matches(image_sub, source_list, matched_idx, colored_idx) buf2 = io.BytesIO() fig2.savefig(buf2, format="png", dpi=150, bbox_inches="tight") plots["image_overlay"] = base64.b64encode(buf2.getvalue()).decode() plt.close(fig2) - # --- Optional FITS --- + # --- 4. Write Calibrated Output FITS --- if out_cfg["write_fits"]: - create_header( - image, wcs_solution, zp, unc, - source_list, matched_idx, colored_idx, - input_fits, - phot_cfg["obs_band"], - phot_cfg["catalog"], - phot_cfg["cal_band"], + write_photometric_calibration_output( + image=image, + wcs_solution=wcs_solution, + source_list=source_list, + matched_idx=matched_idx, + color_corrected_idx=colored_idx, + output_fits=input_fits, + zero_point=zp, + zero_point_uncertainty=unc, + catalog=phot_cfg["catalog"], + obs_band=phot_cfg["obs_band"], + cal_band=phot_cfg["cal_band"], + color_index=color_index, + color_term=color_term ) - cleanup_files(file_base) - results = { "photometry": { "zero_point": float(zp), - "color_term": float(C), + "color_term": float(color_term), "uncertainty": float(unc), }, "sources": { @@ -534,7 +183,6 @@ def validate_and_normalize(body: Dict[str, Any]) -> Dict[str, Any]: Minimal validation + normalization layer. Keeps API robust without adding heavy dependencies. """ - def get_float(name, default=None, required=False): val = body.get(name, default) if val is None: @@ -549,12 +197,10 @@ def get_float(name, default=None, required=False): def get_bool(name, default=False): return bool(body.get(name, default)) - # --- Required --- image_url = body.get("image_url") if not image_url: raise AstrometryValidationError("image_url is required") - # --- WCS --- use_ra_dec = get_bool("use_ra_dec", False) ra = get_float("ra") dec = get_float("dec") @@ -574,16 +220,13 @@ def get_bool(name, default=False): if scale_low > scale_high: raise AstrometryValidationError("scale_low must be <= scale_high") - # --- Detection --- snr = get_float("snr_threshold", 3.0) aperture_radius = get_float("aperture_radius", 7.0) - # --- Photometry --- catalog = body.get("catalog", "PanSTARRS1") obs_band = body.get("obs_band", "g") cal_band = body.get("cal_band", "r") - # --- Output --- return_plot = get_bool("return_plot", False) plot_type = body.get("plot_type", "color_correction") @@ -618,56 +261,29 @@ def get_bool(name, default=False): } -# ## when testing the code locally if __name__ == "__main__": - - # --- Simulated API request body --- + # --- Local Integration Testing --- body = { - "image_url": "Comet_65P_Gunn_LONEOS.fits", # not used in local test + "image_url": "Comet_65P_Gunn_LONEOS.fits", "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", } - # --- Validate + normalize (same as API path) --- cfg = validate_and_normalize(body) - - # --- Use local FITS file instead of downloading --- input_fits = "Comet_65P_Gunn_LONEOS.fits" - # --- Run pipeline --- - results = run_pipeline( - input_fits=input_fits, - user_config=cfg, - ) - - # --- Print summary --- - print("\n=== RESULTS ===") - print(results) - - # --- Optional: visualize plot locally --- - if cfg["output"]["make_plots"]: - import base64 - import matplotlib.pyplot as plt - import io - - plot_type = cfg["meta"]["plot_type"] - img_bytes = base64.b64decode(results["plots"][plot_type]) - - img = plt.imread(io.BytesIO(img_bytes), format="png") - plt.imshow(img) - plt.axis("off") - plt.title(plot_type) - plt.show() + if os.path.exists(input_fits): + results = run_pipeline(input_fits=input_fits, user_config=cfg) + print("\n=== REFACTORED PIPELINE PIPELINE SUCCESS ===") + print(results) + else: + print(f"Skipping local run execution: '{input_fits}' file not found.") \ No newline at end of file diff --git a/catch_analysis_tools/app/services/calibration/__init__.py b/catch_analysis_tools/app/services/calibration/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/catch_analysis_tools/app/services/calibration/astrometry.py b/catch_analysis_tools/app/services/calibration/astrometry.py new file mode 100644 index 0000000..d6033cc --- /dev/null +++ b/catch_analysis_tools/app/services/calibration/astrometry.py @@ -0,0 +1,147 @@ +import os +from typing import Any, Dict +import fitsio +import numpy as np +import sep +import tempfile +import urllib.request +from urllib.parse import urlparse + + +from catch_analysis_tools.calibration.astrometry import run_astrometry_calibration + + +class AstrometryValidationError(ValueError): + pass + + +class AstrometrySolveError(RuntimeError): + pass + + +def _get_float(body: Dict[str, Any], key: str, default=None): + value = body.get(key, default) + + if value is None: + return None + + try: + return float(value) + except (TypeError, ValueError) as exc: + raise AstrometryValidationError(f"{key} must be a number") from exc + + +def _get_bool(body: Dict[str, Any], key: str, default=False): + value = body.get(key, default) + + if isinstance(value, bool): + return value + + if isinstance(value, str): + return value.lower() in {"true", "1", "yes", "y"} + + return bool(value) + +def materialize_input_fits(image_url: str) -> str: + """ + Convert a local FITS path or remote FITS URL into a clean local FITS path. + """ + parsed = urlparse(image_url) + + if parsed.scheme in {"http", "https", "ftp"}: + tmpdir = tempfile.mkdtemp(prefix="catch_astrometry_") + local_path = os.path.join(tmpdir, "input.fits") + + urllib.request.urlretrieve(image_url, local_path) + return local_path + + return image_url + +def validate_and_normalize_astrometry(body: Dict[str, Any]) -> Dict[str, Any]: + image_url = body.get("image_url") + + if not image_url: + raise AstrometryValidationError("image_url is required") + + use_ra_dec = _get_bool(body, "use_ra_dec", False) + + ra = _get_float(body, "ra") + dec = _get_float(body, "dec") + + if use_ra_dec and (ra is None or dec is None): + raise AstrometryValidationError( + "ra and dec are required when use_ra_dec=True" + ) + + pixel_scale = _get_float(body, "pixel_scale", 2.5) + scale_low = _get_float(body, "scale_low") + scale_high = _get_float(body, "scale_high") + + if pixel_scale is None: + if scale_low is None or scale_high is None: + raise AstrometryValidationError( + "Provide pixel_scale or both scale_low and scale_high" + ) + if scale_low > scale_high: + raise AstrometryValidationError("scale_low must be <= scale_high") + + return { + "image_url": image_url, + "ra": ra, + "dec": dec, + "use_ra_dec": use_ra_dec, + "pixel_scale": pixel_scale, + "scale_low": scale_low, + "scale_high": scale_high, + "search_radius": _get_float(body, "search_radius", 2.0), + "snr_threshold": _get_float(body, "snr_threshold", 3.0), + } + + +def run_astrometry_pipeline(body=None, **kwargs) -> Dict[str, Any]: + """ + OpenAPI/Connexion endpoint for /astrometry. + + Compatible with: + catch_analysis_tools.calibration.astrometry.run_astrometry_calibration + """ + if body is None: + body = kwargs + + cfg = validate_and_normalize_astrometry(body) + input_fits = materialize_input_fits(cfg["image_url"]) + + image = fitsio.read(input_fits).astype(np.float32) + bkg = sep.Background(image) + bkg_err = float(bkg.globalrms) + + ra_deg = cfg["ra"] if cfg["use_ra_dec"] else None + dec_deg = cfg["dec"] if cfg["use_ra_dec"] else None + + try: + astrom_res = run_astrometry_calibration( + input_fits=input_fits, + ra_deg=ra_deg, + dec_deg=dec_deg, + bkg_err=bkg_err, + pixel_scale=cfg["pixel_scale"], + snr=cfg["snr_threshold"], + output_fits=None, + ) + except Exception as exc: + raise AstrometrySolveError(f"Astrometry calibration failed: {exc}") from exc + + source_list = astrom_res["source_list"] + wcs_solution = astrom_res["wcs_solution"] + output_fits = astrom_res["output_fits"] + + ny, nx = image.shape + center_world = wcs_solution.pixel_to_world(nx / 2.0, ny / 2.0) + + return { + "wcs_image_url": output_fits, + "sources_detected": int(len(source_list)), + "center_ra_deg": float(center_world.ra.deg), + "center_dec_deg": float(center_world.dec.deg), + "pixel_scale": float(cfg["pixel_scale"]), + } \ No newline at end of file diff --git a/catch_analysis_tools/app/services/calibration/photometry.py b/catch_analysis_tools/app/services/calibration/photometry.py new file mode 100644 index 0000000..afcf86a --- /dev/null +++ b/catch_analysis_tools/app/services/calibration/photometry.py @@ -0,0 +1,426 @@ +import base64 +import io +import os +import tempfile +import urllib.parse +import urllib.request +from copy import deepcopy +from pathlib import Path +from typing import Any, Dict, Optional + +import fitsio +import matplotlib + +matplotlib.use("Agg") + +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 + + +class PhotometryValidationError(ValueError): + pass + + +class PhotometryCalibrationError(RuntimeError): + pass + + +DEFAULT_CONFIG = { + "photometry": { + "catalog": "PanSTARRS1", + "obs_band": "g", + "cal_band": "r", + }, + "detection": { + "snr": 3.0, + "aperture_radius": 7.0, + }, + "output": { + "make_plots": False, + "write_fits": True, + }, + "meta": { + "plot_type": "color_correction", + }, +} + + +def merge_config(user_config: Optional[dict], default_config: dict) -> dict: + result = deepcopy(default_config) + + if not user_config: + return result + + for key, value in user_config.items(): + if isinstance(value, dict) and isinstance(result.get(key), dict): + result[key] = merge_config(value, result[key]) + else: + result[key] = value + + return result + + +def _as_float(body: Dict[str, Any], name: str, default=None, required: bool = False): + value = body.get(name, default) + + if value is None: + if required: + raise PhotometryValidationError(f"{name} is required") + return None + + try: + return float(value) + except (TypeError, ValueError) as exc: + raise PhotometryValidationError(f"{name} must be a number") from exc + + +def _as_bool(body: Dict[str, Any], name: str, default: bool = False) -> bool: + value = body.get(name, default) + + if isinstance(value, bool): + return value + + if isinstance(value, str): + return value.strip().lower() in {"1", "true", "t", "yes", "y"} + + return bool(value) + + +def _safe_filename_from_url(file_url: str) -> str: + parsed = urllib.parse.urlparse(file_url) + name = os.path.basename(parsed.path) or "input.fits" + + if not name.lower().endswith((".fits", ".fit", ".fts")): + name = f"{name}.fits" + + safe = "".join( + ch if ch.isalnum() or ch in {".", "_", "-"} else "_" + for ch in name + ) + + return safe or "input.fits" + + +def _resolve_input_fits(file_url: str) -> str: + parsed = urllib.parse.urlparse(file_url) + + if parsed.scheme in {"http", "https"}: + out_dir = Path(tempfile.gettempdir()) / "catch_analysis_tools" + out_dir.mkdir(parents=True, exist_ok=True) + + out_path = out_dir / _safe_filename_from_url(file_url) + urllib.request.urlretrieve(file_url, out_path) + + return str(out_path) + + if parsed.scheme == "file": + local_path = urllib.request.url2pathname(parsed.path) + else: + local_path = file_url + + if not os.path.exists(local_path): + raise PhotometryValidationError( + f"WCS-solved FITS file not found: {file_url}" + ) + + return local_path + + +def validate_and_normalize(body: Dict[str, Any]) -> Dict[str, Any]: + if body is None: + raise PhotometryValidationError("Request body is required") + + wcs_image_url = body.get("wcs_image_url") + + if not wcs_image_url: + raise PhotometryValidationError("wcs_image_url is required") + + color_index = body.get("color_term") + obs_band = body.get("obs_band") + cal_band = body.get("cal_band", "r") + + if obs_band is None and isinstance(color_index, str) and "-" in color_index: + obs_band = color_index.split("-", 1)[0] + + if obs_band is None: + obs_band = "g" + + return { + "wcs_image_url": wcs_image_url, + "photometry": { + "catalog": body.get("catalog", "PanSTARRS1"), + "obs_band": obs_band, + "cal_band": cal_band, + }, + "detection": { + "snr": _as_float(body, "snr_threshold", 3.0), + "aperture_radius": _as_float(body, "aperture_radius", 7.0), + }, + "output": { + "make_plots": _as_bool(body, "return_plot", False), + "write_fits": _as_bool(body, "write_fits", True), + }, + "meta": { + "plot_type": body.get("plot_type", "color_correction"), + }, + } + + +def _detect_sources_on_background_subtracted_image( + image_sub: np.ndarray, + bkg_err: float, + snr: float, + aperture_radius: float, +) -> pd.DataFrame: + sep.set_sub_object_limit(500) + + sources = sep.extract( + image_sub, + thresh=snr, + err=bkg_err, + deblend_nthresh=16, + ) + + source_list = pd.DataFrame(sources) + + if source_list.empty: + raise PhotometryValidationError( + "No sources were detected in the WCS-solved image." + ) + + flux, flux_err, _ = sep.sum_circle( + 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 + ) + + if source_list.empty: + raise PhotometryValidationError( + "No detected sources have positive aperture_sum." + ) + + return source_list + + +def _load_wcs_from_fits(input_fits: str) -> WCS: + with fits.open(input_fits) as hdul: + wcs_solution = WCS(hdul[0].header, relax=True) + + if not wcs_solution.has_celestial: + raise PhotometryValidationError( + "Input FITS does not contain a celestial WCS." + ) + + return wcs_solution + + +def _attach_sky_coordinates( + source_list: pd.DataFrame, + wcs_solution: WCS, +) -> tuple[pd.DataFrame, SkyCoord]: + x = source_list["x"].to_numpy(dtype=float) + y = source_list["y"].to_numpy(dtype=float) + + world = wcs_solution.pixel_to_world(x, y) + + source_list = source_list.copy() + source_list["RA"] = world.ra.deg + source_list["Dec"] = world.dec.deg + + sky_coords = SkyCoord(source_list["RA"].to_numpy(), source_list["Dec"].to_numpy(), unit="deg") + + return source_list, sky_coords + + +def _load_photometric_assets(input_fits: str, det_cfg: dict) -> dict: + image = fitsio.read(input_fits).astype(np.float32) + + bkg = sep.Background(image) + image_sub = image - bkg.back() + + source_list = _detect_sources_on_background_subtracted_image( + image_sub=image_sub, + bkg_err=float(bkg.globalrms), + snr=det_cfg.get("snr", 3.0), + aperture_radius=det_cfg.get("aperture_radius", 7.0), + ) + + wcs_solution = _load_wcs_from_fits(input_fits) + source_list, sky_coords = _attach_sky_coordinates(source_list, wcs_solution) + + return { + "image": image, + "image_sub": image_sub, + "bkg": bkg, + "source_list": source_list, + "sky_coords": sky_coords, + "wcs_solution": wcs_solution, + } + + +def _encode_figure(fig) -> str: + buf = io.BytesIO() + fig.savefig(buf, format="png", dpi=150, bbox_inches="tight") + return base64.b64encode(buf.getvalue()).decode() + + +def _public_file_reference(path: str) -> str: + return path + + +def run_photometry_from_config(input_fits: str, user_config: dict) -> dict: + from catch_analysis_tools.calibration.photometry import ( + calibrate_photometric_zero_point, + get_color_corrected_indices, + get_matched_indices, + plot_color_correction, + plot_photometric_matches, + write_photometric_calibration_output, + ) + + config = merge_config(user_config, DEFAULT_CONFIG) + + phot_cfg = config["photometry"] + det_cfg = config["detection"] + out_cfg = config["output"] + meta_cfg = config["meta"] + + assets = _load_photometric_assets(input_fits, det_cfg) + + image = assets["image"] + image_sub = assets["image_sub"] + source_list = assets["source_list"] + sky_coords = assets["sky_coords"] + wcs_solution = assets["wcs_solution"] + + try: + calibration = calibrate_photometric_zero_point( + sky_coords=sky_coords, + source_list=source_list, + catalog=phot_cfg["catalog"], + obs_band=phot_cfg["obs_band"], + cal_band=phot_cfg["cal_band"], + ) + except Exception as exc: + raise PhotometryCalibrationError( + f"Photometric calibration failed: {exc}" + ) from exc + + zp = calibration["zp"] + color_term = calibration.get("color_term", calibration.get("C")) + unc = calibration["unc"] + m = calibration["m"] + m_inst = calibration["m_inst"] + color_mags = calibration["color_mags"] + color_index = calibration["color_index"] + objids = calibration["objids"] + + matched_idx = get_matched_indices(objids, len(source_list)) + colored_idx = get_color_corrected_indices(color_mags, len(source_list)) + + plots = {} + + if out_cfg.get("make_plots", False): + requested_plot = meta_cfg.get("plot_type", "color_correction") + + make_color = requested_plot in {"color_correction", "all"} + make_overlay = requested_plot in {"image_overlay", "all"} + + if make_color: + fig, _ = plot_color_correction( + color_mags, + m, + m_inst, + color_term, + zp, + color_index, + ) + plots["color_correction"] = _encode_figure(fig) + plt.close(fig) + + if make_overlay: + fig, _ = plot_photometric_matches( + image_sub, + source_list, + matched_idx, + colored_idx, + ) + plots["image_overlay"] = _encode_figure(fig) + plt.close(fig) + + photometry_image_url = None + + if out_cfg.get("write_fits", True): + file_base = os.path.splitext(input_fits)[0] + output_fits = f"{file_base}_photometry.fits" + + write_photometric_calibration_output( + image=image, + wcs_solution=wcs_solution, + source_list=source_list, + matched_idx=matched_idx, + color_corrected_idx=colored_idx, + output_fits=output_fits, + zero_point=zp, + zero_point_uncertainty=unc, + catalog=phot_cfg["catalog"], + obs_band=phot_cfg["obs_band"], + cal_band=phot_cfg["cal_band"], + color_index=color_index, + color_term=color_term, + ) + + photometry_image_url = _public_file_reference(output_fits) + + result = { + "photometry": { + "zero_point": float(zp), + "color_term": float(color_term), + "uncertainty": float(unc), + }, + "sources": { + "detected": int(len(source_list)), + "matched": int(len(matched_idx)), + }, + } + + if photometry_image_url is not None: + result["photometry_image_url"] = photometry_image_url + + if plots: + result["plots"] = plots + + return result + + +def run_photometry_pipeline(body=None, **kwargs): + if body is None: + body = kwargs + + config = validate_and_normalize(body) + input_fits = _resolve_input_fits(config["wcs_image_url"]) + + result = run_photometry_from_config(input_fits, config) + + if config["output"]["make_plots"]: + plot_type = config["meta"]["plot_type"] + png_b64 = result["plots"][plot_type] + png_bytes = base64.b64decode(png_b64) + + return png_bytes, 200, {"Content-Type": "image/png"} + + return result, 200, {"Content-Type": "application/json"} \ No newline at end of file diff --git a/catch_analysis_tools/astrometry.py b/catch_analysis_tools/astrometry.py deleted file mode 100644 index 706821e..0000000 --- a/catch_analysis_tools/astrometry.py +++ /dev/null @@ -1,474 +0,0 @@ -import os -import subprocess -import argparse -import numpy as np -import pandas as pd -import sep -import fitsio -import matplotlib.pyplot as plt -from astropy.io import fits -from astropy.wcs import WCS -from astropy.table import Table -from astropy.coordinates import SkyCoord -import calviacat as cvc - - -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. - - 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.") - 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." - ) - - 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", - input_fits, - ] - - try: - subprocess.run(command, check=True) - return True - except subprocess.CalledProcessError as e: - raise RuntimeError(f"solve-field failed: {e}") - - -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) - flux, flux_err, _ = sep.sum_circle( - 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) - return source_list, image_sub - - -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: - wcs_solution = WCS(hdul[0].header) - return wcs_solution - - -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] - 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'): - """ - 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'). - cal_band : str, optional - Reference catalog filter for color term (e.g. 'g', 'r', 'i'; default 'g'). - - Returns - ------- - calibration : dict - Dictionary with keys: - - zp : float, zero-point magnitude - - C : float, color-term coefficient - - unc : float, uncertainty of zero-point - - m : array_like, calibrated magnitudes in the observed band - - m_inst : array_like, instrumental magnitudes - - obs_band : str, label of the observed band - - cal_band : str, same as input cal_band - - color_mags : array_like, color indices (obs_band - cal_band) - - 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: - CatalogClass = getattr(cvc, catalog) - except AttributeError: - raise ValueError(f"Catalog '{catalog}' not found in calviacat") - - 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) - zp, C, unc, m_cal, color_mags, _ = ref.cal_color( - objids, - m_inst, - cal_band, - 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, - } - - -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='.') - 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)') - plt.tight_layout() - return fig, ax - - -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) - 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.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): - """ - 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()) - 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') - else: - 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) - - -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', '.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) - - -if __name__ == "__main__": - - parser = argparse.ArgumentParser( - 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("--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 - file_base = os.path.splitext(input_fits)[0] - image = fitsio.read(input_fits).astype(np.float32) - Ra = args.Ra - Dec = args.Dec - bkg_err = args.bkg_err - pixel_scale = args.pixel_scale - snr = args.snr - catalog = args.catalog - obs_band = args.obs_band - cal_band = args.cal_band - - output_wcs = f"{file_base}.wcs" - try: - if run_solve_field(input_fits, output_wcs, pixel_scale, Ra, Dec): - wcs_solution = load_wcs(output_wcs) - else: - raise RuntimeError("solve-field did not produce a WCS solution") - except Exception as e: - raise SystemExit(f"WCS calibration failed: {e}") - - source_list, telescope_image_sub = find_sources(image, bkg_err, snr) - source_list, sky_coords = retrieve_sources(source_list, wcs_solution) - calibration = calibrate_photometry( - sky_coords, - source_list, - catalog=catalog, - obs_band=obs_band, - cal_band=cal_band, - ) - zp = calibration["zp"] - C = calibration["C"] - unc = calibration["unc"] - m = calibration["m"] - m_inst = calibration["m_inst"] - color_mags = calibration['color_mags'] - color_index = calibration['color_index'] - objids = calibration['objids'] - - if hasattr(objids, "mask"): - matched_idx = np.where(~objids.mask)[0] - else: - matched_idx = np.arange(len(source_list)) - if hasattr(color_mags, "mask"): - colored_idx = np.where(~color_mags.mask)[0] - 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) - plt.show() - - create_header( - image, - wcs_solution, - zp, - unc, - source_list, - matched_idx, - colored_idx, - input_fits, - obs_band, - catalog, - cal_band, - ) - - cleanup_files(file_base) diff --git a/catch_analysis_tools/calibration/astrometry.py b/catch_analysis_tools/calibration/astrometry.py index d8f8dca..dc231eb 100644 --- a/catch_analysis_tools/calibration/astrometry.py +++ b/catch_analysis_tools/calibration/astrometry.py @@ -47,16 +47,21 @@ def run_solve_field(input_fits, output_wcs, pixel_scale, Ra_deg, Dec_deg, scale_ "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", - input_fits, ] + if Ra_deg is not None and Dec_deg is not None: + command.extend([ + "--ra", str(Ra_deg), + "--dec", str(Dec_deg), + "--radius", "2", + ]) + + command.append(input_fits) + try: subprocess.run(command, check=True) return True @@ -146,9 +151,12 @@ 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] + world = wcs_solution.pixel_to_world( + np.asarray(source_list["x"], dtype=float), + np.asarray(source_list["y"], dtype=float), + ) + source_list["RA"] = world.ra.deg + source_list["Dec"] = world.dec.deg sky_coords = SkyCoord(source_list['RA'], source_list['Dec'], unit='deg') return source_list, sky_coords @@ -232,7 +240,9 @@ def run_astrometry_calibration( else: raise RuntimeError("solve-field did not produce a WCS solution.") - source_list, image_sub = find_sources(image, bkg_err, snr) + bkg = sep.Background(image) + image_sub = image - bkg.back() + source_list, image_sub = find_sources(image_sub, bkg.globalrms, snr) source_list, sky_coords = retrieve_sources(source_list, wcs_solution) write_astrometry_output( @@ -321,3 +331,77 @@ def run_astrometry_calibration( except Exception as exc: raise SystemExit(f"Astrometric calibration failed: {exc}") +import os +from typing import Dict, Any +import numpy as np +import fitsio +import sep +from astropy.io import fits +from catch_analysis_tools.calibration.astrometry import run_astrometry_calibration + +class AstrometryValidationError(Exception): + pass + +def validate_and_normalize_astrometry(body: Dict[str, Any]) -> Dict[str, Any]: + def get_float(key: str, default: Any = None) -> Any: + val = body.get(key, default) + return float(val) if val is not None else None + + def get_bool(key: str, default: bool = False) -> bool: + val = body.get(key, default) + if isinstance(val, str): + return val.lower() in ("true", "1", "yes") + return bool(val) + + image_url = body.get("image_url") + if not image_url: + raise AstrometryValidationError("image_url is required") + + return { + "image_url": image_url, + "ra": get_float("ra"), + "dec": get_float("dec"), + "pixel_scale": get_float("pixel_scale", 2.5), + "scale_low": get_float("scale_low"), + "scale_high": get_float("scale_high"), + "search_radius": get_float("search_radius", 2.0), + "use_ra_dec": get_bool("use_ra_dec", False), + "snr_threshold": get_float("snr_threshold", 3.0), + } + +def run_pipeline(body: dict) -> dict: + cfg = validate_and_normalize_astrometry(body) + input_fits = cfg["image_url"] + + image = fitsio.read(input_fits).astype(np.float32) + bkg = sep.Background(image) + + astrom_res = run_astrometry_calibration( + input_fits=input_fits, + ra_deg=cfg["ra"] if cfg["use_ra_dec"] else None, + dec_deg=cfg["dec"] if cfg["use_ra_dec"] else None, + bkg_err=bkg.globalrms, + pixel_scale=cfg["pixel_scale"], + snr=cfg["snr_threshold"], + output_fits=None + ) + + source_list = astrom_res["source_list"] + wcs_solution = astrom_res["wcs_solution"] + + ny, nx = image.shape + center_world = wcs_solution.pixel_to_world(nx / 2.0, ny / 2.0) + + output_filename = input_fits.replace(".fits", ".wcs.fits") + with fits.open(input_fits) as hdul: + header = hdul[0].header + header.update(wcs_solution.to_header()) + fits.writeto(output_filename, hdul[0].data, header, overwrite=True) + + return { + "wcs_image_url": output_filename, + "sources_detected": int(len(source_list)), + "center_ra_deg": float(center_world.ra.deg), + "center_dec_deg": float(center_world.dec.deg), + "pixel_scale": cfg["pixel_scale"], + } \ No newline at end of file diff --git a/catch_analysis_tools/calibration/photometry.py b/catch_analysis_tools/calibration/photometry.py index 6ccb709..f1d8695 100644 --- a/catch_analysis_tools/calibration/photometry.py +++ b/catch_analysis_tools/calibration/photometry.py @@ -11,7 +11,7 @@ def calibrate_photometric_zero_point( source_list: pd.DataFrame, catalog: str = "PanSTARRS1", obs_band: str = "obs_band", - cal_band: str = "g", + cal_band: str = "r", catalog_db: str = "cat.db", ): """ @@ -73,10 +73,33 @@ def calibrate_photometric_zero_point( if len(results[0]) < 500: ref.fetch_field(sky_coords) - objids, distances = ref.xmatch(sky_coords) + xmatch_result = ref.xmatch(sky_coords) + + if xmatch_result is None: + raise RuntimeError("Photometric calibration failed: fewer than 10 catalog matches.") + + objids, distances = xmatch_result + + aperture_sum = np.asarray(source_list["aperture_sum"].values, dtype=float) + valid_flux = np.isfinite(aperture_sum) & (aperture_sum > 0) + + objids = objids[valid_flux] + distances = distances[valid_flux] + aperture_sum = aperture_sum[valid_flux] m_inst = -2.5 * np.log10(aperture_sum) + valid_m_inst = np.isfinite(m_inst) + + objids = objids[valid_m_inst] + distances = distances[valid_m_inst] + m_inst = m_inst[valid_m_inst] + + if len(m_inst) < 10: + raise RuntimeError( + f"Photometric calibration failed: only {len(m_inst)} valid matched sources remain after filtering." + ) + zp, color_term, zp_unc, m_cal, color_mags, _ = ref.cal_color( objids, m_inst,