diff --git a/CHANGELOG.md b/CHANGELOG.md index 32e1e2b..47c5125 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Added final SPICE kernels for the WISE mission, it now contains all positions from all phases of operation. There is a gap for the years it was not operating. +- Updated WISE mission phases to reflect the final data products about to be released. +- Updated ZTF for the current release 22. - Added `kete.RectangleFOV.from_wcs`, allowing the construction of a FOV from a given Astropy WCS object. - Added `kete.conversion.bin_data`, which allows for binning matrix data such as images. @@ -21,11 +23,22 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed +- Time scaling bugs when loading times from both Horizons and the MPC were fixed, these + were causing offsets of about a minute in the loaded states. However it was shifting + both the perihelion time and the epoch time by the same amount, leading to only minor + effective errors. - Fixed a time offset in the FOV's downloaded from IRSA WISE/NEOWISE. They were offset by 4.4 seconds. - Fixed rotation approximation in WISE/NEOWISE field of views which was causing a small percentage of objects to not be found during FOV checks when they were close to the edge of the field. +- Constant for the sqrt of GMS was incorrect by a small amount, this value was fixed. +- IRSA username/password options are now being passed through correctly in WISE. + +### Removed + +- Removed `plot_frame` from ZTF, as a better version of this is available in kete.irsa. +- Removed `cache_WISE_frame` and `fetch_WISE_frame` deprecated functions in WISE. ## [0.3.0] - 2024 - 8 - 28 diff --git a/src/examples/plot_comet.py b/src/examples/plot_comet.py index c5d13da..dd55c7b 100644 --- a/src/examples/plot_comet.py +++ b/src/examples/plot_comet.py @@ -173,7 +173,8 @@ def plot_synchrone( # Plot the final results plt.figure(dpi=200) -wcs = kete.ztf.plot_frame(vis.fov) +frame = kete.ztf.fetch_frame(vis.fov) +wcs = kete.irsa.plot_fits_image(frame, percentiles=None) plt.title("Comet NEOWISE - C/2020 F3\n") # plot syndynes diff --git a/src/kete/horizons.py b/src/kete/horizons.py index 286a830..4011972 100644 --- a/src/kete/horizons.py +++ b/src/kete/horizons.py @@ -91,8 +91,12 @@ def fetch(name, update_name=True, cache=True, update_cache=False, exact_name=Fal ] if len(val) == 0: continue - phys[kete_v] = val[0] - phys["epoch"] = float(props["orbit"]["epoch"]) + if kete_v == "peri_time": + phys[kete_v] = Time(val[0], scaling="utc").jd + else: + phys[kete_v] = val[0] + + phys["epoch"] = Time(float(props["orbit"]["epoch"]), scaling="utc").jd if "moid" in props["orbit"] and props["orbit"]["moid"] is not None: phys["moid"] = float(props["orbit"]["moid"]) if "data_arc" in props["orbit"] and props["orbit"]["data_arc"] is not None: diff --git a/src/kete/irsa.py b/src/kete/irsa.py index 53f94f5..422cf4c 100644 --- a/src/kete/irsa.py +++ b/src/kete/irsa.py @@ -146,7 +146,7 @@ def query_irsa_tap( phase_url = url.replace("results/result", "phase") - status = requests.get(phase_url, timeout=timeout) + status = requests.get(phase_url, timeout=timeout, auth=auth) status.raise_for_status() # Status results can have one of 4 outcomes: @@ -164,7 +164,7 @@ def query_irsa_tap( status.content.decode(), ) time.sleep(delay) - status = requests.get(phase_url, timeout=timeout) + status = requests.get(phase_url, timeout=timeout, auth=auth) status.raise_for_status() # Increase time between queries until there is 30 seconds between. @@ -175,7 +175,7 @@ def query_irsa_tap( if status.content.decode().upper() != "COMPLETED": raise ValueError("Job Failed: ", status.content.decode()) - result = requests.get(url, timeout=timeout) + result = requests.get(url, timeout=timeout, auth=auth) result.raise_for_status() return pd.read_csv(io.StringIO(result.text)) diff --git a/src/kete/mpc.py b/src/kete/mpc.py index 049d926..06d5f57 100644 --- a/src/kete/mpc.py +++ b/src/kete/mpc.py @@ -673,8 +673,8 @@ def fetch_known_orbit_data(url=None, force_download=False): incl=obj["i"], lon_node=obj["Node"], peri_arg=obj["Peri"], - peri_time=obj["Tp"], - epoch=obj["Epoch"], + peri_time=Time(obj["Tp"], scaling="utc").jd, + epoch=Time(obj["Epoch"], scaling="utc").jd, arc_len=arc_len, name=obj.get("Name", None), ) diff --git a/src/kete/wise.py b/src/kete/wise.py index c157fdc..9e3da8e 100644 --- a/src/kete/wise.py +++ b/src/kete/wise.py @@ -2,8 +2,9 @@ import os +import logging from collections import namedtuple -from functools import lru_cache, partial +from functools import lru_cache from typing import Optional, Union import matplotlib.pyplot as plt import numpy as np @@ -15,7 +16,6 @@ from .time import Time from .vector import Vector, Frames from .irsa import IRSA_URL, query_irsa_tap, plot_fits_image, zoom_plot, annotate_plot -from .deprecation import rename from ._core import ( WiseCmos, @@ -27,11 +27,11 @@ ) __all__ = [ + "fetch_frame", + "plot_frames", "MISSION_PHASES", "mission_phase_from_jd", "mission_phase_from_scan", - "plot_frames", - "fetch_WISE_frame", "MissionPhase", "w1_color_correction", "w2_color_correction", @@ -39,6 +39,8 @@ "w4_color_correction", ] +logger = logging.getLogger(__name__) + # All constants below are indexed as follows W1 = [0], W2 = [1]. W3 = [2], W4 = [3] _COLOR_CORR = np.array( @@ -88,29 +90,20 @@ # W4 is almost a constant. # These fits perform much better than linear interpolation except near 100k. # -# The values above were fit to polynomial equations as defined below, the results of -# which were hard-coded into the rust backend to allow for fast computation. -# These functions should not be referenced directly, as they also exist in the rust -# code and are only left here for reference to how the rust was constructed. -# -# from numpy.polynomial import Polynomial -# _COLOR_FITS = [ -# Polynomial.fit(_COLOR_CORR[:, 0], 1 / _COLOR_CORR[:, 1], 4), -# Polynomial.fit(_COLOR_CORR[:, 0], 1 / _COLOR_CORR[:, 2], 4), -# Polynomial.fit(_COLOR_CORR[:, 0], 1 / _COLOR_CORR[:, 3], 4), -# Polynomial.fit(_COLOR_CORR[:, 0], 1 / _COLOR_CORR[:, 4], 4), -# ] +# The values above were fit to polynomial equations, the results of which were +# hard-coded into the rust backend to allow for fast computation. SUN_COLOR_CORRECTION: list[float] = [1.0049, 1.0193, 1.0024, 1.0012] """ Flux in the reflected light model should be scaled by these values. + +This corrects for the Sun's spectral difference to the calibrator. """ ZERO_MAGS: list[float] = [306.681, 170.663, 29.0448, 8.2839] """ Non-color corrected values for zero mag corrections in Janskys. -Magnitude can then be computed via -2.5 log10(flux Jy / zero_point) """ ZERO_MAGS_COLOR_CORRECTED: list[float] = [ @@ -119,7 +112,7 @@ 1.08 * ZERO_MAGS[2], 0.96 * ZERO_MAGS[3], ] -"""Color Corrected Zero Mags in units of Jy""" +"""Color Corrected Zero Mags in units of Jy.""" BAND_WAVELENGTHS: list[float] = [3352.6, 4602.8, 11560.8, 22088.3] """Non-color corrected values for the effective central wavelength of the bands (nm)""" @@ -212,10 +205,10 @@ frame_meta_table="neowiser_p1bs_frm", source_table="neowiser_p1bs_psd", ), - "Reactivation_2023": MissionPhase( - name="Reactivation_2023", - jd_start=Time.from_ymd(2023, 1, 1).jd, - jd_end=Time.from_ymd(2023, 12, 13.121151733).jd, + "Reactivation_2024": MissionPhase( + name="Reactivation_2024", + jd_start=Time.from_ymd(2024, 1, 1).jd, + jd_end=Time.from_ymd(2024, 8, 1.291525).jd, bands=(1, 2), frame_url=IRSA_URL + "/ibe/data/wise/neowiser/p1bm_frm/", frame_meta_table="neowiser_p1bs_frm", @@ -224,7 +217,7 @@ } """Public released mission phases of WISE.""" -for year in range(2015, 2023): +for year in range(2015, 2024): MISSION_PHASES[f"Reactivation_{year}"] = MissionPhase( name=f"Reactivation_{year}", jd_start=Time.from_ymd(year, 1, 1).jd, @@ -340,6 +333,10 @@ def mission_phase_from_scan(scan_id: str) -> Optional[MissionPhase]: return MISSION_PHASES["Reactivation_2022"] elif scan_num <= 57041: return MISSION_PHASES["Reactivation_2023"] + elif scan_num <= 57626: + return MISSION_PHASES["Reactivation_2023"] + elif scan_num <= 64272: + return MISSION_PHASES["Reactivation_2024"] return None elif letter == "s": if scan_num <= 1615: @@ -352,8 +349,10 @@ def mission_phase_from_scan(scan_id: str) -> Optional[MissionPhase]: return MISSION_PHASES["Reactivation_2021"] elif scan_num <= 46369: return MISSION_PHASES["Reactivation_2022"] - elif scan_num <= 56807: + elif scan_num <= 57519: return MISSION_PHASES["Reactivation_2023"] + elif scan_num <= 64267: + return MISSION_PHASES["Reactivation_2024"] return None return None @@ -376,7 +375,9 @@ def _scan_frame(scan_id, frame_num=None, band=None): return scan_id, frame_num, band -def fetch_frame(scan_id, frame_num=None, band=None, as_fits=True, im_type="int"): +def fetch_frame( + scan_id, frame_num=None, band=None, as_fits=True, im_type="int", retry=2 +): """ Fetch the WISE FITs frame, if it is not present in the cache, download it first. @@ -413,49 +414,27 @@ def fetch_frame(scan_id, frame_num=None, band=None, as_fits=True, im_type="int") # format the url frame_num = f"{frame_num:03d}" scan_group = str(scan_id[-2:]) - band = f"w{band:1d}" + band_str = f"w{band:1d}" ext = "fits" if im_type == "int" else "fits.gz" - filename = f"{scan_id}{frame_num}-{band}-{im_type}-1b.{ext}" + filename = f"{scan_id}{frame_num}-{band_str}-{im_type}-1b.{ext}" url = f"{phase.frame_url}{scan_group}/{scan_id}/{frame_num}/{filename}" subfolder = os.path.join("wise_frames", scan_group) file_path = download_file(url, auto_zip=True, subfolder=subfolder) if as_fits: - return fits.open(file_path)[0] + try: + return fits.open(file_path)[0] + except OSError as exc: + if retry == 0: + raise ValueError("Failed to fetch WISE frame.") from exc + logger.info("WISE file appears corrupted, attempting to fetch again.") + os.remove(file_path) + return fetch_frame(scan_id, frame_num, band, as_fits, im_type, retry - 1) else: return file_path -_frame_path = cache_path("wise_frames") -_reorg_msg = ( - " The organization of wise cached frames has changed, it is recommended" - " to delete the old cache files with the shell command: \n" - f"rm {_frame_path}/*.fits" - "\nThis only has to be done once. Frames are now stored in subfolders " - "based on the scan group number. This allows for many more files to be" - " saved locally before there are filesystem issues." -) -cache_WISE_frame = partial( - rename( - fetch_frame, - additional_msg=( - "Use `as_fits=False` with as an argument in the new function." + _reorg_msg - ), - deprecated_version="v0.2.4", - old_name="cache_WISE_frame", - ), - as_fits=False, -) - -fetch_WISE_frame = rename( - fetch_frame, - additional_msg=_reorg_msg, - deprecated_version="v0.2.4", - old_name="fetch_WISE_frame", -) - - def plot_frames( scan_id, frame_num=None, @@ -547,7 +526,7 @@ def plot_frames( plt.title(f"W{band:1d}") -@lru_cache() +@lru_cache(maxsize=2) def fetch_WISE_fovs(phase): """ Load all FOVs taken during the specified mission phase of WISE. diff --git a/src/kete/ztf.py b/src/kete/ztf.py index 0040ca8..c5fb1b8 100644 --- a/src/kete/ztf.py +++ b/src/kete/ztf.py @@ -2,20 +2,11 @@ ZTF Related Functions and Data. """ -from functools import lru_cache import os +import logging +from functools import lru_cache from collections import defaultdict -import matplotlib.pyplot as plt -import warnings - from astropy.io import fits -from astropy.wcs import WCS -from astropy.visualization import ( - PowerDistStretch, - AsymmetricPercentileInterval, - ImageNormalize, -) - from .cache import download_file, cache_path from .fov import ZtfCcdQuad, ZtfField, FOVList @@ -26,31 +17,15 @@ from . import spice -__all__ = ["fetch_ZTF_file", "fetch_ZTF_fovs"] +__all__ = ["fetch_ZTF_file", "fetch_ZTF_fovs", "fetch_frame"] SURVEY_START_JD = Time.from_ymd(2018, 3, 20).jd """First image in ZTF dataset.""" -ZTF_IRSA_TABLES = { - "ztf.ztf_current_meta_cal": "ZTF Calibration Metadata Table", - "ztf.ztf_current_meta_deep": "ZTF Deep Reference Images", - "ztf.ztf_current_meta_raw": "ZTF Raw Metadata Table", - "ztf.ztf_current_meta_ref": "ZTF Reference (coadd) Images", - "ztf.ztf_current_meta_sci": "ZTF Science Exposure Images", - "ztf.ztf_current_path_cal": "ZTF Calibration Product Paths", - "ztf.ztf_current_path_deep": "ZTF Deep Reference Product Paths", - "ztf.ztf_current_path_raw": "ZTF Raw Product Paths", - "ztf.ztf_current_path_ref": "ZTF Reference Product Paths", - "ztf.ztf_current_path_sci": "ZTF Science Product Paths", - "ztf_objects_dr16": "ZTF Data Release 16 Objects", - "ztf_objects_dr17": "ZTF Data Release 17 Objects", - "ztf_objects_dr18": "ZTF Data Release 18 Objects", - "ztf_objects_dr19": "ZTF Data Release 19 Objects", - "ztf_objects_dr20": "ZTF Data Release 20 Objects", -} - - -@lru_cache(maxsize=3) +logger = logging.getLogger(__name__) + + +@lru_cache(maxsize=2) def fetch_ZTF_fovs(year: int): """ Load all FOVs taken during the specified mission year of ZTF. @@ -65,7 +40,7 @@ def fetch_ZTF_fovs(year: int): Which year of ZTF, 2018 through 2024. """ year = int(year) - if year not in [2018, 2019, 2020, 2021, 2022, 2023, 2024]: + if year not in range(2018, 2025): raise ValueError("Year must only be in the range 2018-2024") cache_dir = cache_path() dir_path = os.path.join(cache_dir, "fovs") @@ -172,63 +147,12 @@ def file_frac_day_split(filefracday): return (year, month, day, frac_day) -def plot_frame( - fov: ZtfCcdQuad, - cmap="grey", - products="sci", - im_type="sciimg.fits", - force_download=False, -): - """ - Given a ztf FOV, plot the associated frame. - - This returns the associate WCS which is constructed. - - Parameters - ---------- - fov : - A single CCD Quad FOV. - cmap : - Colormap of the plot. - products : - Which data product to fetch. - im_type : - Image extension, this must match the products variable. - force_download : - Optionally force a re-download if the file already exists in the cache. - """ - - frame = fetch_frame_from_fov( - fov, products=products, im_type=im_type, force_download=force_download - ) - - with warnings.catch_warnings(): - warnings.filterwarnings("ignore") - wcs = WCS(frame.header) - - if not plt.get_fignums(): - plt.figure(dpi=300, figsize=(6, 6), facecolor="w") - - ax = plt.subplot(projection=wcs) - - norm = ImageNormalize( - frame.data, - interval=AsymmetricPercentileInterval(10, 99.5), - stretch=PowerDistStretch(0.25), - ) - - ax.imshow(frame.data, origin="lower", norm=norm, cmap=cmap) - ax.set_xlabel("RA") - ax.set_ylabel("DEC") - ax.set_aspect("equal", adjustable="box") - return wcs - - -def fetch_frame_from_fov( +def fetch_frame( fov: ZtfCcdQuad, products="sci", im_type="sciimg.fits", force_download=False, + retry=2, ) -> fits.hdu.image.PrimaryHDU: """ Given a ztf FOV, return the FITs file associated with it. @@ -257,7 +181,14 @@ def fetch_frame_from_fov( im_type, force_download, ) - return fits.open(file)[0] + try: + return fits.open(file)[0] + except OSError as exc: + if retry == 0: + raise ValueError("Failed to fetch ZTF frame.") from exc + logger.info("ZTF file appears corrupted, attempting to fetch again.") + os.remove(file) + return fetch_frame(fov, products, im_type, force_download, retry - 1) def fetch_ZTF_file( diff --git a/src/kete_core/src/constants/gravity.rs b/src/kete_core/src/constants/gravity.rs index 99920e6..e84fb5f 100644 --- a/src/kete_core/src/constants/gravity.rs +++ b/src/kete_core/src/constants/gravity.rs @@ -10,7 +10,7 @@ pub const GMS: f64 = 0.00029591220828411956; /// Gaussian gravitational constant, equivalent to sqrt of GMS. /// AU^(3/2) per (Day sqrt(Solar Mass)) -pub const GMS_SQRT: f64 = 0.01720209894846; +pub const GMS_SQRT: f64 = 0.01720209894996; /// Sun J2 Parameter /// diff --git a/src/kete_core/src/elements.rs b/src/kete_core/src/elements.rs index cf1445f..14eaa08 100644 --- a/src/kete_core/src/elements.rs +++ b/src/kete_core/src/elements.rs @@ -214,32 +214,32 @@ impl CometElements { if elliptical { let (sin_e, cos_e) = ecc_anom.sin_cos(); - let e_dot = (semi_major.powf(1.5) * (1.0 - self.eccentricity * cos_e)).recip(); + let e_dot = semi_major.powf(1.5) * (1.0 - self.eccentricity * cos_e); let b = semi_major * (1.0 - self.eccentricity.powi(2)).sqrt(); x = semi_major * (cos_e - self.eccentricity); y = b * sin_e; - x_dot = -semi_major * e_dot * sin_e * GMS_SQRT; - y_dot = b * e_dot * cos_e * GMS_SQRT; + x_dot = -semi_major / e_dot * sin_e * GMS_SQRT; + y_dot = b / e_dot * cos_e * GMS_SQRT; } else if hyperbolic { let sinh_h = ecc_anom.sinh(); let cosh_h = ecc_anom.cosh(); let b = -semi_major * (self.eccentricity.powi(2) - 1.0).sqrt(); - let h_dot = (semi_major.abs().powf(1.5) * (1.0 - self.eccentricity * cosh_h)).recip(); + let h_dot = semi_major.abs().powf(1.5) * (1.0 - self.eccentricity * cosh_h); x = semi_major * (cosh_h - self.eccentricity); y = b * sinh_h; - x_dot = -semi_major * h_dot * sinh_h * GMS_SQRT; - y_dot = -b * h_dot * cosh_h * GMS_SQRT; + x_dot = -semi_major / h_dot * sinh_h * GMS_SQRT; + y_dot = -b / h_dot * cosh_h * GMS_SQRT; } else { // Parabolic - let d_dot = (self.peri_dist + ecc_anom.powi(2) / 2.0).recip(); + let d_dot = self.peri_dist + ecc_anom.powi(2) / 2.0; x = self.peri_dist - ecc_anom.powi(2) / 2.0; y = (2.0 * self.peri_dist).sqrt() * ecc_anom; - x_dot = -d_dot * ecc_anom * GMS_SQRT; - y_dot = d_dot * (2.0 * self.peri_dist).sqrt() * GMS_SQRT + x_dot = -ecc_anom / d_dot * GMS_SQRT; + y_dot = (2.0 * self.peri_dist).sqrt() / d_dot * GMS_SQRT; } let (s_w, c_w) = self.peri_arg.sin_cos(); diff --git a/src/tests/test_wise.py b/src/tests/test_wise.py index d1d95b5..7ec5da7 100644 --- a/src/tests/test_wise.py +++ b/src/tests/test_wise.py @@ -19,7 +19,7 @@ def test_mission_phase_from_scan(): for letter in "rs": scan_id = "10000" + letter assert mission_phase_from_scan(scan_id) == MISSION_PHASES["Reactivation_2019"] - scan_id = "57042" + letter + scan_id = "64273" + letter assert mission_phase_from_scan(scan_id) is None scan_id = "01000a"