fixing aperture photometry bug for DES images

CHANGELOG.md

@@ -13,6 +13,28 @@ Types of changes:
 - `Fixed`: for any bug fixes.
 - `Security`: in case of vulnerabilities.
 
+## [1.13.1]
+
+### Changed
+
+- Now using astropy version 7.2.0, prospector 1.4.1, photutils dev version (commit 2573ef4, until
+  latest release)
+- Included legacy survey data release [DES DR10](https://www.legacysurvey.org/dr10/description/)
+- Updated how legacy survey downloads works, resulting in fewer download timeouts.
+- Improved efficiency of cutout unit test by using ThreadPoolExecutor to test downloads in parallel.
+
+### Fixed
+
+- Incorporated photutils fix for misaligned apertures in DES legacy survey images that caused offset
+  apertures in the plotting and photometry.
+
+## [1.13.0]
+
+### Changed
+
+- Reintroduce the formerly excluded `host` value in `TransientSerializer` that defines the data
+  object returned by the `/api/transient/` API endpoint.
+
 ## [1.12.0]
 
 ### Added

app/app/settings.py

@@ -4,7 +4,7 @@
 ######################################################################
 # Blast application config
 #
-APP_VERSION = '1.12.0'
+APP_VERSION = '1.13.1'
 # Data paths
 DUSTMAPS_DATA_ROOT = os.environ.get("DUSTMAPS_DATA_ROOT", "/data/dustmaps")
 CUTOUT_ROOT = os.environ.get("CUTOUT_ROOT", "/data/cutout_cdn")

app/host/SBI/run_sbi_blast.py

@@ -232,13 +232,13 @@ def fit_sbi_pp(observations, n_filt_cuts=True, fit_type="global"):
     for f in all_filters:
         if f.name in observations["filternames"]:
             iflt = np.array(observations["filternames"]) == f.name
-            mags = np.append(mags, maggies_to_asinh(observations["maggies"][iflt]))
+            mags = np.append(mags, maggies_to_asinh(observations["maggies"][iflt][0]))
             mags_unc = np.append(
                 mags_unc,
                 2.5
                 / np.log(10)
-                * observations["maggies_unc"][iflt]
-                / observations["maggies"][iflt],
+                * observations["maggies_unc"][iflt][0]
+                / observations["maggies"][iflt][0],
             )
             if f.name in uv_filters:
                 has_uv = True

app/host/cutouts.py

@@ -449,10 +449,40 @@ def DES_cutout(position, image_size=None, filter=None):
     :cutout : :class:`~astropy.io.fits.HDUList` or None
     """
 
-    DEF_ACCESS_URL = "https://datalab.noirlab.edu/sia/ls_dr9"
-    svc_ls_dr9 = sia.SIAService(DEF_ACCESS_URL)
+    # Try DR10 first, fall back to DR9 if data is not returned.
+    for url in ["https://datalab.noirlab.edu/sia/ls_dr10",
+                "https://datalab.noirlab.edu/sia/ls_dr9"]:
+        fits_image = DES_cutout_single_version(
+            position, image_size=image_size, filter=filter,
+            access_url=url
+        )
+        if fits_image is not None:
+            break
+
+    return fits_image
+
+
+def DES_cutout_single_version(
+    position, image_size=None, filter=None,
+    access_url="https://datalab.noirlab.edu/sia/ls_dr10"
+):
+    """
+    Download DES image cutout from NOIRLab
+
+    Parameters
+    ----------
+    :position : :class:`~astropy.coordinates.SkyCoord`
+        Target centre position of the cutout image to be downloaded.
+    :image_size: int: size of cutout image in pixels
+    :filter: str: Panstarrs filter (g r i z y)
+    Returns
+    -------
+    :cutout : :class:`~astropy.io.fits.HDUList` or None
+    """
 
-    imgTable = svc_ls_dr9.search(
+    svc_ls = sia.SIAService(access_url)
+
+    imgTable = svc_ls.search(
         (position.ra.deg, position.dec.deg),
         (image_size / np.cos(position.dec.deg * np.pi / 180), image_size),
         verbosity=2,
@@ -464,44 +494,47 @@ def DES_cutout(position, image_size=None, filter=None):
             valid_urls += [img["access_url"]]
             logger.debug(f'''DES image URL: {valid_urls[-1]}''')
 
-    if len(valid_urls):
-        # we need both the depth and the image
-        time.sleep(1)
-        try:
-            fits_image = fits.open(
-                valid_urls[0].replace("-depth-", "-image-"), cache=None
-            )
-        except Exception as e:
-            ### found some bad links...
-            return None
-        if np.shape(fits_image[0].data)[0] == 1 or np.shape(fits_image[0].data)[1] == 1:
-            # no idea what's happening here but this is a mess
-            return None
-        try:
-            depth_image = fits.open(valid_urls[0])
-        except Exception as e:
-            # wonder if there's some issue with other tasks clearing the cache
-            time.sleep(5)
-            depth_image = fits.open(valid_urls[0])
-        wcs_depth = WCS(depth_image[0].header)
-        xc, yc = wcs_depth.wcs_world2pix(position.ra.deg, position.dec.deg, 0)
-
-        # this is ugly - just assuming the exposure time at the
-        # location of interest is uniform across the image
-        if np.shape(depth_image[0].data) == (1, 1):
-            exptime = depth_image[0].data[0][0]
-        else:
-            exptime = depth_image[0].data[int(yc), int(xc)]
-        if exptime == 0:
-            fits_image = None
-        else:
-            wcs = WCS(fits_image[0].header)
-            cutout = Cutout2D(fits_image[0].data, position, image_size, wcs=wcs)
-            fits_image[0].data = cutout.data
-            fits_image[0].header.update(cutout.wcs.to_header())
-            fits_image[0].header["EXPTIME"] = exptime
+    if not len(valid_urls):
+        return None
+    # we need both the depth and the image
+    time.sleep(1)
+    try:
+        req = requests.get(valid_urls[0].replace("-depth-", "-image-"), stream=True)
+        fits_image = fits.open(BytesIO(req.content))
+    except Exception:
+        print(f'opening the URL {valid_urls[0].replace("-depth-", "-image-")} failed')
+        # found some bad links...
+        return None
+
+    if np.shape(fits_image[0].data)[0] == 1 or np.shape(fits_image[0].data)[1] == 1:
+        # no idea what's happening here but this is a mess
+        return None
+    try:
+        req = requests.get(valid_urls[0], stream=True)
+        depth_image = fits.open(BytesIO(req.content))
+    except Exception:
+        # wonder if there's some issue with other tasks clearing the cache
+        time.sleep(5)
+        req = requests.get(valid_urls[0], stream=True)
+        depth_image = fits.open(BytesIO(req.content))
+
+    wcs_depth = WCS(depth_image[0].header)
+    xc, yc = wcs_depth.wcs_world2pix(position.ra.deg, position.dec.deg, 0)
+
+    # this is ugly - just assuming the exposure time at the
+    # location of interest is uniform across the image
+    if np.shape(depth_image[0].data) == (1, 1):
+        exptime = depth_image[0].data[0][0]
     else:
+        exptime = depth_image[0].data[int(yc), int(xc)]
+    if exptime == 0:
         fits_image = None
+    else:
+        wcs = WCS(fits_image[0].header)
+        cutout = Cutout2D(fits_image[0].data, position, image_size, wcs=wcs)
+        fits_image[0].data = cutout.data
+        fits_image[0].header.update(cutout.wcs.to_header())
+        fits_image[0].header["EXPTIME"] = exptime
 
     return fits_image
 
@@ -521,7 +554,7 @@ def TWOMASS_cutout(position, image_size=None, filter=None):
     :cutout : :class:`~astropy.io.fits.HDUList` or None
     """
 
-    irsaquery = f"https://irsa.ipac.caltech.edu/cgi-bin/2MASS/IM/nph-im_sia?POS={position.ra.deg},{position.dec.deg}&SIZE=0.01"
+    irsaquery = f"https://irsa.ipac.caltech.edu/cgi-bin/2MASS/IM/nph-im_sia?POS={position.ra.deg},{position.dec.deg}&SIZE=0.01"  # noqa
     response = requests.get(url=irsaquery)
 
     fits_image = None

app/host/host_utils.py

@@ -430,7 +430,7 @@ def check_global_contamination(global_aperture_phot, aperture_primary):
                 .to_mask()
                 .to_image(np.shape(image[0].data))
             )
