Devbugs

astrophot/image/image_object.py

@@ -448,13 +448,16 @@ def save(self, filename: str):
         hdulist = fits.HDUList(self.fits_images())
         hdulist.writeto(filename, overwrite=True)
 
-    def load(self, filename: str, hduext: int = 0):
+    def load(self, filename: Union[str, fits.HDUList], hduext: int = 0):
         """Load an image from a FITS file. This will load the primary HDU
         and set the data, CD, crpix, crval, and crtan attributes
         accordingly. If the WCS is not tangent plane, it will warn the user.
 
         """
-        hdulist = fits.open(filename)
+        if isinstance(filename, str):
+            hdulist = fits.open(filename)
+        else:
+            hdulist = filename
         self.data = np.array(hdulist[hduext].data, dtype=np.float64)
 
         self.CD = (

astrophot/models/func/gaussian_ellipsoid.py

@@ -17,7 +17,7 @@ def euler_rotation_matrix(alpha: ArrayLike, beta: ArrayLike, gamma: ArrayLike) -
         (
             backend.stack((ca * cg - cb * sa * sg, -ca * sg - cb * cg * sa, sb * sa)),
             backend.stack((cg * sa + ca * cb * sg, ca * cb * cg - sa * sg, -ca * sb)),
-            backend.stack((sb * cg, sb * cg, cb)),
+            backend.stack((sb * sg, sb * cg, cb)),
         ),
         dim=-1,
     )

astrophot/models/gaussian_ellipsoid.py

@@ -130,6 +130,6 @@ def brightness(
         v = backend.stack(self.transform_coordinates(x, y), dim=0).reshape(2, -1)
         return (
             flux
-            * backend.sum(backend.exp(-0.5 * (v * (inv_Sigma @ v))), dim=0)
+            * backend.exp(-0.5 * backend.sum(v * (inv_Sigma @ v), dim=0))
             / (2 * np.pi * backend.sqrt(backend.linalg.det(Sigma2D)))
         ).reshape(x.shape)

astrophot/utils/__init__.py

@@ -6,6 +6,7 @@
     interpolate,
     parametric_profiles,
 )
+from .fitsopen import ls_open
 
 __all__ = [
     "decorators",
@@ -14,4 +15,5 @@
     "parametric_profiles",
     "initialize",
     "conversions",
+    "ls_open",
 ]

astrophot/utils/fitsopen.py

