From 835f6f3970c11e26e74df9dbc6bac8e4bfc2d3f9 Mon Sep 17 00:00:00 2001 From: Ola Wenda Date: Sat, 23 Aug 2025 15:48:21 +0200 Subject: [PATCH] Fix SDSS catalog implementation --- catalog.py | 170 +++++++++++++++++++++++++++++++++++++-------- config.py | 3 +- dophot3.py | 64 ++++++++--------- filter_matching.py | 3 +- 4 files changed, 177 insertions(+), 63 deletions(-) diff --git a/catalog.py b/catalog.py index bfc477b..ec8191d 100644 --- a/catalog.py +++ b/catalog.py @@ -12,7 +12,13 @@ from astropy.coordinates import SkyCoord import astropy.units as u import logging +import subprocess +import tempfile +import time +from sklearn.neighbors import KDTree +from dataclasses import asdict, dataclass +from typing import Any, Dict, Optional, Tuple, Type, TypeVar, cast # Type aliases TableType = TypeVar("TableType", bound=astropy.table.Table) @@ -86,14 +92,13 @@ class CatalogFilters: 'R2': CatalogFilter('R2', 6400, 'AB', 'e_R2mag'), 'I': CatalogFilter('I', 8100, 'AB', 'e_Imag'), } - - # SDSS filters (using Sloan naming convention) + # SDSS filter SDSS = { - 'Sloan_u': CatalogFilter('Sloan_u', 3551, 'AB', 'Sloan_u_err'), - 'Sloan_g': CatalogFilter('Sloan_g', 4686, 'AB', 'Sloan_g_err'), - 'Sloan_r': CatalogFilter('Sloan_r', 6166, 'AB', 'Sloan_r_err'), - 'Sloan_i': CatalogFilter('Sloan_i', 7480, 'AB', 'Sloan_i_err'), - 'Sloan_z': CatalogFilter('Sloan_z', 8932, 'AB', 'Sloan_z_err'), + "upmag": CatalogFilter("upmag", 3551, "AB", "e_upmag"), + "gpmag": CatalogFilter("gpmag", 4686, "AB", "e_gpmag"), + "rpmag": CatalogFilter("rpmag", 6166, "AB", "e_rpmag"), + "ipmag": CatalogFilter("ipmag", 7480, "AB", "e_ipmag"), + "zpmag": CatalogFilter("zpmag", 8932, "AB", "e_zpmag"), } class Catalog(astropy.table.Table): @@ -109,7 +114,7 @@ class Catalog(astropy.table.Table): GAIA: str = 'gaia' MAKAK: str = 'makak' USNOB: str = 'usno' - SDSS: str = 'sdss' + SDSS: str = "SDSS" # Define available catalogs with their properties KNOWN_CATALOGS: Dict[str, CatalogConfig] @@ -224,28 +229,29 @@ class Catalog(astropy.table.Table): 'pmDE': 'pmdec' } }, + SDSS: { - 'description': 'SDSS Data Release 16', - 'filters': CatalogFilters.SDSS, - 'epoch': 2000.0, - 'local': False, - 'service': 'VizieR', - 'catalog_id': 'V/154/sdss16', - 'column_mapping': { - 'RA_ICRS': 'radeg', - 'DE_ICRS': 'decdeg', - 'upmag': 'Sloan_u', - 'gpmag': 'Sloan_g', - 'rpmag': 'Sloan_r', - 'ipmag': 'Sloan_i', - 'zpmag': 'Sloan_z', - 'e_upmag': 'Sloan_u_err', - 'e_gpmag': 'Sloan_g_err', - 'e_rpmag': 'Sloan_r_err', - 'e_ipmag': 'Sloan_i_err', - 'e_zpmag': 'Sloan_z_err', - 'pmRA': 'pmra', - 'pmDE': 'pmdec', + "description": "SDSS Data Release 16", + "filters": CatalogFilters.SDSS, + "epoch": 2000.0, + "local": False, + "service": "VizieR", + "catalog_id": "V/154/sdss16", + "column_mapping": { + "RA_ICRS": "radeg", + "DE_ICRS": "decdeg", + "upmag": "upmag", + "gpmag": "gpmag", + "rpmag": "rpmag", + "ipmag": "ipmag", + "zpmag": "zpmag", + "e_upmag": "e_upmag", + "e_gpmag": "e_gpmag", + "e_rpmag": "e_rpmag", + "e_ipmag": "e_ipmag", + "e_zpmag": "e_zpmag", + "pmRA": "pmra", + "pmDE": "pmdec", } } } @@ -323,7 +329,7 @@ def _fetch_catalog_data(self) -> Optional[astropy.table.Table]: elif self._catalog_name == self.USNOB: result = self._get_usnob_data() elif self._catalog_name == self.SDSS: - result = self._get_sdss_data() + result = self._get_sdss_online() else: raise ValueError(f"Unknown catalog: {self._catalog_name}") @@ -594,6 +600,110 @@ def _get_gaia_data(self) -> Optional[astropy.table.Table]: except Exception as e: raise ValueError(f"Gaia query failed: {str(e)}") + + def _get_atlas_vizier(self) -> Optional[astropy.table.Table]: + """Get ATLAS RefCat2 data from VizieR with updated column mapping.""" + from astroquery.vizier import Vizier + + # Configure Vizier with correct column names + column_mapping = self.KNOWN_CATALOGS[self.ATLAS_VIZIER]["column_mapping"] + vizier = Vizier( + columns=list(column_mapping.keys()), #Określa, że VizieR powinien zwrócić tylko te kolumny, których nazwy są kluczami w słowniku column_mapping + column_filters={ + "rmag": f"<{self._query_params.mlim}" # Magnitude limit in r-band + }, + row_limit=-1, #Określa, że VizieR ma zwrócić wszystkie pasujące wiersze (obiekty), a nie tylko ograniczoną liczbę. + ) + + catalog=self.KNOWN_CATALOGS[self.ATLAS_VIZIER]["catalog_id"] + + cat = self._get_vizier_data(catalog, column_mapping, vizier) + + # Add computed Johnson magnitudes + self._add_transformed_magnitudes(cat) + + return cat + + # except Exception as e: + # warnings.warn(f"VizieR ATLAS query failed: {e}") + # return None + + + def _get_sdss_online(self): + """Get SDSS DR16 data from VizieR.""" + try: + from astroquery.vizier import Vizier + + config = self.KNOWN_CATALOGS[self.SDSS] + column_mapping = config["column_mapping"] +# catalog=config["catalog_id"] + + # Configure Vizier + vizier = Vizier( + columns=list(column_mapping.keys()), + column_filters={ + "rpmag": f"<{self._query_params.mlim}" # Magnitude limit in R1 + }, + row_limit=-1, # Get all matching objects + ) + + # Create coordinate object + coords = SkyCoord( + ra=self._query_params.ra *u.deg, + dec=self._query_params.dec *u.deg, + frame='icrs' + ) + + + # Query VizieR + result = vizier.query_region( + coords, + width=2*self._query_params.width * u.deg, + height=2*self._query_params.height * u.deg, + catalog=config['catalog_id'] + ) + + if not result or len(result) == 0: + print("No SDSS data found") + return None + + sdss = result[0] + + # Create output catalog + cat = astropy.table.Table() + + # Initialize mapped columns + our_columns = set(column_mapping.values()) + for col in our_columns: + cat[col] = np.zeros(len(sdss), dtype=np.float64) + + # Map columns according to configuration + for vizier_name, our_name in column_mapping.items(): + if vizier_name in sdss.columns: + # Convert proper motions from mas/yr to deg/yr if needed + if vizier_name in ['pmRA', 'pmDE']: + cat[our_name] = sdss[vizier_name] / (3.6e6) + else: + cat[our_name] = sdss[vizier_name] + + + # Handle quality flags and uncertainties + for band in ["u", "g", "r", "i", "z"]: + mag_col = f"{band}pmag" + err_col = f"e_{band}pmag" + if mag_col in cat.columns: + # Set typical errors if not provided + if err_col not in cat.columns or np.all(cat[err_col] == 0): + cat[err_col] = np.where( + cat[mag_col] < 19, + 0.1, # Brighter stars + 0.2, # Fainter stars + ) + + return cat + + except Exception as e: + raise ValueError(f"SDSS query failed: {str(e)}") def _get_usnob_data(self) -> Optional[astropy.table.Table]: """Get USNO-B1.0 data from VizieR""" diff --git a/config.py b/config.py index c89ef37..9bd26bd 100755 --- a/config.py +++ b/config.py @@ -28,9 +28,10 @@ SloanU = Sloan_u, Sloan_g, Sloan_r, Sloan_i, Sloan_z SloanJ = Sloan_g, Sloan_r, Sloan_i, Sloan_z, J Gaia = BP, G, RP, BP, RP -GaiaJ = Johnson_B, Johnson_V, Johnson_R, Johnson_I, Johnson_I # Duplicate last filter to create zero color4, avoiding mixed AB/Vega systems +GaiaJ = Johnson_B, Johnson_V, Johnson_R, Johnson_I, Johnson_I PS = g, r, i, z, y USNO = B2, R2, I, R1, B1 +SDSS = upmag, gpmag, rpmag, ipmag, zpmag """ DEFAULT_CONFIG_FILE = '~/.config/dophot3/config' diff --git a/dophot3.py b/dophot3.py index 00f0340..442ac1e 100755 --- a/dophot3.py +++ b/dophot3.py @@ -604,47 +604,49 @@ def main(): imgno=0 for arg in options.files: + try: + det = process_input_file(arg, options) + if det is None: + logging.warning(f"I do not know what to do with {arg}") + continue - det = process_input_file(arg, options) - if det is None: - logging.warning(f"I do not know what to do with {arg}") - continue - - catalog_name = 'makak' if options.makak else (options.catalog or 'atlas@localhost') - determine_filter(det, options, catalog_name) - logging.info(f'Reference filter is {det.meta["PHFILTER"]}, ' - f'Schema: {det.meta["PHSCHEMA"]}, ' - f'System: {det.meta["PHSYSTEM"]}') - - start = time.time() + catalog_name = 'makak' if options.makak else (options.catalog or 'atlas@localhost') + determine_filter(det, options, catalog_name) + logging.info(f'Reference filter is {det.meta["PHFILTER"]}, ' + f'Schema: {det.meta["PHSCHEMA"]}, ' + f'System: {det.meta["PHSYSTEM"]}') + + start = time.time() - cat, matches, imgwcs, target_match = process_image_with_dynamic_limits(det, options) - if cat is None: - logging.warning(f"Failed to process {arg}, skipping") - continue + cat, matches, imgwcs, target_match = process_image_with_dynamic_limits(det, options) + if cat is None: + logging.warning(f"Failed to process {arg}, skipping") + continue - logging.info(f"Catalog processing took {time.time()-start:.3f}s") + logging.info(f"Catalog processing took {time.time()-start:.3f}s") - det['ALPHA_J2000'], det['DELTA_J2000'] = imgwcs.all_pix2world( [det['X_IMAGE']], [det['Y_IMAGE']], 1) + det['ALPHA_J2000'], det['DELTA_J2000'] = imgwcs.all_pix2world( [det['X_IMAGE']], [det['Y_IMAGE']], 1) - logging.info("Reference filter is %s, Photometric schema: %s, Photometric system: %s"\ - %(det.meta['PHFILTER'],det.meta['PHSCHEMA'],det.meta['PHSYSTEM'])) + logging.info("Reference filter is %s, Photometric schema: %s, Photometric system: %s"\ + %(det.meta['PHFILTER'],det.meta['PHSCHEMA'],det.meta['PHSYSTEM'])) - # make pairs to be fitted - det.meta['IMGNO'] = imgno - n_matched_stars = make_pairs_to_fit(det, cat, matches, imgwcs, options, data, None) - det.meta['IDNUM'] = n_matched_stars + # make pairs to be fitted + det.meta['IMGNO'] = imgno + n_matched_stars = make_pairs_to_fit(det, cat, matches, imgwcs, options, data, None) + det.meta['IDNUM'] = n_matched_stars - if n_matched_stars == 0: - logging.warning(f"No matched stars in {det.meta['FITSFILE']}, skipping image") - continue + if n_matched_stars == 0: + logging.warning(f"No matched stars in {det.meta['FITSFILE']}, skipping image") + continue - metadata.append(det.meta) - alldet.append(det) - target.append(target_match) + metadata.append(det.meta) + alldet.append(det) + target.append(target_match) - imgno += 1 + imgno += 1 + except Exception as e: + logging.error(f"Error processing {arg}: {e}") data.finalize() diff --git a/filter_matching.py b/filter_matching.py index 475e396..bce621a 100644 --- a/filter_matching.py +++ b/filter_matching.py @@ -106,7 +106,8 @@ def get_base_filter(det, options, catalog_name): 'G': 5890, 'BP': 5050, 'RP': 7730, 'Sloan_g': 4810, 'Sloan_r': 6170, 'Sloan_i': 7520, 'Sloan_z': 8660, 'Johnson_U': 3600, 'Johnson_B': 4353, 'Johnson_V': 5477, - 'Johnson_R': 6349, 'Johnson_I': 8797, 'N': 6000 + 'Johnson_R': 6349, 'Johnson_I': 8797, 'N': 6000, + 'upmag': 3551, 'gpmag': 4686, 'rpmag': 6166, 'ipmag': 7480, 'zpmag': 8932, } target_wavelength = FILTER_WAVELENGTHS.get(filter_name)