-            obj_ids = catalog._segment_img.data[np.where(mask_image == True)]  # noqa: E712
+            obj_ids = catalog._segmentation_image.data[np.where(mask_image == True)]  # noqa: E712
             source_obj = source_data._labels
 
             # let's look for contaminants
@@ -555,7 +555,6 @@ def get_source_data(threshhold_sigma):
     # make sure we know this failed
     if source_separation_arcsec > 5:
         return None
-
     return elliptical_sky_aperture(source_data, wcs)
 
 

app/host/plotting_utils.py

@@ -93,14 +93,14 @@ def plot_position(object, wcs, plotting_kwargs=None, plotting_func=None):
 
 def plot_aperture(figure, aperture, wcs, plotting_kwargs=None):
     aperture = aperture.to_pixel(wcs)
-    theta_rad = aperture.theta
+    theta_rad = aperture.theta.rad
     x, y = aperture.positions
     plot_dict = {
         "x": x,
         "y": y,
         "width": aperture.a * 2,
         "height": aperture.b * 2,
-        "angle": theta_rad.value,
+        "angle": theta_rad,
         "fill_color": "#cab2d6",
         "fill_alpha": 0.1,
         "line_width": 4,

app/host/tests/test_cutouts.py

@@ -4,39 +4,67 @@
 from astropy.io import fits
 from django.test import TestCase
 from django.test import tag
+from concurrent.futures import ThreadPoolExecutor, as_completed
 
 from ..cutouts import cutout
 from ..models import Filter
 
+from host.log import get_logger
+logger = get_logger(__name__)
+
 
 class CutoutDownloadTest(TestCase):
 
     def setUp(self):
-        self.transient_name = "2010ag"
-        self.transient_ra = 355.53628555
-        self.transient_dec = 48.70907059166666
+        self.transient_name = "2026dgt"
+        self.transient_ra = 132.3563
+        self.transient_dec = 29.5105
 
     @tag('download')
     def test_cutout_download(self):
         """ "
         Test that cutout data can be downloaded.
         """
+        status_failed = "failed"
+        status_succeeded = "processed"
         download_dir = '/tmp/'
         os.makedirs(download_dir, exist_ok=True)
-        for filter in Filter.objects.all():
+
+        def process_filter(filter):
             print(f'''Processing filter "{filter.name}"...''')
-            position = SkyCoord(
-                ra=self.transient_ra, dec=self.transient_dec, unit="deg"
-            )
+            position = SkyCoord(ra=self.transient_ra, dec=self.transient_dec, unit="deg")
             save_dir = os.path.join(download_dir, f"{self.transient_name}/{filter.name}/")
             os.makedirs(save_dir, exist_ok=True)
             path_to_fits = os.path.join(save_dir, f"{filter.name}.fits")
-            if os.path.exists(path_to_fits):
-                print(f'''Skipping file already downloaded: "{path_to_fits}"...''')
-                continue
+            try:
+                # Delete file if for some reason it already exists
+                os.remove(path_to_fits)
+            except FileNotFoundError:
+                pass
+            cutout_data, status, error = cutout(position, filter)
+            if cutout_data:
+                logger.info(f'[{filter}] Cutout data downloaded.')
+                return status_succeeded
+            elif status == 0 and not error:
+                logger.info(f'[{filter}] No cutout data available, but no errors.')
+                return status_succeeded
             else:
-                cutout_data = cutout(position, filter)[0]
-                if cutout_data:
-                    fits.writeto(path_to_fits, cutout_data[0].data, overwrite=True)
+                logger.error(f'[{filter}] Error downloading cutout data: {error}')
+                return status_failed
+
+        results = []
+        filter_set = Filter.objects.all()
+        with ThreadPoolExecutor(max_workers=len(filter_set)) as executor:
+            future_to_filter = {executor.submit(process_filter, filter): filter for filter in filter_set}
+            for future in as_completed(future_to_filter):
+                filter = future_to_filter[future]
+                try:
+                    result = future.result()
+                except Exception as exc:
+                    logger.error('%r generated an exception: %s' % (filter, exc))
+                    result = status_failed
+                logger.debug(f'''"{filter}" result: "{result}"''')
+                results.append(result)
 
-        self.assertTrue(1 == 1)
+        # If any of the cutout downloads failed, the test fails.
+        self.assertFalse(not results or [result for result in results if result != status_succeeded])

app/host/transient_tasks.py

@@ -518,7 +518,7 @@ def _run_process(self, transient):
         data = {
             "name": f"{aperture_cutout[0].name}_global",
             "cutout": aperture_cutout[0],
-            "orientation_deg": (180 / np.pi) * aperture.theta.value,
+            "orientation_deg": (180 / np.pi) * aperture.theta.rad,
             "ra_deg": aperture.positions.ra.degree,
             "dec_deg": aperture.positions.dec.degree,
             "semi_major_axis_arcsec": aperture.a.value,

app/requirements.txt

@@ -1,7 +1,7 @@
-astro-datalab==2.24.0
-astro-prospector==1.4.0
+astro-datalab==2.22.1
+astro-prospector==1.4.1
 astro-prost==1.2.13
-astro-sedpy==0.4.0
+astro-sedpy==0.4.1
 astroquery==0.4.7
 bokeh==3.7.3
 coverage==7.9.0
@@ -26,10 +26,12 @@ minio==7.2.15
 mozilla-django-oidc==4.0.1
 mysqlclient==2.2.7  # TODO: This can be removed if using PostgreSQL
 psycopg2-binary==2.9.11
-photutils==2.2.0
+git+https://github.com/astropy/photutils.git@2573ef4
 redis==7.0.1
 sbi==0.22.0
 scikit-image==0.25.2
 selenium==4.35.0
 watchdog==6.0.0
 influxdb-client==1.50.0
+pandas==2.3.3
+ipython==9.14.0