diff --git a/astrophot/image/image_object.py b/astrophot/image/image_object.py index b8b05e26..6926877a 100644 --- a/astrophot/image/image_object.py +++ b/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 = ( diff --git a/astrophot/models/func/gaussian_ellipsoid.py b/astrophot/models/func/gaussian_ellipsoid.py index 2a989f61..4b07e9cf 100644 --- a/astrophot/models/func/gaussian_ellipsoid.py +++ b/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, ) diff --git a/astrophot/models/gaussian_ellipsoid.py b/astrophot/models/gaussian_ellipsoid.py index 02cd54bd..2d14da51 100644 --- a/astrophot/models/gaussian_ellipsoid.py +++ b/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) diff --git a/astrophot/utils/__init__.py b/astrophot/utils/__init__.py index b66971a3..33925367 100644 --- a/astrophot/utils/__init__.py +++ b/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", ] diff --git a/astrophot/utils/fitsopen.py b/astrophot/utils/fitsopen.py new file mode 100644 index 00000000..5c9a8d70 --- /dev/null +++ b/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 diff --git a/astrophot/utils/initialize/segmentation_map.py b/astrophot/utils/initialize/segmentation_map.py index 6364ba40..053d257b 100644 --- a/astrophot/utils/initialize/segmentation_map.py +++ b/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], ] diff --git a/docs/requirements.txt b/docs/requirements.txt index d2b0ecc0..07b09906 100644 --- a/docs/requirements.txt +++ b/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 diff --git a/docs/source/tutorials/GroupModels.ipynb b/docs/source/tutorials/GroupModels.ipynb index b4a719eb..f442b811 100644 --- a/docs/source/tutorials/GroupModels.ipynb +++ b/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", diff --git a/pyproject.toml b/pyproject.toml index 7bf255b0..faaf81cf 100644 --- a/pyproject.toml +++ b/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"