@@ -0,0 +1,117 @@
+import numpy as np
+import warnings
+from astropy.utils.data import download_file
+from astropy.io import fits
+from astropy.utils.exceptions import AstropyWarning
+from numpy.core.defchararray import startswith
+
+try:
+    from pyvo.dal import sia
+except:
+    sia = None
+import os
+
+# Suppress common Astropy warnings that can clutter CI logs
+warnings.simplefilter("ignore", category=AstropyWarning)
+
+
+def flip_hdu(hdu):
+    """
+    Flips the image data in the FITS HDU on the RA axis to match the expected orientation.
+
+    Args:
+        hdu (astropy.io.fits.HDUList): The FITS HDU to be flipped.
+
+    Returns:
+        astropy.io.fits.HDUList: The flipped FITS HDU.
+    """
+    assert "CD1_1" in hdu[0].header, "HDU does not contain WCS information."
+    assert "CD2_1" in hdu[0].header, "HDU does not contain WCS information."
+    assert "CRPIX1" in hdu[0].header, "HDU does not contain WCS information."
+    assert "NAXIS1" in hdu[0].header, "HDU does not contain WCS information."
+    hdu[0].data = hdu[0].data[:, ::-1].copy()
+    hdu[0].header["CD1_1"] = -hdu[0].header["CD1_1"]
+    hdu[0].header["CD2_1"] = -hdu[0].header["CD2_1"]
+    hdu[0].header["CRPIX1"] = int(hdu[0].header["NAXIS1"] / 2) + 1
+    hdu[0].header["CRPIX2"] = int(hdu[0].header["NAXIS2"] / 2) + 1
+    return hdu
+
+
+def ls_open(ra, dec, size_arcsec, band="r", release="ls_dr9"):
+    """
+    Retrieves and opens a FITS cutout from the deepest stacked image in the
+    specified Legacy Survey data release using the Astro Data Lab SIA service.
+
+    Args:
+        ra (float): Right Ascension in decimal degrees.
+        dec (float): Declination in decimal degrees.
+        size_arcsec (float): Size of the square cutout (side length) in arcseconds.
+        band (str): The filter band (e.g., 'g', 'r', 'z'). Case-insensitive.
+        release (str): The Legacy Survey Data Release (e.g., 'DR9').
+
+    Returns:
+        astropy.io.fits.HDUList: The opened FITS file object.
+    """
+
+    if sia is None:
+        raise ImportError(
+            "Cannot use ls_open without pyvo. Please install pyvo (pip install pyvo) before continuing."
+        )
+
+    # 1. Set the specific SIA service endpoint for the desired release
+    # SIA endpoints for specific surveys are listed in the notebook.
+    service_url = f"https://datalab.noirlab.edu/sia/{release.lower()}"
+    svc = sia.SIAService(service_url)
+
+    # 2. Convert size from arcseconds to degrees (FOV) for the SIA query
+    # and apply the cosine correction for RA.
+    fov_deg = size_arcsec / 3600.0
+
+    # The search method takes the position (RA, Dec) and the square FOV.
+    imgTable = svc.search(
+        (ra, dec), (fov_deg / np.cos(dec * np.pi / 180.0), fov_deg), verbosity=2
+    ).to_table()
+
+    # 3. Filter the table for stacked images in the specified band
+    target_band = band.lower()
+
+    sel = (
+        (imgTable["proctype"] == "Stack")
+        & (imgTable["prodtype"] == "image")
+        & (startswith(imgTable["obs_bandpass"].astype(str), target_band))
+    )
+
+    Table = imgTable[sel]
+
+    if len(Table) == 0:
+        raise ValueError(
+            f"No stacked FITS image found for {release} band '{band}' at the requested RA {ra} and Dec {dec}."
+        )
+
+    # 4. Pick the "deepest" image (longest exposure time)
+    # Note: 'exptime' data needs explicit float conversion for np.argmax
+    max_exptime_index = np.argmax(Table["exptime"].data.data.astype("float"))
+    row = Table[max_exptime_index]
+
+    # 5. Download the file and open it
+    url = row["access_url"]  # get the download URL
+
+    # Use astropy's download_file, which handles the large data transfer
+    # and automatically uses a long timeout (120s in the notebook example)
+    filename = download_file(url, cache=False, show_progress=False, timeout=120)
+
+    # Open the downloaded FITS file
+    hdu = fits.open(filename)
+
+    try:
+        hdu = flip_hdu(hdu)
+    except AssertionError:
+        pass  # If WCS info is missing, skip flipping
+
+    # Clean up the temporary file created by download_file
+    try:
+        os.remove(filename)
+    except OSError:
+        pass  # Ignore if cleanup fails
+
+    return hdu

astrophot/utils/initialize/segmentation_map.py

