diff --git a/python/lsst/pipe/tasks/finalizeCharacterization.py b/python/lsst/pipe/tasks/finalizeCharacterization.py index 1b6f2a415..21be6e65e 100644 --- a/python/lsst/pipe/tasks/finalizeCharacterization.py +++ b/python/lsst/pipe/tasks/finalizeCharacterization.py @@ -38,6 +38,7 @@ import astropy.units as u import numpy as np import esutil +import galsim from smatch.matcher import Matcher @@ -50,7 +51,10 @@ from lsst.meas.algorithms import MeasureApCorrTask from lsst.meas.base import SingleFrameMeasurementTask, ApplyApCorrTask from lsst.meas.algorithms.sourceSelector import sourceSelectorRegistry +from lsst.afw.cameraGeom import PIXELS, FOCAL_PLANE +from lsst.geom import Point2D +from .lsstCamFocalPlaneHeight import make_metrology_table, get_height_interpolator from .reserveIsolatedStars import ReserveIsolatedStarsTask from lsst.obs.base.utils import TableVStack @@ -439,6 +443,11 @@ def _make_output_schema_mapper(self, input_schema): doc="Maximum value in the star image used to train PSF.", doReplace=True, ) + output_schema.addField( + 'sensor_height', + type=np.float32, + doc="sensor height in PSF fit." + ) alias_map = input_schema.getAliasMap() alias_map_output = afwTable.AliasMap() @@ -489,7 +498,11 @@ def _make_selection_schema_mapper(self, input_schema): doc="Maximum value in the star image used to train PSF.", doReplace=True, ) - + selection_schema.addField( + 'sensor_height', + type=np.float32, + doc="Color used in PSF fit." + ) return mapper, selection_schema def concat_isolated_star_cats(self, band, isolated_star_cat_dict, isolated_star_source_dict): @@ -621,6 +634,26 @@ def add_src_colors(self, srcCat, fgcmCat, band): srcCat[idSrc]['psf_color_value'] = colors[idColor] srcCat[idSrc]['psf_color_type'] = f"{magStr1}-{magStr2}" + def add_src_itl_height(self, srcCat, exposure): + + #'slot_Centroid_x', 'slot_Centroid_y', + #'sensor_height' + + metrology_table = make_metrology_table() + height_interpolator = get_height_interpolator(metrology_table) + + for i in range(len(srcCat)): + x = srcCat[i]['slot_Centroid_x'] + y = srcCat[i]['slot_Centroid_y'] + + detector = exposure.getDetector() + + image_pos = Point2D(x, y) #galsim.PositionD(x, y) # In CCD pixel coordinates + focal_pos = detector.transform(image_pos, PIXELS, FOCAL_PLANE) + dz = height_interpolator(focal_pos.getX(), focal_pos.getY()) + srcCat[i]['psf_color_value'] = dz + srcCat[i]['psf_color_type'] = 'dz' + def compute_psf_and_ap_corr_map(self, visit, detector, exposure, src, isolated_source_table, fgcm_standard_star_cat): """Compute psf model and aperture correction map for a single exposure. @@ -709,6 +742,9 @@ def compute_psf_and_ap_corr_map(self, visit, detector, exposure, src, self.add_src_colors(selected_src, fgcm_standard_star_cat, band) self.add_src_colors(measured_src, fgcm_standard_star_cat, band) + self.add_src_itl_height(selected_src, exposure) + self.add_src_itl_height(measured_src, exposure) + # Select the psf candidates from the selection catalog try: psf_selection_result = self.make_psf_candidates.run(selected_src, exposure=exposure) diff --git a/python/lsst/pipe/tasks/lsstCamFocalPlaneHeight.py b/python/lsst/pipe/tasks/lsstCamFocalPlaneHeight.py new file mode 100644 index 000000000..5b336dad5 --- /dev/null +++ b/python/lsst/pipe/tasks/lsstCamFocalPlaneHeight.py @@ -0,0 +1,143 @@ +import numpy as np +from astropy.table import Table +import re +from sklearn.neighbors import KNeighborsRegressor +import tqdm +from lsst.daf.butler import Butler +import lsst.afw.cameraGeom as cameraGeom +import lsst.geom as geom +from lsst.obs.lsst import LsstCam +import lsst.afw.image as afwImage +from lsst.geom import Point2D +from lsst.afw.cameraGeom import PIXELS, FOCAL_PLANE + +from astropy.io import fits +from astropy.table import Table + + +def make_metrology_table(file="/sdf/home/l/laliotis/rubin-user/lsst_dev/meas_extensions_piff/python/lsst/meas/extensions/piff/LSST_FP_cold_b_measurement_4col_bysurface.fits", rsid=None, write=False): + """ + Make an astropy table of the height measurement data. + Inputs: + file: string, file path for measurement file + rsid: string (optional) like R##_S## if you want data for just one sensor + write: bool (default False), whether to write out the table as a fits file + Outputs: + bigtable: One large astropy table with focal plane x and y coordinates, modeled and measured z values, and the RSID for which detector each fpx,fpy coord pair is on + """ + + rows = [] + with fits.open(file) as hdulist: + for hdu in tqdm.tqdm(hdulist): + if isinstance(hdu, fits.BinTableHDU): + table = Table(hdu.data) + extname = hdu.header['EXTNAME'] + if rsid is not None: + if extname == rsid: # filter to the single det , 172 + extname = re.sub(r'(R\d\d)(S\d\d)', r'\1_\2', extname) + for x, y, z_mod, z_meas in zip(table['X_CCS'], table['Y_CCS'], table['Z_CCS_MODEL'], table['Z_CCS_MEASURED']): + fpx = y + fpy = x + rows.append([fpx, fpy, z_mod, z_meas, extname]) + else: + if re.fullmatch(r'R\d\dS\d\d', extname): + extname = re.sub(r'(R\d\d)(S\d\d)', r'\1_\2', extname) + for x, y, z_mod, z_meas in zip(table['X_CCS'], table['Y_CCS'], table['Z_CCS_MODEL'], table['Z_CCS_MEASURED']): + fpx = y + fpy = x + rows.append([fpx, fpy, z_mod, z_meas, extname]) + + bigtable = Table(rows=rows, names=['fpx', 'fpy', 'z_mod', 'z_meas', 'det']) + if write: bigtable.write('metrology_fp.fits', format='fits', overwrite=True) + + return bigtable + + + +def get_height_interpolator(metrology_table, k=3, weight_type='distance', ztype='measured'): + """ + Create a KNN interpolator for height given the metrology table created above + Inputs: + table: astropy.table.Table, Table with columns 'fpx', 'fpy', and 'z_meas'. + k : int, Number of neighbors to use. + weight_type: str, 'uniform' or 'distance' weighting for KNN. + Returns: + interp_func : function, function that takes (fpx, fpy) and returns interpolated z_meas. + """ + + x=np.column_stack((metrology_table['fpx'], metrology_table['fpy'])) + if ztype=='measured': y=np.array(metrology_table['z_meas']) + elif ztype=='model': y=np.array(metrology_table['z_mod']) + else: raise ValueError('Specify a valid z type (measured, model)') + + knn = KNeighborsRegressor(n_neighbors=k, weights=weight_type) + knn.fit(x, y) + + def interp_func(x, y): + x = np.atleast_1d(x) + y = np.atleast_1d(y) + point = np.column_stack((x, y)) + z = knn.predict(point) + return z[0] #because z is an array + + + return interp_func + + + +def detector_fromcoords(x, y, camera): + """ + Function to get the detector ID for a given coordinate pair + Inputs: + x: float, x-coordinate + y: float, y-coordinate + camera: lsst.obs.lsst object, the camera being used (lsstcam or comcam) + """ + + cam = camera.getCamera() + point = geom.Point2D(x, y) + + for det in cam: + bbox = det.getBBox() + transform = det.getTransform(FOCAL_PLANE, det.PIXELS) + try: + pixel_point = transform.applyForward(point) + if bbox.contains(pixel_point): + return det.getId() + except Exception: + continue + + raise ValueError(f"No detector contains point ({x},{y}) in focal plane coords.") + + +def pixel_to_focal(x_arr, y_arr, det_arr): + """ + Convert pixel coordinates to focal plane coordinates for a given detector. + + Inputs: + x: float, x-coordinate in pixel space + y: float, y-coordinate in pixel space + det: int, detector ID + + Returns: + focal_point: lsst.geom.Point2D, focal plane coordinates + """ + + cam = LsstCam.getCamera() + + # Cache transforms to avoid redundant lookups + transform_cache = {} + fpx_list = [] + fpy_list = [] + + for x, y, d in zip(x_arr, y_arr, det_arr): + if d not in transform_cache: + det = cam[d] + transform_cache[d] = det.getTransform(PIXELS, FOCAL_PLANE) + + tx = transform_cache[d] + pt = tx.applyForward(Point2D(x, y)) + fpx_list.append(pt.getX()) + fpy_list.append(pt.getY()) + + return np.array(fpx_list), np.array(fpy_list)