Merge branch 'dev'

Author: temuller

README.md

@@ -170,13 +170,15 @@ If you make use of HostPhot, please cite the following [paper](https://joss.theo
 v3.0.0
 * Major internal restructuring of the code
     * Convert the cutouts into a module with each survey as a different submodule
-    * Storage of input parameters
+    * Storage of input parameters (with timestamps)
     * Storage of output parameters in a more friendly format
     * Restructuring of output files
 * Fix the VISTA SCIENCE ARCHIVE url for downloading VISTA images
+* Fix issue with VISTA image size when downloading cuouts
 * The download of VISTA images now automatically chooses the survey with overlap for the given coordinates, prioritising the survey with deeper coverage.
 * The user experience is kept approximately consistent with previous versions, but there are some changes
-
+* Download of HST images is now working!
+* Adding download of JWST images!
 ___
 v2.10.0
 * Adding JWST! (few NIRCam filters for now - thanks to Lisa Kelsey for this addition)

src/hostphot/_version.py

@@ -1,3 +1,3 @@
 # -*- coding: utf-8 -*-
 
-__version__ = "3.0.0rc0"
+__version__ = "3.0.0rc1"

src/hostphot/cutouts/downloads.py

@@ -5,7 +5,7 @@
 from typing import Optional
 
 from hostphot._constants import workdir
-from hostphot.utils import check_work_dir
+from hostphot.utils import check_work_dir, store_input
 from hostphot.surveys_utils import (
     get_survey_filters,
     check_survey_validity,
@@ -64,8 +64,8 @@ def download_images(
 
     # save input parameters
     if save_input is True:
-        inputs_df = pd.DataFrame({key: [value] for key, value in input_params.items()})
-        inputs_df.to_csv(survey_dir / "input_cutouts.csv", index=False)
+        inputs_file = survey_dir / "input_cutouts.csv"
+        store_input(input_params, inputs_file)
 
     # check existing images
     if overwrite is False:
@@ -105,8 +105,9 @@ def download_images(
             continue  # skip missing filter/image
         # get output file name
         if survey == "HST":
-            inst = version.replace("/", "_")
-            outfile = Path(survey_dir, f"{survey}_{inst}_{filters}.fits")
+            #inst = version.replace("/", "_")
+            #outfile = Path(survey_dir, f"{survey}_{inst}_{filters}.fits")
+            outfile = Path(survey_dir, f"{survey}_{filt}.fits")
         else:
             outfile = Path(survey_dir, f"{survey}_{filt}.fits")
         if overwrite is True or outfile.is_file() is False:

src/hostphot/cutouts/hst.py

@@ -1,20 +1,17 @@
-import shutil
-import zipfile
 import numpy as np
-import pandas as pd
 from pathlib import Path
 
-from astropy import wcs
 import astropy.units as u
 from astropy.io import fits
+from astropy.wcs import WCS
+from astropy.nddata import Cutout2D
 from astropy.coordinates import SkyCoord
 
-from pyvo.dal import sia
 from astroquery.esa.hubble import ESAHubble  # HST
 
 from hostphot._constants import workdir
-from hostphot.utils import check_work_dir
-from hostphot.surveys_utils import check_HST_filters
+from hostphot.utils import check_work_dir, suppress_stdout
+from hostphot.surveys_utils import check_HST_filters, survey_pixel_scale
 
 import warnings
 from astropy.utils.exceptions import AstropyWarning
@@ -34,7 +31,7 @@ def update_HST_header(hdu: fits.hdu.ImageHDU) -> None:
         warnings.simplefilter("ignore", AstropyWarning)
         if "PHOTPLAM" not in hdu[0].header:
             # MAST image: move things to hdu[0] to homogenise
-            img_wcs = wcs.WCS(hdu[1].header)
+            img_wcs = WCS(hdu[1].header)
             hdu[0].header.update(img_wcs.to_header())
             hdu[0].header["PHOTFLAM"] = hdu[1].header["PHOTFLAM"]
             hdu[0].header["PHOTPLAM"] = hdu[1].header["PHOTPLAM"]
@@ -66,6 +63,7 @@ def set_HST_image(file: str, filt: str, name: str) -> None:
     """
     # check output directory
     check_work_dir(workdir)
+    check_HST_filters(filt)
     obj_dir = Path(workdir, name, "HST")
     if obj_dir.is_dir() is False:
         obj_dir.mkdir(parents=True, exist_ok=True)
@@ -76,36 +74,24 @@ def set_HST_image(file: str, filt: str, name: str) -> None:
     hdu.writeto(outfile, overwrite=True)
 
 def get_HST_images(ra: float, dec: float, size: float | u.Quantity = 3, 
-                        filt: str = "WFC3_UVIS_F225W") -> list[fits.ImageHDU]:
+                        filters: list = ["WFC3_UVIS_F225W"]) -> list[fits.ImageHDU]:
     """Downloads a set of HST fits images for a given set
-    of coordinates and filters using the MAST archive.
+    of coordinates and filters using astroquery.
 
     Parameters
     ----------
     ra: Right ascension in degrees.
     dec: Declination in degrees.
     size: Image size. If a float is given, the units are assumed to be arcmin.
-    filt: Filter to use, e.g. ``WFC3_UVIS_F225W``.
+    filters: Filters to use, e.g. ``WFC3_UVIS_F225W``.
 
     Return
     ------
     hdu_list: List with fits image for the given filter. ``None`` is returned if no image is found.
     """
     esahubble = ESAHubble()
     esahubble.get_status_messages()
-    check_HST_filters(filt)
-
-    # separate the instrument name from the actual filter
-    split_filt = filt.split("_")
-    if len(split_filt) == 2:
-        filt = split_filt[-1]
-        instrument = split_filt[0]
-    elif len(split_filt) == 3:
-        filt = split_filt[-1]
-        instrument = f"{split_filt[0]}/{split_filt[1]}"
-    else:
-        raise ValueError(f"Incorrect filter name: {filt}")
-
+    
     if isinstance(size, (float, int)):
         size_arcsec = (size * u.arcmin).to(u.arcsec)
     else:
@@ -116,43 +102,43 @@ def get_HST_images(ra: float, dec: float, size: float | u.Quantity = 3,
         coords = SkyCoord(
             ra=ra, dec=dec, unit=(u.degree, u.degree), frame="icrs"
         )
-
-    version = None
-    if version == "HLA":
-        # This does not seem to be faster
-        fov = 0.2  # field-of-view in degrees
-        access_url = " https://hla.stsci.edu/cgi-bin/hlaSIAP.cgi"
-        svc = sia.SIAService(access_url)
-        imgs_table = svc.search(
-            (ra, dec), (fov / np.cos(dec * np.pi / 180), fov), verbosity=2
-        )
-        obs_df = pd.DataFrame(imgs_table)
-        obs_df = obs_df[obs_df.Mode == "IMAGE"]
-        obs_df = obs_df[obs_df.Format.str.endswith("fits")]
-        obs_df = obs_df[obs_df.Detector == instrument]
-        obs_df = obs_df[obs_df.Spectral_Elt == filt]
-        obs_df = obs_df[obs_df.ExpTime == obs_df.ExpTime.max()]
-        hdu = fits.open(obs_df.URL.values[0])
-    else:
+    
+    # query observations at the given coordinates
+    with suppress_stdout():
         result = esahubble.cone_search_criteria(
             radius=3,
             coordinates=coords,
             calibration_level="PRODUCT",
             data_product_type="image",
-            instrument_name=instrument,
-            filters=filt,
+            #instrument_name=instrument,
+            #filters=filt,
             async_job=True,
         )
-
-        obs_df = result.to_pandas()
-        obs_df = obs_df[obs_df["filter"] == filt]
+    results_df = result.to_pandas()
+    
+    hdu_list = []
+    for filt in filters:
+        check_HST_filters(filt)
+        # separate the instrument name from the actual filter
+        split_filt = filt.split("_")
+        if len(split_filt) == 2:
+            filt = split_filt[-1]
+            instrument = split_filt[0]
+        elif len(split_filt) == 3:
+            filt = split_filt[-1]
+            instrument = f"{split_filt[0]}/{split_filt[1]}"
+        else:
+            raise ValueError(f"Incorrect filter name: {filt}")
+        # filter by filter and instrument
+        obs_df = results_df[results_df["filter"] == filt]
+        obs_df = obs_df[obs_df.instrument_name == instrument]
         # get only exposures shorter than one hour
         obs_df = obs_df[obs_df.exposure_duration < 3600]
         obs_df.sort_values(
             by=["exposure_duration"], ascending=False, inplace=True
         )
 
-        print("Looking for HST images...")
+        # start download 
         filename = f"HST_tmp_{ra}_{dec}"  # the extension is added below
         for obs_id in obs_df.observation_id:
             try:
@@ -167,14 +153,24 @@ def get_HST_images(ra: float, dec: float, size: float | u.Quantity = 3,
                 pass
         temp_file = Path(f"{filename}.fits.gz")
         if temp_file.is_file() is False:
-            return None
-        temp_dir = Path(filename)  # same name as the extensionless file above
-        with zipfile.ZipFile(temp_file, "r") as zip_ref:
-            zip_ref.extractall(temp_dir)
-            fits_file = [file for file in temp_dir.rglob("*.gz")][0]
-        hdu = fits.open(fits_file)
-        # remove the temporary files and directory
+            hdu_list.append(None)
+            continue
+        hdu = fits.open(temp_file)
+        # remove the temporary files
         temp_file.unlink()
-        shutil.rmtree(temp_dir, ignore_errors=True)
-    update_HST_header(hdu)
-    return [hdu]
+        # add necessary information to the header
+        update_HST_header(hdu)
+        hdu_list.append(hdu)
+    # HST images are large so need to be trimmed
+    for hdu, filt in zip(hdu_list, filters):
+        if hdu is None:
+            continue
+        pixel_scale = survey_pixel_scale("HST", filt)  # same pixel scale for all filters
+        size_pixels = int(size_arcsec / pixel_scale)
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore", AstropyWarning)
+            img_wcs = WCS(hdu[0].header)
+        trimmed_data = Cutout2D(hdu[0].data, coords, size_pixels, img_wcs)
+        hdu[0].data = trimmed_data.data
+        hdu[0].header.update(trimmed_data.wcs.to_header())
+    return hdu_list

src/hostphot/cutouts/jwst.py

@@ -1,12 +1,17 @@
 import numpy as np
 from pathlib import Path
 
-from astropy import wcs
+import astropy.units as u
 from astropy.io import fits
+from astropy.wcs import WCS
+from astropy.nddata import Cutout2D
+from astropy.coordinates import SkyCoord
+
+from astroquery.esa.jwst import Jwst
 
 from hostphot._constants import workdir
-from hostphot.utils import check_work_dir
-from hostphot.surveys_utils import check_JWST_filters
+from hostphot.utils import check_work_dir, suppress_stdout
+from hostphot.surveys_utils import check_JWST_filters, survey_pixel_scale
 
 import warnings
 from astropy.utils.exceptions import AstropyWarning
@@ -22,7 +27,7 @@ def update_JWST_header(hdu: fits.hdu.ImageHDU) -> None:
     # get WCS - SCI is hdu[1] - hostphot always assumes hdu[0] is used, so need to move them over
     with warnings.catch_warnings():
         warnings.simplefilter("ignore", AstropyWarning)
-        img_wcs = wcs.WCS(hdu[1].header)
+        img_wcs = WCS(hdu[1].header)
     hdu[0].header.update(img_wcs.to_header())
     hdu[0].header["PIXAR_SR"] = hdu[1].header["PIXAR_SR"]
     hdu[0].data = hdu[1].data
@@ -52,11 +57,101 @@ def set_JWST_image(file: str, filt: str, name: str) -> None:
     """
     # check output directory
     check_work_dir(workdir)
+    check_JWST_filters(filt)
     obj_dir = Path(workdir, name, "JWST")
     if obj_dir.is_dir() is False:
         obj_dir.mkdir(parents=True, exist_ok=True)
     # update header and save file
     hdu = fits.open(file)
     update_JWST_header(hdu)
     outfile = obj_dir / f"JWST_{filt}.fits"
-    hdu.writeto(outfile, overwrite=True)
+    hdu.writeto(outfile, overwrite=True)
+    
+def get_JWST_images(ra: float, dec: float, size: float | u.Quantity = 3, 
+                        filters: list = ["WFC3_UVIS_F225W"]) -> list[fits.ImageHDU]:
+    """Downloads a set of JWST fits images for a given set
+    of coordinates and filters using astroquery.
+
+    Parameters
+    ----------
+    ra: Right ascension in degrees.
+    dec: Declination in degrees.
+    size: Image size. If a float is given, the units are assumed to be arcmin.
+    filters: Filters to use, e.g. ``WFC3_UVIS_F225W``.
+
+    Return
+    ------
+    hdu_list: List with fits image for the given filter. ``None`` is returned if no image is found.
+    """
+    Jwst.get_status_messages()
+    
+    if isinstance(size, (float, int)):
+        size_arcsec = (size * u.arcmin).to(u.arcsec)
+    else:
+        size_arcsec = size.to(u.arcsec)
+    size_arcsec = size_arcsec.value
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", AstropyWarning)
+        coords = SkyCoord(
+            ra=ra, dec=dec, unit=(u.degree, u.degree), frame="icrs"
+        )
+
+    # query observations at the given coordinates
+    with suppress_stdout():
+        result = Jwst.cone_search(
+            radius=3*u.arcsec,
+            coordinate=coords,
+            cal_level="Top",
+            prod_type="image",
+            #instrument_name=instrument,
+            #filter_name=filt,
+            only_public=True,
+            async_job=True,
+        ).get_data()
+    results_df = result.to_pandas()
+            
+    hdu_list = []
+    for filt in filters:
+        check_JWST_filters(filt)
+        # separate the instrument name from the actual filter
+        split_filt = filt.split("_")
+        filt = split_filt[-1]
+        instrument = split_filt[0].upper() + "/IMAGE"
+        # filter by filter and instrument
+        obs_df = results_df[results_df["instrument_name"] == instrument]
+        obs_df = obs_df[obs_df["energy_bandpassname"] == filt]
+        obs_df = obs_df[obs_df["calibrationlevel"] == 3]
+
+        # start download
+        temp_file = None
+        for obs_id in obs_df.observationid:
+            try:
+                product_list = Jwst.get_product_list(observation_id=obs_id, product_type='science', cal_level=3).to_pandas()
+                # choose first file as there es no exptime info and the images take too lonk to download multiple ones
+                file_name = [file for file in product_list.filename if file.endswith(".fits")][0] 
+                temp_file = Jwst.get_product(file_name=file_name)
+                break
+            except:
+                pass
+        if temp_file is None:
+            hdu_list.append(None)
+            continue
+        hdu = fits.open(temp_file)
+        # remove the temporary files
+        Path(temp_file).unlink()
+        # add necessary information to the header
+        update_JWST_header(hdu)
+        hdu_list.append(hdu)
+    # JWST images are large so need to be trimmed
+    for hdu, filt in zip(hdu_list, filters):
+        if hdu is None:
+            continue
+        pixel_scale = survey_pixel_scale("JWST", filt)  # same pixel scale for all filters
+        size_pixels = int(size_arcsec / pixel_scale)
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore", AstropyWarning)
+            img_wcs = WCS(hdu[0].header)
+        trimmed_data = Cutout2D(hdu[0].data, coords, size_pixels, img_wcs)
+        hdu[0].data = trimmed_data.data
+        hdu[0].header.update(trimmed_data.wcs.to_header())
+    return hdu_list