@@ -56,7 +56,7 @@ def centroids_from_segmentation_map(
     if sky_level is None:
         sky_level = np.nanmedian(backend.to_numpy(image.data))
 
-    data = backend.to_numpy(image.data) - sky_level
+    data = backend.to_numpy(image._data) - sky_level
     centroids = {}
 
     II, JJ = np.meshgrid(np.arange(seg_map.shape[0]), np.arange(seg_map.shape[1]), indexing="ij")
@@ -94,7 +94,7 @@ def PA_from_segmentation_map(
     if sky_level is None:
         sky_level = np.nanmedian(backend.to_numpy(image.data))
 
-    data = backend.to_numpy(image.data) - sky_level
+    data = backend.to_numpy(image._data) - sky_level
 
     if centroids is None:
         centroids = centroids_from_segmentation_map(
@@ -141,7 +141,7 @@ def q_from_segmentation_map(
     if sky_level is None:
         sky_level = np.nanmedian(backend.to_numpy(image.data))
 
-    data = backend.to_numpy(image.data) - sky_level
+    data = backend.to_numpy(image._data) - sky_level
 
     if centroids is None:
         centroids = centroids_from_segmentation_map(
@@ -232,8 +232,8 @@ def scale_windows(windows, image: "Image" = None, expand_scale=1.0, expand_borde
             new_window = [
                 [max(0, new_window[0][0]), max(0, new_window[0][1])],
                 [
-                    min(image.data.shape[0], new_window[1][0]),
-                    min(image.data.shape[1], new_window[1][1]),
+                    min(image._data.shape[0], new_window[1][0]),
+                    min(image._data.shape[1], new_window[1][1]),
                 ],
             ]
         new_windows[index] = new_window
@@ -296,7 +296,7 @@ def filter_windows(
         if min_flux is not None:
             if (
                 np.sum(
-                    backend.to_numpy(image.data)[
+                    backend.to_numpy(image._data)[
                         windows[w][0][0] : windows[w][1][0],
                         windows[w][0][1] : windows[w][1][1],
                     ]
@@ -307,7 +307,7 @@ def filter_windows(
         if max_flux is not None:
             if (
                 np.sum(
-                    backend.to_numpy(image.data)[
+                    backend.to_numpy(image._data)[
                         windows[w][0][0] : windows[w][1][0],
                         windows[w][0][1] : windows[w][1][1],
                     ]

docs/requirements.txt

@@ -3,12 +3,13 @@ corner
 emcee
 graphviz
 ipywidgets
-jax
+jax<=0.7.0
 jupyter-book<2.0
 matplotlib
 nbformat
 nbsphinx
 photutils
+pyvo
 scikit-image
 sphinx
 sphinx-rtd-theme

docs/source/tutorials/GroupModels.ipynb

@@ -37,11 +37,8 @@
    "outputs": [],
    "source": [
     "# first let's download an image to play with\n",
-    "hdu = fits.open(\n",
-    "    \"https://www.legacysurvey.org/viewer/fits-cutout?ra=155.7720&dec=15.1494&size=150&layer=ls-dr9&pixscale=0.262&bands=r\"\n",
-    ")\n",
+    "hdu = ap.utils.ls_open(155.7720, 15.1494, 150 * 0.262, band=\"r\")\n",
     "target_data = np.array(hdu[0].data, dtype=np.float64)\n",
-    "\n",
     "fig1, ax1 = plt.subplots(figsize=(8, 8))\n",
     "plt.imshow(np.arctan(target_data / 0.05), origin=\"lower\", cmap=\"inferno\")\n",
     "plt.axis(\"off\")\n",
@@ -61,7 +58,7 @@
     "#########################################\n",
     "from photutils.segmentation import detect_sources, deblend_sources\n",
     "\n",
-    "initsegmap = detect_sources(target_data, threshold=0.02, npixels=5)\n",
+    "initsegmap = detect_sources(target_data, threshold=0.02, npixels=6)\n",
     "segmap = deblend_sources(target_data, initsegmap, npixels=5).data\n",
     "fig8, ax8 = plt.subplots(figsize=(8, 8))\n",
     "ax8.imshow(segmap, origin=\"lower\", cmap=\"inferno\")\n",

pyproject.toml

@@ -39,7 +39,7 @@ Repository = "https://github.com/Autostronomy/AstroPhot"
 Issues = "https://github.com/Autostronomy/AstroPhot/issues"
 
 [project.optional-dependencies]
-dev = ["pre-commit", "nbval", "nbconvert", "graphviz", "ipywidgets", "jupyter-book", "matplotlib", "photutils", "scikit-image", "caustics", "emcee", "corner", "jax"]
+dev = ["pre-commit", "nbval", "nbconvert", "graphviz", "ipywidgets", "jupyter-book", "matplotlib", "photutils", "scikit-image", "caustics", "emcee", "corner", "jax<=0.7.0", "pyvo"]
 
 [project.scripts]
 astrophot = "astrophot:run_from_terminal"