Fix SDSS catalog implementation

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"""

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'

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()
 

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)